Let’s reexamine the dental data set.
library(readr)
## Warning: package 'readr' was built under R version 3.4.1
dental <- read_csv("~/BiostatCourses/PublicHealthComputing/Lectures/LongitudinalData/dental.csv")
## Parsed with column specification:
## cols(
## ID = col_integer(),
## gender = col_character(),
## age = col_integer(),
## Dist = col_double()
## )
dental$gender <- factor(dental$gender,levels=c("M","F"), labels=c("Male","Female"))
dental$ID <- factor(dental$ID)
library(lattice)
xyplot(Dist~age|gender,groups=ID,data=dental,type="l")
First we fit a model with a random intercept term and fixed effects for age, gender and age*gender.
library(lme4)
## Loading required package: Matrix
dental.lm1 <- lmer(Dist ~ age*gender + (1|ID),data=dental)
summary(dental.lm1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Dist ~ age * gender + (1 | ID)
## Data: dental
##
## REML criterion at convergence: 433.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5980 -0.4546 0.0158 0.5024 3.6862
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 3.299 1.816
## Residual 1.922 1.386
## Number of obs: 108, groups: ID, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.3406 0.9813 16.652
## age 0.7844 0.0775 10.121
## genderFemale 1.0321 1.5374 0.671
## age:genderFemale -0.3048 0.1214 -2.511
##
## Correlation of Fixed Effects:
## (Intr) age gndrFm
## age -0.869
## genderFemal -0.638 0.555
## age:gndrFml 0.555 -0.638 -0.869
ranef(dental.lm1)
## $ID
## (Intercept)
## 1 -1.11090166
## 2 0.30748171
## 3 0.96212019
## 4 1.94407791
## 5 -0.01983753
## 6 -1.32911449
## 7 0.30748171
## 8 0.63480095
## 9 -1.32911449
## 10 -3.62034917
## 11 3.25335487
## 12 2.42761770
## 13 -1.39110677
## 14 -0.62736188
## 15 1.44565998
## 16 -1.71842601
## 17 1.22744715
## 18 -1.06378753
## 19 -0.95468112
## 20 0.13638302
## 21 3.95510748
## 22 -1.17289394
## 23 -0.62736188
## 24 -0.62736188
## 25 -0.08182981
## 26 0.79102150
## 27 -1.71842601
dental.pred <- predict(dental.lm1)
xyplot(dental.pred ~ age|gender,groups=ID,data=dental,type="l")
Now with a random intercept for age as well.
dental.lm2 <- lmer(Dist ~ age*gender + (1 + age|ID),data=dental)
summary(dental.lm2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Dist ~ age * gender + (1 + age | ID)
## Data: dental
##
## REML criterion at convergence: 432.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1681 -0.3859 0.0071 0.4452 3.8495
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 5.78642 2.4055
## age 0.03252 0.1803 -0.67
## Residual 1.71620 1.3100
## Number of obs: 108, groups: ID, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.3406 1.0185 16.043
## age 0.7844 0.0860 9.121
## genderFemale 1.0321 1.5957 0.647
## age:genderFemale -0.3048 0.1347 -2.262
##
## Correlation of Fixed Effects:
## (Intr) age gndrFm
## age -0.880
## genderFemal -0.638 0.562
## age:gndrFml 0.562 -0.638 -0.880
ranef(dental.lm2)
## $ID
## (Intercept) age
## 1 -0.64132078 -0.044754794
## 2 -0.66019894 0.090293455
## 3 -0.24892341 0.113564894
## 4 1.66111242 0.028212917
## 5 0.57096613 -0.054963525
## 6 -0.82630604 -0.048057899
## 7 0.05820248 0.023482825
## 8 1.41328251 -0.071778465
## 9 -0.53894548 -0.074782151
## 10 -2.98417179 -0.062697304
## 11 2.19630290 0.101480047
## 12 1.58202018 0.081009078
## 13 -1.15234103 -0.023562688
## 14 -0.43305245 -0.018682886
## 15 2.97663085 -0.140967846
## 16 -1.64533936 -0.008474155
## 17 1.35484276 -0.010649689
## 18 -0.94670327 -0.011926969
## 19 0.36707177 -0.123853489
## 20 -0.43216529 0.053007546
## 21 3.45163819 0.050682303
## 22 0.32576678 -0.140518719
## 23 -1.15145387 0.048127744
## 24 -3.88137925 0.302008140
## 25 0.67597199 -0.070554693
## 26 -0.30825033 0.103003237
## 27 -0.78325766 -0.088646912
dental.pred2 <- predict(dental.lm2)
xyplot(dental.pred2 ~ age|gender,groups=ID,data=dental,type="l")