Skip to content

Question about misspecified multilevel model

3 messages · Phillip Alday, Christopher David Desjardins

Original
Christopher David Desjardins
# 
Hi,

I have a question that's potentially off-topic but that I'm hoping that
someone here can shed some insight on.

Assume that I know that I know my true model and that my true is a
three-level model. My observations are such that I have a measurement on a
student nested within a classroom nested within a school. The true model
would be:

Y_ijk = pi_0jk + e_ijk  # student within classroom within schools (1st
level)

pi_0jk = beta_j0k + r_p0k  # classroom within schools (2nd level)

beta_j0k = gamma_pq0 + u_pqk  # schools (3rd level)

The model in lmer would be:
classroom <- read.csv("http://www-personal.umich.edu/~bwest/classroom.csv")
library("lme4")
correct.mod <- lmer(mathgain ~ (1 | schoolid/classid), data = classroom)

What I am wondering about is, if I were to omit that second level, the
whole classroom within schools equation, where would that variance that
would end up as the random intercept go? Would it go to the random
intercept for school or would it go down to the residual?  The model I am
referring to is below:

misspecified.mod <- lmer(mathgain ~ (1 | schoolid), data = classroom)
summary(correct.mod); summary(misspecified.mod)

It looks like the variances for both the residual and the random intercept
for school change. But maybe they do in a predictable way?

If someone could suggest a paper has an answer or better that would be very
helpful.

Chris
Phillip Alday
# 
Hi Chris,

Since you have a concrete example why not just calculate the models and see what happens? The intercept-only model that you proposed is quite simple and fast to compute:

Random effects:
 Groups           Name        Variance Std.Dev.
 classid:schoolid (Intercept)   99.23   9.961  
 schoolid         (Intercept)   77.49   8.803  
 Residual                     1028.23  32.066  
Number of obs: 1190, groups:  classid:schoolid, 312; schoolid, 107

Now, if you do a model with just classroom, then you want this model:

class.mod <- lmer(mathgain ~ (1 | schoolid:classid), data = classroom)

(I'm assuming that classid isn't globally unique, but rather only unique within schools, otherwise the / and : operators add a bit of unnecessary overhead, but that's a different issue)

which has the following random effects: 

Random effects:
 Groups           Name        Variance Std.Dev.
 schoolid:classid (Intercept)  180     13.42   
 Residual                     1025     32.02   
Number of obs: 1190, groups:  schoolid:classid, 312

As you can see, the classid absorbed most of the variance. If you look at the fixed effects, then you'll see slight differences in the error estimates there as well -- models with incomplete random-effect structures typically have error estimates that are two small and are thus anti-conversative. This can be attributed to the variance-bias tradeoff or equivalently overfitting (see e.g. the new book Statistical Rethinking or any of the other classics mentioned here such as Pinheiro & Bates or Gelman & Hill). Roughly, leaving out random effects structure implies that the data are independent in in a way that they are not, which leads you to overestimate how much information you actually have and thus how little error.

Now, if you're curious you can also compare log likelihood if you refit the models via ML*:
refitting model(s) with ML (instead of REML)
Data: classroom
Models:
class.mod: mathgain ~ (1 | schoolid:classid)
correct.mod: mathgain ~ (1 | schoolid/classid)
            Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
class.mod    3 11785 11800 -5889.5    11779                            
correct.mod  4 11779 11800 -5885.7    11771 7.5906      1   0.005867 **
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

And you see that the model fits while tantalisingly close are nonetheless significantly different. 

Now, I'm not sure what will happen with more complex fixed-effects or random-effects structures, but my guess is that generally the variance will be allocated to whatever variable is "close enough" (here classroom is close enough to school because a collection of classrooms makes up a school) and if there are none close enough, the remaining variance will just contribute to the residual variance (which didn't happen here).

* The fixed-effects model matrix is the same for both models so that you could potentially compare REML-fitted models, but that also gets you into all sorts of fun debate about the pros and cons of REML and it's just easier to avoid all the issues with comparisons of REML-fit

Best,
Phillip
Christopher David Desjardins
# 
Thanks, Phillip. The concrete example was meant to be just an example, so
you would know what I am talking about. I am thinking about
misspecification in general. Thanks for the detailed response to that. I am
cutting your response below.

On Saturday, August 6, 2016, Phillip Alday <Phillip.Alday at unisa.edu.au>
wrote:
This is what I am wondering about out and this is what I expected happens
but wasn't sure if it was just for this example and wasn't sure if there
was a mathematical reason why.