Growth curve data on an orthdontic measurement example
Investigators at the University of North Carolina Dental School followed the growth of 27 children (16 males, 11 females) from age 8 until age 14. Every two years they measured the distance between the pituitary and the pterygomaxillary fissure, two points that are easily identified on x-ray exposures of the side of the head.
distance – distance from the pituitary to the pterygomaxillary fissure (mm) measured on x-ray images of the skull.
age – age of the subject (yr).
Subject – subject on which the measurement was made. The levels are labelled M01 to M16 for the males and F01 to F13 for the females. The ordering is by increasing average distance within sex.
Sex – sex (Male/Female)
“Strange” behavior for M09 and M13
Random intercept for females, REML
\[ distance_{ij}|\alpha_{0i} \sim N(\beta_0 + \beta_1 age_{ij} + \alpha_{0i}, \sigma^2) \]
\[\alpha_{0i} \sim N(0,\sigma^2_{\alpha_0})\]
fit1 <- lme(distance ~ age, random = ~ 1|Subject, method = "REML", subset = Sex ==
"Female")
print(summary(fit1,cor=F))## Linear mixed-effects model fit by REML
## Data: NULL
## Subset: Sex == "Female"
## AIC BIC logLik
## 149.2183 156.169 -70.60916
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 2.06847 0.7800331
##
## Fixed effects: distance ~ age
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 0.8587419 32 20.230440 0
## age 0.479545 0.0525898 32 9.118598 0
## Correlation:
## (Intr)
## age -0.674
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2736479 -0.7090164 0.1728237 0.4122128 1.6325181
##
## Number of Observations: 44
## Number of Groups: 11Random intercept for females, ML
fit2 <- lme(distance ~ age, random = ~ 1|Subject, method = "ML", subset =
Sex == "Female")
print(summary(fit2,cor=F))## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## Subset: Sex == "Female"
## AIC BIC logLik
## 146.0304 153.1672 -69.0152
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.96987 0.7681235
##
## Fixed effects: distance ~ age
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 0.8506287 32 20.423397 0
## age 0.479545 0.0530056 32 9.047078 0
## Correlation:
## (Intr)
## age -0.685
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.3056159 -0.7192392 0.1763611 0.4257994 1.6689361
##
## Number of Observations: 44
## Number of Groups: 11Differenceis between ML and REML is in estimating random part
\[ \widehat{distance}=\widehat{\beta}^t x=\widehat{\beta}_0 + \widehat{\beta}_1 age \]
\[ \widehat{\alpha}_{0i}= E(\alpha_0|\bar{y}_i)=\frac{\sigma^2_{\alpha_0}}{\sigma^2_{\alpha_0}+\sigma^2/n_i} (\bar{y}_i-\widehat{\beta}^t x_i) \]
\[residual=distance - \widehat{distance} - \widehat{\alpha}_0 \]
par(mfrow = c(1, 2),pty="s")
plot(age[Sex=="Female"], distance[Sex=="Female"])
lines(age[Sex=="Female"],fit2$fitted[,1],col="blue")
plot(fitted(fit2),residuals(fit2))
qqnorm(residuals(fit2),main = "Normal Q-Q plot for residuals")
qqnorm(ranef(fit2)[,1],main = "Normal Q-Q plot for random intercepts")
Random slope and random intercept for females
\[ distance_{ij}|\alpha_{0i},\alpha_{1i} \sim N(\beta_0 + \beta_1 age_{ij} + \alpha_{0i}+\alpha_{1i}~ age_{ij}, \sigma^2)\]
\[(\alpha_{0i},\alpha_{1i})^t \sim N(0,V)\]
fit3 <- lme(distance ~ age, random = ~ age|Subject, method = "ML", subset = Sex ==
"Female")
print(summary(fit3,cor=F))## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## Subset: Sex == "Female"
## AIC BIC logLik
## 146.5093 157.2144 -67.25463
##
## Random effects:
## Formula: ~age | Subject
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.7238448 (Intr)
## age 0.1466746 -0.298
## Residual 0.6682746
##
## Fixed effects: distance ~ age
## Value Std.Error DF t-value p-value
## (Intercept) 17.372727 0.7422722 32 23.404794 0
## age 0.479545 0.0646183 32 7.421205 0
## Correlation:
## (Intr)
## age -0.637
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.93507352 -0.47348303 0.08334371 0.45964963 1.60921704
##
## Number of Observations: 44
## Number of Groups: 11 print(anova(fit2,fit3))## Model df AIC BIC logLik Test L.Ratio p-value
## fit2 1 4 146.0304 153.1671 -69.01520
## fit3 2 6 146.5093 157.2144 -67.25463 1 vs 2 3.521127 0.1719Random intercept, both sexes (with interaction)
\[distance_{ij}|\alpha_{0i} \sim N(\beta_0 + \beta_1 age_{ij} + \beta_2 Sex_i + \beta_3 age_{ij}:Sex_i+ \alpha_{0i}, \sigma^2)\]
\[\alpha_{0i} \sim N(0,\sigma^2_{\alpha_0})\]
fit4 <- lme(distance ~ age * Sex, random = ~ 1|Subject, method = "ML")
print(summary(fit4,cor=F))## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 440.6391 456.7318 -214.3195
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.740851 1.369159
##
## Fixed effects: distance ~ age * Sex
## Value Std.Error DF t-value p-value
## (Intercept) 16.340625 0.9814310 79 16.649795 0.0000
## age 0.784375 0.0779963 79 10.056564 0.0000
## SexFemale 1.032102 1.5376069 25 0.671239 0.5082
## age:SexFemale -0.304830 0.1221968 79 -2.494580 0.0147
## Correlation:
## (Intr) age SexFml
## age -0.874
## SexFemale -0.638 0.558
## age:SexFemale 0.558 -0.638 -0.874
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.64686407 -0.46341443 0.01556892 0.52172245 3.73335102
##
## Number of Observations: 108
## Number of Groups: 27Random incercept with different variances for males and females, both sexes (including interaction)
\[distance_{ij}|\alpha_{0i},\alpha_{1i} \sim N(\beta_0 + \beta_1 age_{ij} + \beta_2 Sex_i + \beta_3 age_{ij}:Sex_i+ \alpha_{0i}+\alpha_{1i}Sex_i, \sigma^2)\]
\[(\alpha_{0i},\alpha_{1i})^t \sim N(0,V)\]
fit5 <- lme(distance ~ age * Sex, random = ~ Sex|Subject, method = "ML")
print(anova(fit4,fit5))## Model df AIC BIC logLik Test L.Ratio p-value
## fit4 1 6 440.6391 456.7318 -214.3195
## fit5 2 8 444.4321 465.8891 -214.2160 1 vs 2 0.2069813 0.9017Random intercpet and random slope (with the same variances for both sexes)
\[distance_{ij}|\alpha_{0i},\alpha_{1i} \sim N(\beta_0 + \beta_1 age_{ij} + \beta_2 Sex_i + \beta_3 age_{ij}:Sex_i+ \alpha_{0i}+\alpha_{1i}age_{ij}, \sigma^2)\]
\[(\alpha_{0i},\alpha_{1i})^t \sim N(0,V)\]
fit6 <- lme(distance~age*Sex, random = ~ age|Subject, method = "ML")
print(anova(fit4,fit6))## Model df AIC BIC logLik Test L.Ratio p-value
## fit4 1 6 440.6391 456.7318 -214.3195
## fit6 2 8 443.8060 465.2630 -213.9030 1 vs 2 0.8331072 0.6593Random intercept without Sex
\[distance_{ij}|\alpha_{0i} \sim N(\beta_0 + \beta_1 age_{ij} + \alpha_{0i}, \sigma^2)\]
\[\alpha_{0i} \sim N(0,\sigma^2_{\alpha_0})\]
fit7<-lme(distance~age,random=~1|Subject,method="ML")Final model (random intercept with the same variance for both sexes, without main effect for Sex)
\[distance_{ij}|\alpha_{0i} \sim N(\beta_0 + \beta_1 age_{ij} + \beta_3 age_{ij}:Sex_i+ \alpha_{0i}, \sigma^2)\]
\[\alpha_{0i} \sim N(0,\sigma^2_{\alpha_0})\]
fit8 <- lme(distance ~ age + age:Sex, random = ~ 1|Subject, method = "ML")
print(anova(fit4,fit8,fit7))## Model df AIC BIC logLik Test L.Ratio p-value
## fit4 1 6 440.6391 456.7318 -214.3195
## fit8 2 5 439.1059 452.5166 -214.5530 1 vs 2 0.46688 0.4944
## fit7 3 4 451.3895 462.1181 -221.6948 2 vs 3 14.28360 0.0002 print(summary(fit8,cor=F))## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 439.1059 452.5166 -214.553
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.744358 1.372178
##
## Fixed effects: distance ~ age + age:Sex
## Value Std.Error DF t-value p-value
## (Intercept) 16.76111 0.7535284 79 22.243501 0e+00
## age 0.75516 0.0645574 79 11.697504 0e+00
## age:SexFemale -0.23312 0.0591758 79 -3.939446 2e-04
## Correlation:
## (Intr) age
## age -0.811
## age:SexFemale 0.000 -0.373
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.71241387 -0.49000265 0.03652699 0.50085138 3.73675494
##
## Number of Observations: 108
## Number of Groups: 27 par(mfrow = c(1, 2),pty="s")
plot(age, distance, type="n")
text(age,distance,c("+","*")[Sex],col=c("red","blue")[Sex])
lines(age[Sex=="Female"],fit8$fitted[,1][Sex=="Female"],col="blue")
lines(age[Sex=="Male"],fit8$fitted[,1][Sex=="Male"],col="red")
plot(fitted(fit8),resid(fit8))
text(fitted(fit8)[Subject=="M09"&age==10],resid(fit8)[Subject=="M09"& age==10],c("M09"),col="blue",pos=1,cex=.75)
text(fitted(fit8)[Subject=="M09"&age==12],resid(fit8)[Subject=="M09"& age==12],c("M09"),col="blue",pos=1,cex=.75)
text(fitted(fit8)[Subject=="M13"&age==8],resid(fit8)[Subject=="M13"& age==8],c("M13"),col="blue",pos=3,cex=.75)
qqnorm(residuals(fit8),main = "Normal Q-Q plot for residuals")
qqnorm(ranef(fit8)[,1],main = "Normal Q-Q plot for random intercepts") 