> # Load required library: > library(nlme) > > ## Read in data from FLW webpage: > MIT = read.table("http://www.hsph.harvard.edu/fitzmaur/ala2e/fat.dat", header=FALSE) > names(MIT) = c("ID","Age","Age.M","Time.M","PBF") > MIT.grp = groupedData(PBF ~ Time.M | ID, data = MIT, + labels=list(x = "Time Relative to Menarche", y = "Percent Body Fat"), units = list(x = "(years)")) > # Attach the long form of the data set: > attach(MIT.grp) > Time.M0 = Time.M*(Time.M>=0) > mod1 = lme(PBF ~ Time.M + Time.M0, random = ~ Time.M + Time.M0 | ID) > summary(mod1) Linear mixed-effects model fit by REML Data: NULL AIC BIC logLik 6082.401 6131.929 -3031.201 Random effects: Formula: ~Time.M + Time.M0 | ID Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 6.777995 (Intr) Time.M Time.M 1.277142 0.292 Time.M0 1.658222 -0.544 -0.827 Residual 3.077859 Fixed effects: PBF ~ Time.M + Time.M0 Value Std.Error DF t-value p-value (Intercept) 21.361383 0.5645551 885 37.83755 0.0000 Time.M 0.417113 0.1571572 885 2.65411 0.0081 Time.M0 2.047139 0.2279685 885 8.97992 0.0000 Correlation: (Intr) Time.M Time.M 0.351 Time.M0 -0.515 -0.872 Standardized Within-Group Residuals: Min Q1 Med Q3 -2.77421523 -0.59003208 -0.03592615 0.59464973 Max 3.37978348 Number of Observations: 1049 Number of Groups: 162 > res.pop = PBF - fitted(mod1, level=0) > res.raw = residuals(mod1, type="response") > res.std = residuals(mod1, type="pearson") > res.norm = residuals(mod1, type="normalized") > ## Plot of standardized residuals vs. Fitted values: > plot(mod1) # Random scatter > plot(res.std ~ fitted(mod1), xlab="Fitted Values", ylab="Std. Residuals") > abline(h=0) > lines(lowess(res.std~fitted(mod1)),col="red",lwd=3) > plot(mod1, resid(.) ~ Time.M, abline=0) > plot(mod1, resid(.) ~ Time.M, abline=0) > plot(res.std~Time.M) > abline(h=0) > lines(lowess(res.std~Time.M),col="red",lwd=3) > plot(res.std~Time.M0) > abline(h=0) > lines(lowess(res.std~Time.M0),col="red",lwd=3) > plot(mod1, resid(.) ~ Time.M, abline=0) > plot(res.std~Time.M) > abline(h=0) > lines(lowess(res.std~Time.M),col="red",lwd=3) > plot(res.std~Time.M0) > abline(h=0) > lines(lowess(res.std~Time.M0),col="red",lwd=3) > hist(res.pop, breaks=20, xlim=c(-20,20), col="gray", freq=F, + xlab="Raw Residuals", ylab="Density") > # Add normal density: > curve(dnorm(x,mean(res.pop),sd(res.pop)),add=T) > qqnorm(res.std) > qqline(res.std) > plot(mod1, resid(.) ~ Time.M, abline=0) > plot(res.std~Time.M) > abline(h=0) > lines(lowess(res.std~Time.M),col="red",lwd=3) > identify(res.std~Time.M) [1] 4 502 531 841 > MIT.grp[c(4,502,531,841)] Error in `[.data.frame`(MIT.grp, c(4, 502, 531, 841)) : undefined columns selected > MIT.grp[c(4,502,531,841),] Grouped Data: PBF ~ Time.M | ID ID Age Age.M Time.M PBF 4 1 12.19 13.19 -1.00 23.23 502 74 12.96 13.80 -0.84 32.75 531 79 11.80 13.53 -1.73 36.94 841 128 11.47 13.31 -1.84 44.14 > MIT.grp[MIT.grp$ID==1,] Grouped Data: PBF ~ Time.M | ID ID Age Age.M Time.M PBF 1 1 9.32 13.19 -3.87 7.94 2 1 10.33 13.19 -2.86 15.65 3 1 11.24 13.19 -1.95 13.51 4 1 12.19 13.19 -1.00 23.23 5 1 13.24 13.19 0.05 10.52 6 1 14.24 13.19 1.05 20.45 > MIT.grp[MIT.grp$ID==402,] Grouped Data: PBF ~ Time.M | ID [1] ID Age Age.M Time.M PBF <0 rows> (or 0-length row.names) > MIT.grp[MIT.grp$ID==502,] Grouped Data: PBF ~ Time.M | ID [1] ID Age Age.M Time.M PBF <0 rows> (or 0-length row.names) > MIT.grp[MIT.grp$ID==74,] Grouped Data: PBF ~ Time.M | ID ID Age Age.M Time.M PBF 499 74 9.89 13.8 -3.91 22.89 500 74 10.88 13.8 -2.92 20.02 501 74 11.96 13.8 -1.84 19.52 502 74 12.96 13.8 -0.84 32.75 503 74 13.88 13.8 0.08 19.95 504 74 14.88 13.8 1.08 26.55 505 74 15.96 13.8 2.16 34.62 506 74 16.79 13.8 2.99 32.56 > plot(res.std~Time.M0) > abline(h=0) > lines(lowess(res.std~Time.M0),col="red",lwd=3) > ?Variogram.lme starting httpd help server ... done > Variogram(mod1, form = ~Time.M | ID) variog dist n.pairs 1 0.8323425 0.860 165 2 0.9310979 0.920 139 3 0.7722093 1.000 246 4 0.6758794 1.010 58 5 0.7628363 1.090 167 6 0.9793767 1.650 134 7 1.0127083 1.910 154 8 1.1234137 2.000 246 9 1.1888176 2.065 86 10 1.1069837 2.140 122 11 1.0970626 2.850 172 12 1.0564168 3.000 239 13 0.8972909 3.080 66 14 0.8479152 3.350 143 15 0.8381721 3.910 136 16 0.7794723 4.000 168 17 0.4991909 4.600 135 18 0.5047039 4.970 153 19 0.5132008 5.770 153 20 0.4511741 6.560 149 > plot(Variogram(mod1, form = ~Time.M | ID)) > Time.M02 = Time.M0^2 > mod2 = lme(PBF ~ Time.M + Time.M0 + Time.M02, random = ~ Time.M + Time.M0 + Time.M02| ID) > plot(Variogram(mod2, form = ~Time.M | ID)) > identify(Variogram(mod2, form = ~Time.M | ID)) warning: no point within 0.25 inches warning: no point within 0.25 inches warning: no point within 0.25 inches integer(0) > res.std2 = residuals(mod2, type="pearson") > > > plot(res.std2~Time.M0) > abline(h=0) > lines(lowess(res.std2~Time.M0),col="red",lwd=3) > > plot(ranef(mod1)) > qqnorm(mod1, ~ranef(.)) >