# How do I setup a model with hierarchical structure using lmer in R?

I am trying to isolate the important predictors for my response variable "Y". I know that "TL" (which is an individual level predictor) affects "Y", and now I want to determine if adding the site level predictors "PC1" and "PC4" (both continuous variables) improve the fit of the model.

The data looks like this:

Y    TL    PC1    PC4    Site
5.6  17    -2.26  1.89   Site A
5.9  19    -2.26  1.89   Site A
5.2  20    -2.26  1.89   Site A
4.0  17    -1.56  2.34   Site B
5.0  19    -1.56  2.34   Site B
5.6  18    -1.56  2.34   Site B


Overall I have approx. 15 samples at thirteen sites. Here are three models I setup:

model1<-lmer(Y~ TL + (1 + TL | Site), data=sauro.all) #(no site level predictors)
model2<-lmer(Y ~ TL + PC1 + (1 + TL | Site), data=sauro.all) #PC1 as site predictor)
model3<-lmer(Y ~ TL + PC4 + (1 + TL | Site), data=sauro.all) #PC4 as site level predictor


Summary of Model 1:

Linear mixed model fit by maximum likelihood
Formula: Y ~ TL + (1 + TL | Site)
Data: sauro.all
AIC   BIC logLik deviance REMLdev
400.5 419.6 -194.3    388.5   395.4
Random effects:
Groups   Name        Variance   Std.Dev. Corr
Site     (Intercept) 0.24335394 0.493309
TL          0.00043485 0.020853 1.000
Residual             0.40595113 0.637143
Number of obs: 178, groups: Site, 13

Fixed effects:
Estimate Std. Error t value
(Intercept) -16.84997    0.41581  -40.52
TL            0.10890    0.02012    5.41

Correlation of Fixed Effects:
(Intr)
TL -0.803


Summary of Model 2:

Linear mixed model fit by maximum likelihood
Formula: Y ~ TL + PC1 + (1 | Site)
Data: sauro.all
AIC   BIC logLik deviance REMLdev
396.3 412.2 -193.1    386.3   395.1
Random effects:
Groups   Name        Variance Std.Dev.
Site     (Intercept) 0.67023  0.81868
Residual             0.40740  0.63828
Number of obs: 178, groups: Site, 13

Fixed effects:
Estimate Std. Error t value
(Intercept) -16.89643    0.45413  -37.21
TL            0.11022    0.01927    5.72
PC1          -0.32754    0.17083   -1.92

Correlation of Fixed Effects:
(Intr) TL
TL  -0.859
PC1  0.054 -0.054


...and summary of model3:

Linear mixed model fit by maximum likelihood
Formula: Y ~ TL + PC4 + (1 + TL | Site)
Data: sauro.all
AIC   BIC logLik deviance REMLdev
399 421.2 -192.5      385   393.3
Random effects:
Groups   Name        Variance   Std.Dev. Corr
Site     (Intercept) 0.13411114 0.366212
TL          0.00044998 0.021213 1.000
Residual             0.40582841 0.637047
Number of obs: 178, groups: Site, 13

Fixed effects:
Estimate Std. Error t value
(Intercept) -16.83238    0.40527  -41.53
TL            0.11003    0.02016    5.46
PC4          -0.45466    0.22571   -2.01

Correlation of Fixed Effects:
(Intr) TL
TL  -0.845
PC4  0.005 -0.052


Model 2 has the lowest AIC score so I think this is the best model. I am having trouble with interpretation. My questions are: *How do I interpret the contribution of PC1 to model2 without p values? *This analysis is telling me how PC1 affects intercepts of (Y~TL) among sites, correct? However, slopes of this relationship could also vary with PC1. How do I check this? *What is the best way to visualize an lmer model?

Thank you!

It makes perfect sense for model 2 to be better. I didn't even have to look at the AIC. Look at the variance of your random effects. model1 and model3 are barely hanging to have impact as high as your residuals variance (measurement error if you like) and model2 has higher variance. It is a no-brainer. Also check out the correlation between your intercept and $TL$, 1 ? That's plain wrong, you definitely need to remove something (ie. give more freedom to your model). At best your design caused lme4 not to optimize properly.

