Monday, September 20, 2010

R: Log-Linear Modelling

# In general, I didn't code this myself.
# Compiled from various sources.

# Data input: generic

# Red indicates stuff you have to edit.

# - Start of Log-linear Modelling Script 1 -

# lib(s): stats

# before using this script:
library(stats)

mytab <- table(Cat.Var1,Cat.Var2)# where the categories of Cat.Var1 will be the rows and Cat.Var2 will be the columns.
loglin(mytab,margin=c(1,2)) # Produces a seemingly coherent output. But i no idea in hell what's going on. WTF is margin???

# description of margin: a list of vectors with the marginal totals to be fit. (Hierarchical) log-linear models can be specified in terms of these marginal totals which give the ‘maximal’ factor subsets contained in the model. For example, in a three-factor model, list(c(1, 2), c(1, 3)) specifies a model which contains parameters for the grand mean, each factor, and the 1-2 and 1-3 interactions, respectively (but no 2-3 or 1-2-3 interaction), i.e., a model where factors 2 and 3 are independent conditional on factor 1 (sometimes represented as ‘[12][13]’).

# - End of Log-linear Modelling Script 1 -

# - Start of Log-linear Modelling Script 2 - alpha tested ok on 14 Sept 2010

# lib(s): MASS

# before using this script:
library(MASS)

loglm(~ Cat.Var1 + Cat.Var2, xtabs(~ Cat.Var1 + Cat.Var2, infile))

# - End of Log-linear Modelling Script 2 -

# - Start of Log-linear Modelling Script 3 - alpha tested ok on 14 Sept 2010.
# (This is a somewhat DIY log-linear modelling script. All it reqires is the input of a table like the one mentioned at the start of the log-linear modelling section.)

# Log-Linear analysis using General Linear Model (http://users.ox.ac.uk/~polf0050/ISS%20Lecture%206.pdf(?))

mytab <- table(infile$CarVar1,infile$CatVar2)

lg.lnr <- function(mytab)
{
r.str <- rep(rownames(mytab),times=length(colnames(mytab)))
c.str <- rep(colnames(mytab),each=length(rownames(mytab)))
freq <- as.vector(mytab)
dat.frm <- data.frame(r.str,c.str,freq)
result <- glm(freq~r.str+c.str,family="poisson",data=dat.frm)
r1 <- anova.glm(result,test="Chisq")
r2 <- anova(result,test="Chisq")
plot(result)
r3 <- summary(result)
to.print <- list(r1,r2,r3)
to.print
}

# - End of Log-linear Modelling Script 3 -

# /* END OF LOG-LINEAR MODELLING */

No comments:

Post a Comment