## ----echo=FALSE---------------------------------------------------------- library("nlme") library("lattice") data(Orthodont, package = "nlme") ## ------------------------------------------------------------------------ str(Orthodont) head(Orthodont) ## ------------------------------------------------------------------------ Ort <- as.data.frame(Orthodont) str(Ort) head(Ort) ## ------------------------------------------------------------------------ Ort2 <- groupedData(formula = distance ~ age | Subject, data = Ort, outer = ~Sex, labels = list(x = "Age", y = "Distance from pituitary to pterygomaxillary fissure"), units = list(x = "[yr]", y = "[mm]")) str(Ort2) head(Ort2) ## ------------------------------------------------------------------------ summary(Ort) ## ------------------------------------------------------------------------ gsummary(Ort2) gsummary(Ort2, FUN = max) ## ------------------------------------------------------------------------ summary(gsummary(Ort2)) ## ------------------------------------------------------------------------ summary(Ort2) ## ----fig.width=15-------------------------------------------------------- plot(Ort2) ## ----fig.width=10-------------------------------------------------------- plot(Ort2, outer = TRUE, aspect = "fill") ## ----fig.width=10-------------------------------------------------------- plot(Ort2, outer = ~Sex, aspect = "fill") ## ------------------------------------------------------------------------ fitlm <- lm(distance ~ I(age - 8) * Sex, data = Ort2) summary(fitlm) ## ----fig.width = 10------------------------------------------------------ par(mfrow = c(2, 2)) plot(fitlm, pch = 23, col = "darkblue", bg = "skyblue") par(mfrow = c(1, 1)) ## ----fig.width = 10------------------------------------------------------ COL <- c("darkblue", "red") BG <- c("skyblue", "pink") PCH <- c(24, 25) names(COL) <- names(BG) <- names(PCH) <- c("Male", "Female") plot(distance ~ age, col = COL[Sex], pch = PCH[Sex], bg = BG[Sex], data = Ort2) abline(lm(distance ~ age, data = subset(Ort2, Sex == "Female")), col="red", lwd = 2) abline(lm(distance ~ age, data = subset(Ort2, Sex == "Male")), col = "blue", lwd = 2) ## ----fig.width = 8------------------------------------------------------- bwplot(Ort2$Subject ~ residuals(fitlm)) ## ------------------------------------------------------------------------ fitslm <- lmList(Ort2) ## ------------------------------------------------------------------------ fitslm <- lmList(distance ~ age | Subject, data = Ort) ## ------------------------------------------------------------------------ summary(fitslm) ## ----fig.width = 8------------------------------------------------------- plot(fitslm) ## ------------------------------------------------------------------------ fitslm2 <- lmList(distance ~ I(age - 8), data = Ort2) ## ------------------------------------------------------------------------ fitslm2 <- lmList(distance ~ I(age - 8) | Subject, data = Ort) ## ------------------------------------------------------------------------ summary(fitslm2) ## ------------------------------------------------------------------------ pairs(fitslm2, id = 0.01, adj = -0.5) ## ----fig.width = 8------------------------------------------------------- pairs(fitslm, id = 0.01, adj = -0.5) ## ----fig.width = 8------------------------------------------------------- plot(intervals(fitslm2)) ## ------------------------------------------------------------------------ Ort <- transform(Ort, resdistance = residuals(fitlm)) ## ------------------------------------------------------------------------ resfitslm <- lmList(resdistance ~ age | Subject, data = Ort) resfitslm2 <- lmList(resdistance ~ I(age - 8) | Subject, data = Ort) ## ----fig.width = 8------------------------------------------------------- plot(resfitslm) ## ----fig.width = 8------------------------------------------------------- pairs(resfitslm, id = 0.01, adj = -0.5) ## ----fig.width = 8------------------------------------------------------- pairs(resfitslm2, id = 0.01, adj = -0.5) ## ----fig.width = 8------------------------------------------------------- plot(intervals(resfitslm2)) ## ------------------------------------------------------------------------ fit.1 <- lme(fixed = distance ~ I(age - 8), random = ~ I(age - 8) | Subject, data = Ort2) ## ------------------------------------------------------------------------ summary(fit.1) ## ------------------------------------------------------------------------ fixef(fit.1) ## ------------------------------------------------------------------------ ranef(fit.1) ## ------------------------------------------------------------------------ coef(fit.1) ## ------------------------------------------------------------------------ intervals(fit.1) ## ------------------------------------------------------------------------ fit.2 <- lme(distance ~ I(age - 8)*Sex, random = ~ I(age - 8) | Subject, data = Ort2) summary(fit.2) ## ------------------------------------------------------------------------ fitted(fit.1, level = 0) ## ------------------------------------------------------------------------ fitted(fit.1, level = 1) ## ------------------------------------------------------------------------ resid(fit.1, level = 1, type = "pearson") ## ------------------------------------------------------------------------ plot(fit.1) ## ------------------------------------------------------------------------ plot(fit.1, resid(., type = "p") ~ fitted(.) | Sex, abline = 0) ## ------------------------------------------------------------------------ plot(fit.1, Subject ~ resid(.)) ## ------------------------------------------------------------------------ plot(fit.1, distance ~ fitted(.) | Subject, abline = c(0, 1)) ## ------------------------------------------------------------------------ pairs(fit.1) ## ------------------------------------------------------------------------ pairs(fit.1, ~ coef(., augFrame = TRUE) | Sex, id = 0.1, adj = -0.5) ## ------------------------------------------------------------------------ pairs(fit.1, ~ranef(.)) ## ------------------------------------------------------------------------ pairs(fit.1, ~ranef(.) | Sex) ## ------------------------------------------------------------------------ plot(ranef(fit.1)) ## ------------------------------------------------------------------------ fit1RE <- ranef(fit.1, aug = TRUE) print(fit1RE) plot(fit1RE, form = ~ Sex) ## ----fig.width = 10------------------------------------------------------ par(mfrow = c(1, 2)) hist(fit1RE[, "(Intercept)"], xlab = expression(hat(b)[0]), prob = TRUE, main = "Intercept", col = "salmon") hist(fit1RE[, "I(age - 8)"], xlab = expression(hat(b)[1]), prob = TRUE, main = "Slope", col = "salmon") par(mfrow = c(1, 1)) ## ----fig.width = 10------------------------------------------------------ par(mfrow = c(1, 2)) qqnorm(fit1RE[, "(Intercept)"], main = "Intercept", col = "red4", bg = "salmon", pch = 21) qqline(fit1RE[, "(Intercept)"], col = "darkblue", lwd = 2) qqnorm(fit1RE[, "I(age - 8)"], main = "Slope", col = "red4", bg = "salmon", pch = 21) qqline(fit1RE[, "I(age - 8)"], col = "darkblue", lwd = 2) par(mfrow = c(1, 1)) ## ------------------------------------------------------------------------ Ort2$age.m.8 <- Ort2$age - 8 fit.1x=lme(fixed= distance ~ age.m.8, random = ~ age.m.8 | Subject, data=Ort2) fit1xRE <- ranef(fit.1x, aug = TRUE) print(fit1xRE) plot(fit1xRE, form = age.m.8 ~ Sex) ## ------------------------------------------------------------------------ cmp.1 <- compareFits(coef(fitslm2), coef(fit.1)) plot(cmp.1, mark = fixef(fit.1)) ## ----fig.width=10-------------------------------------------------------- cmp.2 <- comparePred(fitslm2, fit.1, length.out = 2) plot(cmp.2, layout = c(8, 4), between = list(y = c(0, 0.5))) ## ------------------------------------------------------------------------ fit.1.uncor <- lme(fixed = distance ~ age.m.8, random = pdDiag(~ age.m.8), data = Ort2) summary(fit.1.uncor)