Solved – Specifying group covariates and cross level interactions in lmer

interactionlme4-nlmemixed modelr

I'm new to R, and despite trying to read as much as I can about how lmer works in R, I still don't feel like I know how to correctly specify more complex models using the lmer syntax.

For example, at the moment I want to use lmer for a two-level multilevel model, where the first level is features of a specific course taken by a specific student in a specific semester (with covariates like kind of course, teaching method, etc) and the second level is the student, which also has a number of covariates (e.g. ethnicity, gender, age, gpa at the beginning of the study, and score on a specific instrument). I also want to assess interactions with the teaching method, for example by including cross-level interactions between the second-level covariates ethnicity, gender, age, gpa, and the instrument score and the first-level covariate teaching method.

For the sake of readability, let's limit the equation to teaching method, course type, and score. It seems to me that I have seen three different ways of doing something like what I want to do in R – they are clearly all specifying something somewhat different, but I can't figure out what the underlying math is supposed to be for each case. So, for example, in different online references and the books on MLM that I have on hand, it seems that one of these three models is recommended for what I want to do:

course_outcome ~ course_type*teaching_method + score*teaching_method + (1|student)
course_outcome ~ course_type*score + (teaching_method|student)
course_outcome ~ course_type*teaching_method + score*teaching_method + (teaching_method|student)

I'm a little confused about the differences between how lmer interprets each of these three codings – is there any chance that someone who really understands R better could possibly translate this into the basic regression equation(s) structure that would be calculated in each of these three examples? Or direct me to a reference that might more clearly explain the difference by giving the regression equations for each case?

Thanks for your time!

Best Answer

The transformations from lmer codes to regression equations are listed below, where $u$ denotes the random effect, $i$ and $j$ indicate student and course respectively. Suppose the course features would not change with the students.

course_outcome ~ course_type*teaching_method + score*teaching_method + (1|student)

$$\mathrm{course\_outcome}_{ij}= \beta_0 +\beta_1\mathrm{course\_type}_j +\beta_2\mathrm{teaching\_method}_j +\beta_3\mathrm{score}_i +\beta_4\mathrm{course\_type}_j*\mathrm{teaching\_method}_j +\beta_5\mathrm{score}_i*\mathrm{teaching\_method}_j +u_i+e_{ij}$$
course_outcome ~ course_type*score + (teaching_method|student)

$$\mathrm{course\_outcome}_{ij}= \beta_0 +\beta_1\mathrm{course\_type}_j +\beta_2\mathrm{score}_i +\beta_3\mathrm{course\_type}_j*\mathrm{score}_i +u_{0i} +u_{1i}\mathrm{teaching\_method}_j +e_{ij}$$
course_outcome ~ course_type*teaching_method + score*teaching_method + (teaching_method|student)

$$\mathrm{course\_outcome}_{ij}= \beta_0 +\beta_1\mathrm{course\_type}_j +\beta_2\mathrm{teaching\_method}_j +\beta_3\mathrm{score}_i +\beta_4\mathrm{course\_type}_j*\mathrm{teaching\_method}_j +\beta_5\mathrm{score}_i*\mathrm{teaching\_method}_j +u_{0i}+u_{1i}\mathrm{teaching\_method}_j +e_{ij}$$

It may help you form some general ideas to interpret the codes. Note that $\beta$s across the equations have different meanings. It seems that course_type and teaching_method may be factors. If so, they would have several coefficients corresponding to the number of levels.

Related Solutions

Solved – Am I specifying the lmer model correctly

Your model specification is fine.

The varying intercept for Ward, specified in lmer as you've done with (1 | Ward), is saying that subjects within each ward might be more similar to each other on Selfreject for reasons other than WardSize or Gender, so you are controlling for between-ward heterogeneity.

You can think of the "1" as a column of 1s (i.e., a constant) in the data to which an intercept is fit. Usually the "1" is implied automatically in lm, for instance

lm(Y ~ X1 + X2)

actually specifies

lm(Y ~ 1 + X1 + X2)

Now that you have your basic model, you can start asking further questions like "Does the relationship between BSItotal and Selfreject differ between wards?"

lmer(formula=Selfreject ~ WardSize + WAS + Gender + BSITotal + (1 + BSITotal | Ward))

That is, both the intercept and the slope of BSITotal can differ between wards.

If you haven't picked it up yet, Gelman & Hill's Data Analysis Using Regression and Multilevel Model/Hierarchical Models is a great book that explains fitting models like this with lmer.

Solved – Factors missing in summary and Anova output from glm in R

There are three common causes for the zero degrees of freedom:

Saturated model. These are usually spotted from NAs is summary(model))
Factor level combinations with zero observations.These are spotted as zero dfs in Anova(model) for the interaction effects, but the main effect would still have, e.g., one df.
Perfectly collinear explanatory variables. These are spotted from the zero dfs for the interaction effects in presense of the main effects in Anova(model), or if only main effects are included the model, some of the main effect could have zero dfs.

Your case falls under the third cause. Judging from the specification of the models 1 and 2, some of your explanatory variables you have added to model 2 (ethnicity, gender, age, enrollment, Pell, TANF, GPA, onlineExp) are collinear with level and type.

For finding out the problematic variables, consult the answer on Stackoverflow. Briefly, command model.matrix(model) gives out, well, the model matrix for your model. Check the output! Following the advice George Dontas gave on the linked thread, you can then find out which ones are the problematic predictors. Small eigen values from the following command indicate problems:

eigen(cov(model.matrix(model)[,-1]))$values

For a very small R example, try the following:

set.seed(1)
y<-rnorm(9)
x1<-c(0,0,0,1,1,1,2,2,2)
x2<-c(0,1,2,0,1,2,0,1,2)
x3<-c(1,1,1,2,2,2,3,3,3)
fit1<-lm(y~x1+x2)
fit2<-lm(y~x1+x2+x3)
fit3<-lm(y~x1+x2+x3+x1*x3)
library(car)
Anova(fit1)
Anova(fit2) # main effects only
Anova(fit3) # interactions, also
e<-eigen(cov(model.matrix(fit2)[,-1]))$values
names(e)<-colnames(model.matrix(fit2)[,-1]) # x3 creates the problem, the value is very low

Another possibility is to add the new variables one by one, and see when the problem shows itself. You can then try to rectify the situation, e.g., by not adding the problematic variables to the model. Sometimes combining the variables in a suitable way helps, also.

Best Answer

Related Solutions

Solved – Am I specifying the lmer model correctly

Solved – Factors missing in summary and Anova output from glm in R

Related Question