1 Introduction

Rambling notes on deciphering linear mixed effects models using lme4 in R.

Here are the libraries we will be using.

library(readxl)
library(tidyverse)
library(ggplot2)
library(lme4)
library(lmerTest)

2 Notation & nomenclature

lm=linear model, lmer = linear mixed effect

tilde (~) = “predicted by”. For example, lm(IQ ~ Age) examines the linear relationship of IQ as predicted by Age.

REML = restricted maximum likelihood estimation; the default in lme4

3 Random vs. Fixed Effects

A fixed effect is a variable of interest. If you are interested in modeling a specific variable’s contribution to the model, enter it as a fixed effect.

A variable that is controlled/blocked is a random effect. Its variance will still be computed, but you won’t get a parameter estimate in the summary statistics. So if you wanted to manipulate dose of Tylenol on headache severity while controlling for sex, time of day, etc – Tylenol Dose is the fixed effect.

3.0.1 Random effects

The notation (1|___) denotes a random effect. The “1” means a random intercept. For example, here’s a model with a single fixed effect (Age), controlling for the random effect of where you are from (City):
how.smart <- lmer(IQ ~ Age + (1 | City), data=XX)

Another way of specifying the effects of individual participant(s) over time. (1 + Time | Person)

3.0.2 Fixed effects

Just add the variable name, nothing fancy. Here’s age as the fixed effect.
how.smart <- lmer(IQ ~ Age + (1 | City))

3.0.3 Correlated fixed effects: modeled as correlated random slopes

Let’s say you are interested in Age and Income as fixed effects, but these variables are correlated. You need to build the correlation into the model as a random effect slope since you are not interested specifically in the correlation itself as a meaningful fixed effect but you do need to partial its variance.

how.smart <- lmer(IQ ~ Age + Income + (Age | Income) + (1| City))

3.0.4 Uncorrelated fixed effects: modeled

R assumes that the fixed effects are correlated. If they’re not, you should model as uncorrelated random effects slope using the double pipe ||

Let’s say, for example, that you have two uncorrelated fixed effects – Age, ApoE4 allele presence. Age has nothing to do with APOE4 presence. Model the null correlation as an uncorrelated random slope:
how.demented <- lmer(DementiaSeverity ~ Age + ApoE4 + (Age || ApoE4))

3.0.5 Interactions among fixed effects

m <- lmer(IQ ~ Age * Income + (1 + Age | PersonID), data=XX, REML=F)
Age * Income specifies all main effects and interactions. You could also do this as: Age + Income + Age:Income

4 Summary/output of lmer

summary(m) produces a bunch of useful detail about your model
coef(m) produces parameter estimates/coefficients
fixef(m) produces parameter estimates for the fixed effects
ranef(m) extracts random effect coefficients
confint(m) produces confidence intervals

5 Explore ChickWeight

head(ChickWeight, n = 10)
## Grouped Data: weight ~ Time | Chick
##    weight Time Chick Diet
## 1      42    0     1    1
## 2      51    2     1    1
## 3      59    4     1    1
## 4      64    6     1    1
## 5      76    8     1    1
## 6      93   10     1    1
## 7     106   12     1    1
## 8     125   14     1    1
## 9     149   16     1    1
## 10    171   18     1    1

5.1 Plot chickweights by diet

ggplot(data = ChickWeight, aes(Time, weight, group = Chick, color = as.factor(Chick))) + 
    geom_line() + facet_wrap(~Diet)

5.2 Plot it smooth

ggplot(data = ChickWeight, aes(Time, weight, group = Diet, color = as.factor(Diet))) + 
    geom_smooth(method = "glm")

5.3 Run a crude lmer

chick.notgood <- lmer(weight ~ Diet + (1 | Chick), data = ChickWeight)
summary(chick.notgood)
## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Diet + (1 | Chick)
##    Data: ChickWeight
## 
## REML criterion at convergence: 6506.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8495 -0.7499 -0.1835  0.6363  3.0852 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Chick    (Intercept)  305.8   17.49   
##  Residual             4526.6   67.28   
## Number of obs: 578, groups:  Chick, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   101.63       6.04  16.826
## Diet2          20.98      10.24   2.050
## Diet3          41.32      10.24   4.036
## Diet4          33.40      10.27   3.253
## 
## Correlation of Fixed Effects:
##       (Intr) Diet2  Diet3 
## Diet2 -0.590              
## Diet3 -0.590  0.348       
## Diet4 -0.588  0.347  0.347

5.4 Run a better lmer

This just tests differences at baseline, where t=0

