Week 3 Review Question Constrained models
R version 3.4.4 (2018-03-15) -- "Someone to Lean On"
# Class example, for reference
> summary(sleeplmer)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.9536 -0.4634 0.0231 0.4634 5.1793
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 612.09 24.740
Days 35.07 5.922 0.07
Residual 654.94 25.592
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.405 6.825 36.838
Days 10.467 1.546 6.771
Correlation of Fixed Effects:
(Intr)
Days -0.138
## For sleepstudy the constrained model (refer to the sleeplmer in class example) that constrains
## no individual difference in time gradient (increase in reaction time per day)
> sleeplmerC = lmer(Reaction ~ 1 + Days + (1 |Subject), sleepstudy)
## In terms of our model formulation guide sleeplmerC appears in the random effects to set each individual with a separate flat line
# and then adds in a constant gradient in the fixed effects
#comparison with class examples shows same fixed effects,
#larger residual variance and variability of intercept (initial level)
# and as indicated no variability for gradient (slope) over subjects
> summary(sleeplmerC)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ 1 + Days + (1 | Subject)
Data: sleepstudy
REML criterion at convergence: 1786.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.2257 -0.5529 0.0109 0.5188 4.2506
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.2 37.12
Residual 960.5 30.99
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.4051 9.7467 25.79
Days 10.4673 0.8042 13.02
Correlation of Fixed Effects:
(Intr)
Days -0.371
# compare results by looking at confidence intervals
> confint(sleeplmer, oldNames = FALSE) # to get good labels
Computing profile confidence intervals ...
2.5 % 97.5 %
sd_(Intercept)|Subject 14.3815822 37.715996
cor_Days.(Intercept)|Subject -0.4815007 0.684986
sd_Days|Subject 3.8011641 8.753383
sigma 22.8982669 28.857997
(Intercept) 237.6806955 265.129515
Days 7.3586533 13.575919
> confint(sleeplmerC, oldNames = FALSE) # to get good labels
Computing profile confidence intervals ...
2.5 % 97.5 %
sd_(Intercept)|Subject 26.007120 52.93598
sigma 27.813847 34.59105
(Intercept) 231.992326 270.81788
Days 8.886551 12.04802
Warning message:
In optwrap(optimizer, par = start, fn = function(x) dd(mkpar(npar1, :
convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q
# which model better fits these data? clearly the class example, as our
# plots showed considerable indiviual difference in change
> anova(sleeplmerC, sleeplmer)
refitting model(s) with ML (instead of REML)
Data: sleepstudy
Models:
sleeplmerC: Reaction ~ 1 + Days + (1 | Subject)
sleeplmer: Reaction ~ Days + (1 + Days | Subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
sleeplmerC 4 1802.1 1814.8 -897.04 1794.1
sleeplmer 6 1763.9 1783.1 -875.97 1751.9 42.139 2 7.072e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Further constraint would be to require all individuals to have same growth curve
# which can be done via lm or better geee (which tries to account for correlation)
# average growth curve point estimates identical to fixed effects
> sleeplm = lm(Reaction ~ 1 + Days, sleepstudy)
> summary(sleeplm)
Call:
lm(formula = Reaction ~ 1 + Days, data = sleepstudy)
Residuals:
Min 1Q Median 3Q Max
-110.848 -27.483 1.546 26.142 139.953
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 251.405 6.610 38.033 < 2e-16 ***
Days 10.467 1.238 8.454 9.89e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 47.71 on 178 degrees of freedom
Multiple R-squared: 0.2865, Adjusted R-squared: 0.2825
F-statistic: 71.46 on 1 and 178 DF, p-value: 9.894e-15