How is a covariate added to a log-linear model?

There are a number of good Q/As about log-linear models (i.e., here). The data I am analyzing is off shifts between categories, 0, 1, 2, 3, 4, and 5, within students, in that Var1 represents a student's category at Time 1, and Var2 represents a student's category at Time 2.

I'm trying to use a log-linear model to analyze relations between two variables in the following frequency table, tab:

> tab
Var2
Var1  0  1  2  3  4  5
0  2  5  6  2  9  0
1  2 18 13  4 11  1
2  5 11 34  4 16  0
3  1  4 11  1  5  0
4  1  8 12  8 13  0
5  0  0  1  0  0  0


I'm using the following R syntax:

m1 <- MASS::loglm( ~ Var1 + Var2, tab)


This prints the Likelihood Ratio and Pearson statistics:

summary(m1)


And this obtains the standardized residuals:

as.data.frame(resid(m1))


My question is a bit conceptual and a bit technical: How is a covariate, such as a students' baseline level of some other variable that may influence both Var1 and Var2, added to the model? The reason I ask: I'm exploring whether it is possible to fit log-linear models using a multi-level model, in order to account for student dependencies in the data.

Code for tab is here, if it were helpful at all:

tab <- structure(c(2, 2, 5, 1, 1, 0, 5, 18, 11, 4, 8, 0, 6, 13, 34,
11, 12, 1, 2, 4, 4, 1, 8, 0, 9, 11, 16, 5, 13, 0, 0, 1, 0, 0,
0, 0), .Dim = c(6L, 6L), .Dimnames = structure(list(Var1 = c("0",
"1", "2", "3", "4", "5"), Var2 = c("0", "1", "2", "3", "4", "5"
)), .Names = c("Var1", "Var2")), class = c("xtabs", "table"), call = xtabs(formula = Freq ~ Var1 + Var2, data = x))


Here is the raw data from which tab was constructed using tab <- table(df$Var1, df$Var2):

df <- structure(list(Var1 = c(3, 4, 0, 1, 2, 3, 2, 2, 1, 2, 2, 1, 4,
1, 2, 1, 1, 1, 2, 0, 2, 0, 2, 2, 3, 2, 1, 3, 2, 2, 4, 2, 3, 1,
4, 2, 4, 2, 2, 2, 0, 3, 1, 2, 4, 2, 1, 2, 2, 2, 2, 1, 2, 1, 1,
4, 2, 3, 2, 4, 1, 4, 1, 1, 1, 4, 1, 2, 2, 2, 2, 1, 2, 1, 0, 3,
2, 1, 2, 3, 0, 0, 0, 2, 0, 0, 2, 1, 2, 2, 1, 1, 3, 1, 1, 2, 1,
3, 2, 4, 4, 2, 1, 3, 4, 1, 1, 0, 2, 2, 4, 1, 2, 3, 1, 1, 3, 2,
3, 1, 2, 4, 1, 4, 1, 3, 2, 4, 2, 3, 4, 2, 4, 4, 3, 2, 2, 1, 3,
2, 3, 2, 2, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 4, 2, 1, 1, 1, 1, 1,
1, 1, 1, 4, 1, 1, 4, 1, 4, 1, 4, 1, 1, 4, 1, 1, 2, 2, 1, 4, 2,
3, 2, 4, 2, 4, 2, 1, 1, 1, 1, 2, 4, 4, 1, 4, 1, 2, 4, 4, 4, 1,
4, 4, 4, 4, 4, 1, 4, 3, 4, 4, 1, 1, 4, 4, 1, 1, 2, 4, 2, 1, 4,
1, 2, 3, 4, 2, 1, 1, 4, 2, 2, 1, 1, 1, 2, 1, 1, 4, 1, 2, 4, 3,
1, 0, 2, 2, 4, 2, 2, 2, 0, 2, 4, 2, 4, 2, 1, 2, 0, 4, 4, 4, 2,
2, 1, 3, 2, 4, 4, 1, 4, 2, 2, 2, 2, 1, 4, 2), Var2 = c(1, 1,
3, 1, 1, 4, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 2, 2, 2, 0, 3, 1, 1,
3, 1, 2, 1, 3, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 4, 3, 1, 3, 1, 0,
2, 2, 2, 1, 1, 4, 1, 3, 1, 1, 2, 1, 0, 2, 2, 4, 1, 1, 1, 4, 1,
1, 0, 0, 1, 2, 1, 1, 3, 3, 0, 2, 3, 0, 2, 1, 0, 2, 2, 2, 2, 1,
2, 4, 4, 4, 3, 4, 3, 4, 4, 2, 2, 2, 4, 2, 4, 2, 2, 0, 4, 2, 1,
0, 2, 2, 2, 4, 0, 3, 4, 4, 2, 2, 4, 3, 2, 4, 5, 4, 1, 2, 0, 4,
2, 3, 4, 2, 4, 3, 3, 2, 2, 2, 4, 0, 2, 3, 2, 0, 0, 2, 4, 2, 4,
0, 2, 3, 2, 2, 2, 4, 0, 0, 2, 0, 0, 4, 1, 0, 0, 4, 0, 2, 4, 0,
0, 1, 1, 3, 2, 1, 0, 2, 0, 4, 0, 2, 0, 2, 0, 3, 4, 4, 1, 1, 2,
1, 2, 4, 3, 1, 2, 1, 2, 0, 4, 4, 4, 4, 4, 0, 4, 1, 1, 0, 1, 2,
4, 1, 1, 4, 0, 1, 2, 4, 2, 2, 2, 1, 2, 4, 1, 3, 1, 4, 2, 2, 4,
1, 1, 1, 4, 4, 4, 4, 1, 4, 1, 1, 1, 4, 4, 2, 0, 2, 2, 4, 1, 4,
2, 2, 2, 2, 2, 4, 2, 4, 2, 4, 4, 4, 2, 4, 2, 2, 2, 2, 4, 2, 3,
4, 4, 2, 4, 0)), class = "data.frame", .Names = c("Var1", "Var2"
), row.names = c(NA, -280L))

• Do you not have the data per student? (e.g. student 1 started in category 0 and moved to category 3, student 2 started in category 2 and remained in 2, etc...) – IWS Jun 19 '17 at 12:15
• If you are not specifically looking to use a log-linear model, I would model this using a multinomial regression model, using var1 as predictor and var2 as outcome. If there is enough data to fit this model this would give you the expected probability for each var1 category to end up in a specific var2 category. – IWS Jun 19 '17 at 12:23
• I imagine it would be straightforward to add covariates to that model. Any suggestions on whether it would be possible to integrate the hierarchical structure (i.e., random effects) through such an approach? – Joshua Rosenberg Jun 19 '17 at 12:48