chick.better <- lmer(weight ~ Diet + Time + (Time | Chick), data = ChickWeight)
summary(chick.better)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ Diet + Time + (Time | Chick)
##    Data: ChickWeight
## 
## REML criterion at convergence: 4803.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7617 -0.5758 -0.0353  0.4789  3.5025 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Chick    (Intercept) 153.87   12.40         
##           Time         14.13    3.76    -0.98
##  Residual             163.45   12.78         
## Number of obs: 578, groups:  Chick, 50
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  26.3562     2.2907 62.5949  11.506   <2e-16 ***
## Diet2         2.8386     2.3627 45.6626   1.201   0.2358    
## Diet3         2.0044     2.3627 45.6626   0.848   0.4007    
## Diet4         9.2548     2.3657 45.8859   3.912   0.0003 ***
## Time          8.4438     0.5403 48.6614  15.628   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) Diet2  Diet3  Diet4 
## Diet2 -0.351                     
## Diet3 -0.351  0.344              
## Diet4 -0.350  0.343  0.343       
## Time  -0.796 -0.004 -0.004 -0.005

5.5 An even better lmer

This just tests differences at baseline, where t=0

dat <- lmer(weight ~ Diet + Time + (Diet * Time) + (Time | Chick), data = ChickWeight)
summary(dat)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ Diet + Time + (Diet * Time) + (Time | Chick)
##    Data: ChickWeight
## 
## REML criterion at convergence: 4781.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7223 -0.5672 -0.0343  0.4579  3.5184 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Chick    (Intercept) 116.91   10.812        
##           Time         10.92    3.305   -0.97
##  Residual             163.37   12.782        
## Number of obs: 578, groups:  Chick, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  33.6613     2.9192  47.5374  11.531 2.25e-15 ***
## Diet2        -5.0277     5.0108  46.1545  -1.003 0.320915    
## Diet3       -15.4109     5.0108  46.1545  -3.076 0.003525 ** 
## Diet4        -1.7504     5.0179  46.4057  -0.349 0.728786    
## Time          6.2770     0.7613  46.9901   8.245 1.10e-10 ***
## Diet2:Time    2.3321     1.3044  45.7577   1.788 0.080411 .  
## Diet3:Time    5.1459     1.3044  45.7577   3.945 0.000272 ***
## Diet4:Time    3.2550     1.3051  45.8576   2.494 0.016298 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) Diet2  Diet3  Diet4  Time   Dt2:Tm Dt3:Tm
## Diet2      -0.583                                          
## Diet3      -0.583  0.339                                   
## Diet4      -0.582  0.339  0.339                            
## Time       -0.880  0.513  0.513  0.512                     
## Diet2:Time  0.514 -0.882 -0.299 -0.299 -0.584              
## Diet3:Time  0.514 -0.299 -0.882 -0.299 -0.584  0.341       
## Diet4:Time  0.514 -0.299 -0.299 -0.882 -0.583  0.340  0.340

Here’s what the output means – comparisons relative to baseline. Group 3 - Group 1 (1 is the reference group). The intercept differs across groups because some chicks are fatter to start. There is a slope as outcome - acceleration of weight as a function of time and diet. The rate of change for Diet 1 becomes the baseline. This is the interaction term for Diet*Time. Relationship between diet, time, and weight. Is the diet influencing acceleration of weight gain across time.

5.6 change it around so that we are modeling ‘growth’ backward from the endpoint.

this analysis seems mnuch more informative toward evaluating the superiority of Diet #3

new <- ChickWeight %>% mutate(pt.time = Time - max(Time))
dat.1 <- lmer(weight ~ Diet + pt.time + (Diet * pt.time) + (pt.time | Chick), 
    data = new)
summary(dat.1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ Diet + pt.time + (Diet * pt.time) + (pt.time | Chick)
##    Data: new
## 
## REML criterion at convergence: 4781.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7223 -0.5672 -0.0343  0.4579  3.5184 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Chick    (Intercept) 3469.91  58.906       
##           pt.time       10.92   3.305   1.00
##  Residual              163.37  12.782       
## Number of obs: 578, groups:  Chick, 50
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   165.4782    13.4890  46.8222  12.268 3.13e-16 ***
## Diet2          43.9473    23.0937  45.6848   1.903 0.063354 .  
## Diet3          92.6524    23.0937  45.6848   4.012 0.000221 ***
## Diet4          66.6040    23.1028  45.7554   2.883 0.005984 ** 
## pt.time         6.2770     0.7613  46.9907   8.245 1.10e-10 ***
## Diet2:pt.time   2.3321     1.3044  45.7583   1.788 0.080409 .  
## Diet3:pt.time   5.1459     1.3044  45.7583   3.945 0.000272 ***
## Diet4:pt.time   3.2550     1.3051  45.8582   2.494 0.016298 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Diet2  Diet3  Diet4  pt.tim Dt2:p. Dt3:p.
## Diet2       -0.584                                          
## Diet3       -0.584  0.341                                   
## Diet4       -0.584  0.341  0.341                            
## pt.time      0.995 -0.581 -0.581 -0.581                     
## Diet2:pt.tm -0.581  0.995  0.339  0.339 -0.584              
## Diet3:pt.tm -0.581  0.339  0.995  0.339 -0.584  0.341       
## Diet4:pt.tm -0.580  0.339  0.339  0.995 -0.583  0.340  0.340