> ############## > # Section 1 > ############## > > library(faraway) > data(psid) > psid[1:40,] age educ sex income year person 1 31 12 M 6000 68 1 2 31 12 M 5300 69 1 3 31 12 M 5200 70 1 4 31 12 M 6900 71 1 5 31 12 M 7500 72 1 6 31 12 M 8000 73 1 7 31 12 M 8000 74 1 8 31 12 M 9600 75 1 9 31 12 M 9000 76 1 10 31 12 M 9000 77 1 11 31 12 M 23000 78 1 12 31 12 M 22000 79 1 13 31 12 M 8000 80 1 14 31 12 M 10000 81 1 15 31 12 M 21800 82 1 16 31 12 M 19000 83 1 24 29 16 M 7500 68 2 25 29 16 M 7300 69 2 26 29 16 M 9250 70 2 27 29 16 M 10300 71 2 28 29 16 M 10900 72 2 29 29 16 M 11900 73 2 30 29 16 M 13000 74 2 31 29 16 M 12900 75 2 32 29 16 M 19900 76 2 33 29 16 M 18900 77 2 34 29 16 M 19100 78 2 35 29 16 M 20300 79 2 36 29 16 M 31600 80 2 37 29 16 M 16600 81 2 38 29 16 M 28267 82 2 39 29 16 M 29900 83 2 40 29 16 M 33900 84 2 41 29 16 M 36500 85 2 42 29 16 M 38400 86 2 43 29 16 M 35250 87 2 44 29 16 M 37600 88 2 45 29 16 M 42700 89 2 46 29 16 M 41000 90 2 47 36 9 F 2450 68 3 > > ############ > # data reported income according to > # theire educ (education), gender, age in 1968. > ############ > > library(dplyr) > psid20 <- filter(psid, person <= 20) > library(ggplot2) > ggplot(psid20, aes(x=year, y=income))+geom_line()+facet_wrap(~ person) Warning messages: 1: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 2: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 3: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 4: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 5: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 6: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 7: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 8: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 9: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. > ggplot(psid20, aes(x=year, y=income+100, group=person)) +geom_line() + facet_wrap(~ sex) + scale_y_log10() Warning messages: 1: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 2: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 3: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 4: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 5: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 6: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 7: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 8: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 9: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. > > ########### > > lmod <- lm(log(income) ~ I(year-78), subset=(person==1), psid) > coef(lmod) (Intercept) I(year - 78) 9.3999568 0.0842667 > > ########### > # Linear Model for an individual > # We treat 1978 as year 0. > ########### > > library(lme4) > ml <- lmList(log(income) ~ I(year-78) | person, psid) > intercepts <- sapply(ml,coef)[1,] > slopes <- sapply(ml,coef)[2,] > plot(intercepts,slopes,xlab="Intercept",ylab="Slope") > > ########### > # Look at the plot. > # Intercepts and slopts are both random > ########### > > psex <- psid$sex[match(1:85,psid$person)] > boxplot(split(slopes,psex)) > > ########## > # It seems that they are normal. > ########## > > t.test(slopes[psex=="M"],slopes[psex=="F"]) Welch Two Sample t-test data: slopes[psex == "M"] and slopes[psex == "F"] t = -2.3786, df = 56.736, p-value = 0.02077 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.05916871 -0.00507729 sample estimates: mean of x mean of y 0.05691046 0.08903346 > t.test(intercepts[psex=="M"],intercepts[psex=="F"]) Welch Two Sample t-test data: intercepts[psex == "M"] and intercepts[psex == "F"] t = 8.2199, df = 79.719, p-value = 3.065e-12 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.8738792 1.4322218 sample estimates: mean of x mean of y 9.382325 8.229275 > > ########### > # We test their significance. > ########### > > library(lme4) > psid$cyear <- psid$year-78 > mmod <- lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid) > sumary(mmod, digits=3) Fixed Effects: coef.est coef.se (Intercept) 7.249 0.540 cyear 0.072 0.006 sex1 -0.575 0.061 age 0.011 0.014 educ 0.104 0.021 cyear:sex1 0.013 0.006 Random Effects: Groups Name Std.Dev. Corr person (Intercept) 0.531 cyear 0.049 0.187 Residual 0.684 --- number of obs: 1661, groups: person, 85 AIC = 3842.5, DIC = 3748.4 deviance = 3785.5 > > ############## > # We fit a longitudinal model. > ############## > > library(pbkrtest) > mmod <- lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid, REML=FALSE) > mmodr <- lmer(log(income) ~ cyear + sex +age+educ+(cyear|person),psid, REML=FALSE) > KRmodcomp(mmod,mmodr) F-test with Kenward-Roger approximation; computing time: 2.34 sec. large : log(income) ~ cyear + sex + age + educ + (cyear | person) + cyear:sex small : log(income) ~ cyear + sex + age + educ + (cyear | person) stat ndf ddf F.scaling p.value Ftest 4.5376 1.0000 247.1521 1 0.03415 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Warning message: In deviance.merMod(object, ...) : deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit > confint(mmod, method="boot") Computing bootstrap confidence intervals ... 2.5 % 97.5 % .sig01 0.409613606 0.59116506 .sig02 -0.078707508 0.45849077 .sig03 0.037625486 0.05746306 .sigma 0.660078239 0.71007733 (Intercept) 6.255072093 8.30346349 cyear 0.060669284 0.08443805 sex1 -0.679018676 -0.44637779 age -0.015571949 0.03753770 educ 0.064031957 0.14709779 cyear:sex1 0.001511962 0.02457185 > > ############## > # The above is Kenward-Roger F-test. > # The p-value is 0.034115. > # The bootstrap confidence interval is given by the output. > ############## > > diagd <- fortify(mmod) > ggplot(diagd,aes(sample=.resid))+stat_qq()+facet_grid(~sex) Warning messages: 1: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 2: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 3: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 4: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 5: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 6: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 7: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 8: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 9: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. > diagd$edulevel <- cut(psid$educ,c(0,8.5,12.5,20), labels=c("lessHS","HS","moreHS")) > ggplot(diagd, aes(x=.fitted,y=.resid)) + geom_point(alpha=0.3) + geom_hline(yintercept=0) + facet_grid(~ edulevel) + xlab("Fitted") + ylab("Residuals") Warning messages: 1: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 2: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 3: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 4: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 5: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 6: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 7: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 8: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. 9: In structure(NULL, class = "waiver") : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes. Consider 'structure(list(), *)' instead. > > ############## > # Section 2 > ############## > > library(faraway) > data(vision) > vision acuity power eye subject 1 116 6/6 left 1 2 119 6/18 left 1 3 116 6/36 left 1 4 124 6/60 left 1 5 120 6/6 right 1 6 117 6/18 right 1 7 114 6/36 right 1 8 122 6/60 right 1 9 110 6/6 left 2 10 110 6/18 left 2 11 114 6/36 left 2 12 115 6/60 left 2 13 106 6/6 right 2 14 112 6/18 right 2 15 110 6/36 right 2 16 110 6/60 right 2 17 117 6/6 left 3 18 118 6/18 left 3 19 120 6/36 left 3 20 120 6/60 left 3 21 120 6/6 right 3 22 120 6/18 right 3 23 120 6/36 right 3 24 124 6/60 right 3 25 112 6/6 left 4 26 116 6/18 left 4 27 115 6/36 left 4 28 113 6/60 left 4 29 115 6/6 right 4 30 116 6/18 right 4 31 116 6/36 right 4 32 119 6/60 right 4 33 113 6/6 left 5 34 114 6/18 left 5 35 114 6/36 left 5 36 118 6/60 left 5 37 114 6/6 right 5 38 117 6/18 right 5 39 116 6/36 right 5 40 112 6/60 right 5 41 119 6/6 left 6 42 115 6/18 left 6 43 94 6/36 left 6 44 116 6/60 left 6 45 100 6/6 right 6 46 99 6/18 right 6 47 94 6/36 right 6 48 97 6/60 right 6 49 110 6/6 left 7 50 110 6/18 left 7 51 105 6/36 left 7 52 118 6/60 left 7 53 105 6/6 right 7 54 105 6/18 right 7 55 115 6/36 right 7 56 115 6/60 right 7 > vision$npower <- rep(1:4,14) > vision acuity power eye subject npower 1 116 6/6 left 1 1 2 119 6/18 left 1 2 3 116 6/36 left 1 3 4 124 6/60 left 1 4 5 120 6/6 right 1 1 6 117 6/18 right 1 2 7 114 6/36 right 1 3 8 122 6/60 right 1 4 9 110 6/6 left 2 1 10 110 6/18 left 2 2 11 114 6/36 left 2 3 12 115 6/60 left 2 4 13 106 6/6 right 2 1 14 112 6/18 right 2 2 15 110 6/36 right 2 3 16 110 6/60 right 2 4 17 117 6/6 left 3 1 18 118 6/18 left 3 2 19 120 6/36 left 3 3 20 120 6/60 left 3 4 21 120 6/6 right 3 1 22 120 6/18 right 3 2 23 120 6/36 right 3 3 24 124 6/60 right 3 4 25 112 6/6 left 4 1 26 116 6/18 left 4 2 27 115 6/36 left 4 3 28 113 6/60 left 4 4 29 115 6/6 right 4 1 30 116 6/18 right 4 2 31 116 6/36 right 4 3 32 119 6/60 right 4 4 33 113 6/6 left 5 1 34 114 6/18 left 5 2 35 114 6/36 left 5 3 36 118 6/60 left 5 4 37 114 6/6 right 5 1 38 117 6/18 right 5 2 39 116 6/36 right 5 3 40 112 6/60 right 5 4 41 119 6/6 left 6 1 42 115 6/18 left 6 2 43 94 6/36 left 6 3 44 116 6/60 left 6 4 45 100 6/6 right 6 1 46 99 6/18 right 6 2 47 94 6/36 right 6 3 48 97 6/60 right 6 4 49 110 6/6 left 7 1 50 110 6/18 left 7 2 51 105 6/36 left 7 3 52 118 6/60 left 7 4 53 105 6/6 right 7 1 54 105 6/18 right 7 2 55 115 6/36 right 7 3 56 115 6/60 right 7 4 > > ######### > > ggplot(vision, aes(y=acuity, x=npower, linetype=eye)) + geom_line() + facet_wrap(~ subject, ncol=4) + scale_x_continuous("Power",breaks=1:4,labels=c("6/6","6/18","6/36","6/60")) There were 11 warnings (use warnings() to see them) > > ########## > > mmod <- lmer(acuity~power + (1|subject) + (1|subject:eye),vision) > sumary(mmod) Fixed Effects: coef.est coef.se (Intercept) 113.41 2.03 power1 -0.77 0.94 power2 0.02 0.94 power3 -1.77 0.94 Random Effects: Groups Name Std.Dev. subject:eye (Intercept) 3.21 subject (Intercept) 4.64 Residual 4.07 --- number of obs: 56, groups: subject:eye, 14; subject, 7 AIC = 345.5, DIC = 346.9 deviance = 339.2 > > ########## > # We fit a repeat measurement model. > ########## > > 4.64^2/(4.64^2+3.21^2+4.07^2) [1] 0.4448393 > (4.64^2+3.21^2)/(4.64^2+3.21^2+4.07^2) [1] 0.6577401 > > ############ > # The aboves are variations. > ############ > > library(pbkrtest) > mmod <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,REML=FALSE) > nmod <- lmer(acuity~1+(1|subject)+(1|subject:eye),vision,REML=FALSE) > KRmodcomp(mmod, nmod) F-test with Kenward-Roger approximation; computing time: 0.05 sec. large : acuity ~ power + (1 | subject) + (1 | subject:eye) small : acuity ~ 1 + (1 | subject) + (1 | subject:eye) stat ndf ddf F.scaling p.value Ftest 4.6065 3.0000 6.5234 1 0.04808 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Warning message: In deviance.merMod(object, ...) : deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit > > ########### > > mmodr <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,REML=FALSE, subset=-43) > nmodr <- lmer(acuity~1+(1|subject)+(1|subject:eye),vision,REML=FALSE, subset=-43) > KRmodcomp(mmodr, nmodr) F-test with Kenward-Roger approximation; computing time: 0.05 sec. large : acuity ~ power + (1 | subject) + (1 | subject:eye) small : acuity ~ 1 + (1 | subject) + (1 | subject:eye) stat ndf ddf F.scaling p.value Ftest 3.859 3.000 18.294 1 0.02675 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Warning message: In deviance.merMod(object, ...) : deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit > > ########### > # The above is Kenward-Roger F-test > ########### > > op <- options(contrasts=c("contr.helmert", "contr.poly")) > mmodr <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,subset=-43) > sumary(mmodr) Fixed Effects: coef.est coef.se (Intercept) 113.79 1.76 power1 0.39 0.54 power2 0.04 0.32 power3 0.71 0.22 Random Effects: Groups Name Std.Dev. subject:eye (Intercept) 4.97 subject (Intercept) 2.87 Residual 2.88 --- number of obs: 55, groups: subject:eye, 14; subject, 7 AIC = 319.2, DIC = 308.9 deviance = 307.1 > > ########## > # Fitted by mixed effect models with repeated measurements. > ########## > > options(op) > contr.helmert(4) [,1] [,2] [,3] 1 -1 -1 -1 2 1 -1 -1 3 0 2 -1 4 0 0 3 > plot(resid(mmodr) ~ fitted(mmodr),xlab="Fitted",ylab="Residuals") > abline(h=0) > qqnorm(ranef(mmodr)$"subject:eye"[[1]],main="") >