More specifically:

How do I interpret the contribution of PC1 to model2 without p values? you don't need $p$-values to tell you a model is correct. Interpreter it as it is. And if someone is really picky on that matter (so doesn't really understand what $p$-value is in the context of LME.) just feed him the result of an ANOVA between the two models to get a $p$-value (I wouldn't believe it be true but hey, some guys just love their $p$-values). I would suggest just to say your model came out as the best model based on an $AIC$ selection criterion. (Because really that's why you just said you looked at in the first place and it is a very reasonable to do indeed in a lot of cases.) Just take notice that if you have difference in the AIC scores of the two models that is smaller than 2 you can't reject a model instead of the other. I don't know if it still works but you could look into MCMC generated $p$-values, using mcmcsamp, those are pretty much the golden standard (if the chain mixes well that is). I try to report those if ever possible.

This analysis is telling me how PC1 affects intercepts of (Y~TL) among sites, correct? Correct.

However, slopes of this relationship could also vary with PC1. How do I check this? Use interaction term as ndoogan says. (*)

For "standard visualization" you can check the following thread: What would be an illustrative picture for linear mixed models? but as ndoogan says, if you are slightly more specific you'll get better answers when it comes to your visualization.

(*) I just realized you have a problem. I didn't read your initial model specification R-code but now that I do I see that the summary of Model 2 you are presenting is not the summary of model2 you define. You define model2<-lmer(Y ~ TL + PC1 + (1 + TL | Site), data=sauro.all) #PC1 as site predictor) but your fitted model is Y ~ TL + PC1 + (1 | Site) so you have eaten up your TL slope. Are you sure you don't have a mistake? Both answers you got so far are perfectly valid but it really looks like you are not comparing models with the same random effects structure. (I actually think that your "accidental" structure is the correct one)

• Yes, that was a mistake in the fitted model for model2. However, in comparing the model allowing slopes to vary and one without, it looks like a better fit to not let slopes vary. – Laura May 2 '13 at 20:06
• There you have it then; do not allow for different slopes. If, based on modelling assumptions and your metric, "no slopes" give you a "better fit" that is the correct thing to do. No reason to be shy about it but importantly see what kind of limitations your theoretical choices entail so you can appreciate what you can say and what you can't. – usεr11852 says Reinstate Monic May 3 '13 at 0:43
• I forgot to add that you should compare linear mixed models fitted using ML and not ReML. (REML=F in R) You appear to use the default option (ReML=T). – usεr11852 says Reinstate Monic May 3 '13 at 0:47

Preface

There is a discrepancy in your output. You do not allow $TL$ slope to vary according to the output for model 2, which could explain the difference in AIC if that additional error term doesn't account for much of the variation in the data.

Title

To answer the question in your title (which appears to be somewhat of a mismatch to the rest of what you wrote), I think you've done a fine job of setting the model up properly for hierarchically structured data. I say this because you've included two error terms to capture random variation in the intercept and the $TL$ slope across groups defined by the data collection site.

Other

Model 2 tells you that $PC1$ is linearly and negatively predictive of $Y$. To assess statistical significance, you could take a look at the $t$-value supplied by lmer. Very basically, $p$ values are not given because it's unclear in many cases what the degrees of freedom should be. But here, if you assume infinite degrees of freedom (i.e. do a liberal z-test) you still would not be making it to two-sided significance at $\alpha = .05$.

Nevertheless, there could still be an interaction between the site's $PC1$ value and individual level $TL$. You assess this with an interaction term which you may call a cross level interaction. You can place it in your model by adding a term, or by changing your existing terms.

You could add TL:PC1 or you could change TL + PC1 to TL * PC1. Both modifications would produce the same model. The resulting additional coefficient would tell you how $TL$'s effect on $Y$ changes with $PC1$.

You have to be more specific about what you mean when you ask how to visualize the model.