Jamie’s Linear Mixed Effects Model Notes
Jamie Reilly, Ph.D.
2018-11-04
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.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 15.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.3475.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.0055.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.340Here’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