Courses Notes: Applied Multilevel Modeling

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Applied Multilevel Modeling

This course covers the theory and application of multilevel statistical models. Research data in the social sciences are often grouped in ways that impact our statistical analyses (e.g., in marital studies, spouses are more similar to one another than other study participants) and lead to interesting and substantive hypotheses (e.g., how do qualities of the relationship interact with an individual’s personality?). We will focus on why these types of data are problematic for classical statistics and the advantages of multilevel approaches. The seminar will heavily emphasize the practical application of multilevel models and rely on examples to demonstrate their need and application. We will cover the general aspects of multilevel models as well as their extension to longitudinal and multivariate data.

Note that the R code below is a modified version of the R code posted at UCLA’s ATS website for Singer and Willett’s Applied Longitudinal Data Analysis book:http://www.ats.ucla.edu/stat/R/examples/alda/

本帖隐藏的内容

[Syllabus]
[Lecture Notes]
[Lecture Slides]
[R code]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Multilevel Modeling Applied Courses Course individual another similar course impact

1. ### R code to accompany Singer & Willett's ALDA book (12/27/06)
   
2. #
   
3. ################################################################################
   
4. ###
   
5. ###                         Chapter 2
   
6. ###
   
7. ################################################################################
   
8. #
   
9. ### NOTE: datasets can be downloaded from:
   
10. #
    
11. ### http://www.ats.ucla.edu/stat/examples/alda.htm
    
12. #
    
13. ### The R code below uses the "comma separated text files" linked toward the top
    
14. ### of the page.
    
15. #
    
16. #inputting and printing the tolerance "person-level" data
    
17. tolerance <- read.table(file = file.choose(), sep=",", header=TRUE)
    
18. names(tolerance) <- tolower(names(tolerance))
    
19. summary(tolerance)
    
20. str(tolerance)
    
21. print(tolerance)
    
22. 
    
23. ### NOTE: perhaps something on converting wide to long formats using reshape lib
    
24. #
    
25. tol.df2 <- reshape(tolerance, idvar = "id", direction = "long",
    
26.                   varying = list(names(tolerance)[2:6]), v.name = "tolerance")
    
27. tol.df2 <- reshape(tolerance, idvar = "id", direction = "long",
    
28.                   varying = list(c("tol11","tol12","tol13","tol14","tol15")), v.name = "tolerance")
    
29. tol.df2$age <- tol.df2$time + 10
    
30. ### see reshape library as well (and functions melt() and cast()
    
31. #
    
32. ### re-order data.frame by id and time
    
33. tol.df2 <- with(tol.df2, tol.df2[order(id, time),])
    
34. 
    
35. ### also see sort_df in reshape library
    
36. 
    
37. #inputting and printing the tolerance "person-level" data
    
38. tolerance.pp <- read.table(file = file.choose(), sep=",", header=T)
    
39. names(tolerance.pp) <- tolower(names(tolerance.pp))
    
40. summary(tolerance.pp)
    
41. 
    
42. print(tolerance.pp[c(1:9, 76:80), ])
    
43. 
    
44. #    Table 2.1, p. 20.
    
45. #    Bivariate correlation among tolerance scores assessed on five occasions.
    
46. 
    
47. round(cor(tolerance[ , 2:6]), 2)
    
48. 
    
49. #    Fig. 2.2, p. 25.
    
50. #    Empirical growth plots.
    
51. 
    
52. library(lattice)  # needed for xyplot()
    
53. tolerance.pp$id <- factor(tolerance.pp$id)
    
54. ### NOTE: changing id to factor will print id levels in panels
    
55. 
    
56. trellis.device(theme = col.whitebg)
    
57. xyplot(tolerance ~ age | id, data = tolerance.pp, as.table = TRUE, ,
    
58.       xlab = "AGE", ylab = "TOL", pch = 16, ylim = c(0, 4))
    
59. 
    
60. #    Fig. 2.3, p. 27.
    
61. #    Smooth nonparametric trajectories superimposed on empirical growth plots.
    
62. 
    
63. ### NOTE: following won't work with panel.loess due to small n
    
64. xyplot(tolerance ~ age | id, data = tolerance.pp,
    
65.   panel = function(x, y){
    
66.     panel.xyplot(x, y, pch = 16)
    
67.     panel.lines(smooth.spline(x, y, df=4), type = "l", col = "red", lwd = 2)
    
68.   }, xlab = "AGE", ylab = "TOL", ylim=c(0, 4), as.table=T)
    
69. 
    
70. xyplot(tolerance ~ age | male, data = tolerance.pp, as.table = TRUE, ,
    
71.       xlab = "AGE", ylab = "TOL", pch = 16, ylim = c(0, 4),
    
72.       type = c("p","smooth"))
    
73. #    Table 2.2, p. 30.
    
74. #    Fitting separate within person OLS regression models.
    
75. 
    
76. ### alternative using lmList
    
77. library(nlme)
    
78. 
    
79. ### fit regressions by id
    
80. tol.lis <- lmList(tolerance ~ time | id, data = tolerance.pp)
    
81. 
    
82. ### get individuals coefficients "augmented" by other variables in data.frame
    
83. tol.lis.aug <- coef(tol.lis, augFrame = TRUE)
    
84. names(tol.lis.aug)[1] <- "intercept"
    
85. print(tol.lis.aug)
    
86. 
    
87. ### summary of individual fits
    
88. tol.lis.sum <- summary(tol.lis)
    
89. print(tol.lis.sum)
    
90. str(tol.lis.sum)
    
91. 
    
92. ### plot 95% CI around intercepts and slopes
    
93. plot(intervals(tol.lis))
    
94. 
    
95. #    Fig. 2.4, p. 31.
    
96. 
    
97. ### stem and leaf of intercepts
    
98. stem(coef(tol.lis)[[1]])
    
99. 
    
100. ### stem and leaf of slopes
     
101. stem(coef(tol.lis)[[2]])
     
102. 
     
103. # stem plot for R sq
     
104. names(tol.lis.sum)
     
105. stem(tol.lis.sum$r.squared)
     
106. 
     
107. #    Fig. 2.5, p. 32.
     
108. #    Fitted OLS trajectories superimposed on the empirical growth plots.
     
109. 
     
110. trellis.device(theme = col.whitebg)
     
111. xyplot(tolerance ~ age | id, data = tolerance.pp,
     
112.   panel = function(x, y){
     
113.     panel.xyplot(x, y)
     
114.     panel.lmline(x, y)
     
115.   }, xlab = "AGE", ylab = "TOL", ylim=c(0, 4), as.table=T)
     
116. 
     
117. #    Fig. 2.6, p. 34.
     
118. #    The collection of the trajectories of the raw data in the top panel; in the bottom panel is the collection of fitted OLS trajectories.
     
119. 
     
120. #plot of the raw data
     
121. trellis.device(theme = col.whitebg)
     
122. xyplot(tolerance ~ age, data = tolerance.pp, groups = id,
     
123.         panel = panel.superpose, panel.groups = panel.lmline,
     
124.         xlab = "AGE", ylab = "TOL", ylim=c(0, 4))
     
125. 
     
126. ### alternative with mean line
     
127. trellis.device(theme = col.whitebg)
     
128. xyplot(tolerance ~ age, data = tolerance.pp, groups = id,
     
129.         panel = function(x, y, subscripts, groups){
     
130.             panel.grid()
     
131.             panel.superpose(x, y, subscripts, groups, panel.groups = "panel.lmline")
     
132.             panel.lmline(x,y, lwd=3)
     
133.             },
     
134.         xlab = "AGE", ylab = "TOL", ylim=c(0, 4)) # auto.key = TRUE adds key but too big
     
135. 
     
136. #    Table 2.3, p. 37
     
137. #    Descriptive statistics of the estimates obtained by fitting the linear model by id.
     
138. 
     
139. #obtaining the intercepts from linear model by id
     
140. mean(coef(tol.lis)[[1]])
     
141. sd(coef(tol.lis)[[1]])
     
142. 
     
143. #obtaining the slopes from linear model by id
     
144. mean(coef(tol.lis)[[2]])
     
145. sd(coef(tol.lis)[[2]])
     
146. cor(coef(tol.lis))
     
147. 
     
148. #    Fig. 2.7, 38
     
149. #    OLS fitted trajectories separated by levels of selected predictors.
     
150. #    The two first panels are separated by gender; the two last are separated by levels of exposure.
     
151. 
     
152. ### convert to categorical variables (ie, factors) to use in trellis plots
     
153. tolerance.pp$male <- factor(tolerance.pp$male, 0:1, c("Male","Female"))
     
154. tolerance.pp$exp.cut <- tolerance.pp$exposure < 1.145
     
155. tolerance.pp$exp.cut <- factor(tolerance.pp$exp.cut, c(FALSE, TRUE), c("Low","High"))
     
156. 
     
157. trellis.device(theme = col.whitebg)
     
158. xyplot(tolerance ~ age | male, data = tolerance.pp, groups = id,
     
159.       xlab = "AGE", ylab = "TOL", type = "r", lwd = 2)
     
160. 
     
161. ### with mean line
     
162. xyplot(tolerance ~ age | male, data = tolerance.pp, groups = id,
     
163.         panel = function(x, y, subscripts, groups){
     
164.             panel.grid()
     
165.             panel.superpose(x, y, subscripts, groups,
     
166.                             panel.groups = "panel.lmline")
     
167.             panel.lmline(x,y, lwd=3)
     
168.             }, xlab = "AGE", ylab = "TOL", ylim=c(0, 4))
     
169. 
     
170. trellis.device(theme = col.whitebg)
     
171. xyplot(tolerance ~ age | exp.cut, data = tolerance.pp, groups = id,
     
172.       xlab = "AGE", ylab = "TOL", type = "r", lwd = 2)
     
173. 
     
174. ### with mean line
     
175. xyplot(tolerance ~ age | exp.cut, data = tolerance.pp, groups = id,
     
176.         panel = function(x, y, subscripts, groups){
     
177.             panel.grid()
     
178.             panel.superpose(x, y, subscripts, groups,
     
179.                             panel.groups = "panel.lmline")
     
180.             panel.lmline(x,y, lwd=3)
     
181.             }, xlab = "AGE", ylab = "TOL", ylim=c(0, 4))
     
182. 
     
183. #    Fig. 2.8, p. 40.
     
184. #    OLS estimates plotted against the predictors male and exposure.
     
185. 
     
186. par(mfrow=c(2,2))
     
187. plot(intercept ~ exposure, data = tol.lis.aug, pch = 16)
     
188. plot(intercept ~ male, data = tol.lis.aug, col = "yellow") # boxplot by default
     
189. plot(time ~ exposure, data = tol.lis.aug, pch = 16)
     
190. plot(time ~ male, data = tol.lis.aug, col = "yellow") # boxplot by default
     
191. 
     
192. ### correlations between OLS intecepts/slopes and predictors
     
193. with(tol.lis.aug, cor(as.numeric(male), intercept))
     
194. with(tol.lis.aug, cor(as.numeric(male), time))
     
195. with(tol.lis.aug, cor(exposure, intercept))
     
196. with(tol.lis.aug, cor(exposure, time))
     
197. 
     
198. ### NOTE: something on reliability, precision, and variance?

复制代码

1. ### R code to accompany Singer & Willett's ALDA book (12/27/06)
   
2. #
   
3. ################################################################################
   
4. ###
   
5. ###                         Chapter 3
   
6. ###
   
7. ################################################################################
   
8. 
   
9. #    Inputting and printing the early intervention data set, table 3.1, p. 48.
   
10. 
    
11. # reading in the opposites data
    
12. #
    
13. ### NOTE: the early intervention data used in Chapter 3 was not "released" to the
    
14. ###       public by the researchers; thus, the data below use the "opposites" data
    
15. ###       from Chapter 7, but the analyses map on directly to those done
    
16. opposites <- read.table(file.choose(), header = TRUE, sep=",")
    
17. names(opposites) <- tolower(names(opposites))
    
18. summary(opposites)
    
19. 
    
20. ### "cut" cog into low/high to mimic program variable in early intervention data
    
21. opposites$cog2 <- cut(opposites$cog, breaks = 2, labels = c("Low","High"))
    
22. 
    
23. ### change id to factor for correct labels in xyplot strips
    
24. opposites$id <- factor(opposites$id)
    
25. 
    
26. #Table 3.1, p. 48
    
27. print(opposites[1:12,])
    
28. 
    
29. #    Fig. 3.1, p. 50.
    
30. #    OLS trajectories superimposed on the empirical growth plots.
    
31. 
    
32. library(lattice)
    
33. 
    
34. trellis.device(device = pdf, height = 8.5, width = 11,
    
35.     file = "My plots.pdf", theme = col.whitebg)
    
36. xyplot(opp ~ time | id, data = opposites, type = c("g","p","r"),
    
37.         as.table = TRUE,
    
38.         layout = c(4,3,3))
    
39. ### insert more graphs...
    
40. dev.off()
    
41. 
    
42. ### for a more useful ordering, order by fitted intercept
    
43. 
    
44. trellis.device(theme = col.whitebg)
    
45. xyplot(opp ~ time | id, data = opposites, type = c("g","p","r"),
    
46.         as.table = TRUE, index.cond = function(x, y) coef(lm(y ~ x))[1])
    
47. 
    
48. ### same, but order by slopes
    
49. xyplot(opp ~ time | id, data = opposites, type = c("g","p","r"),
    
50.         as.table = TRUE, index.cond = function(x, y) coef(lm(y ~ x))[2])
    
51. 
    
52. #    Fig. 3.3, p. 57.
    
53. #    Fitted OLS trajectories and stem plots of fitted initial status and fitted rate of change by id.
    
54. 
    
55. # fitting the linear model by id
    
56. 
    
57. library(nlme)
    
58. lis.3.3 <- lmList(opp ~ time | id, data = opposites)
    
59. coef.3.3 <- coef(lis.3.3, augFrame = TRUE)
    
60. sum.3.3 <- summary(lis.3.3)
    
61. 
    
62. #plotting the linear fit by id
    
63. trellis.device(theme = col.whitebg)
    
64. xyplot(opp ~ time, data = opposites, groups = id,
    
65.         panel = function(x, y, subscripts, groups){
    
66.             panel.grid()
    
67.             panel.superpose(x, y, subscripts, groups, panel.groups = "panel.lmline")
    
68.             panel.lmline(x,y, lwd=3)
    
69.             },
    
70.         xlab = "OPP", ylab = "TIME") # auto.key = TRUE adds key but too big
    
71. 
    
72. # stem plot for fitted initial value
    
73. stem(coef.3.3[,1])
    
74. 
    
75. #stem plot for fitted rate of change
    
76. stem(coef.3.3[,2])
    
77. 
    
78. #stem plot for sigma.sq
    
79. stem(sum.3.3$r.squared)
    
80. 
    
81. #    Fig. 3.4, p. 59.
    
82. #    The top panel represents fitted OLS trajectories for program=0;
    
83. #    the bottom panel represents fitted OLS trajectories for program=1.
    
84. 
    
85. xyplot(opp ~ time | cog2, data = opposites, groups = id,
    
86.        type = c("g","r"))
    
87. xyplot(opp ~ time | cog2, data = opposites, groups = id,
    
88.         panel = function(x, y, subscripts, groups){
    
89.             panel.grid()
    
90.             panel.superpose(x, y, subscripts, groups, panel.groups = "panel.lmline")
    
91.             panel.loess(x,y, lwd=3)
    
92.             },
    
93.         xlab = "OPP", ylab = "TIME") # auto.key = TRUE adds key but too big
    
94. 
    
95. #    Table 3.3, p. 69
    
96. 
    
97. model.3.3 <- lme(opp ~ cog2*time, data = opposites,
    
98.                 random = ~ time | id, method = "ML",
    
99.                 control = list(msVerbose = TRUE, niterEM = 0, opt = "nlminb"))
    
100. summary(model.3.3)
     
101. 
     
102. #    Fig. 3.5 on page 71.
     
103. 
     
104. pred.3.3 <- expand.grid(time = 0:3, cog2 = c("Low","High"))
     
105. pred.3.3$pred <- predict(model.3.3, pred.3.3, level = 0)
     
106. 
     
107. library(lattice)
     
108. trellis.device(theme = col.whitebg)
     
109. xyplot(pred ~ time, data = pred.3.3, groups = cog2,
     
110.       type = c("g","l"), lwd = 2, auto.key = list(points = FALSE, lines = TRUE))
     
111. 
     
112. 
     
113. model.3.3.1 <- update(model.3.3, opp ~ -1 + cog2/time)
     
114. summary(model.3.3.1)
     
115. 
     
116. library(gmodels)
     
117. fixef(model.3.3.1)
     
118. c.mat <- rbind(
     
119.             'H vs. L at T0'   = c(1, -1, 0, 0),
     
120.             'H vs. L at T1'   = c(1, -1, 1, -1),
     
121.             'H vs. L at T2'   = c(1, -1, 2, -2),
     
122.             'H vs. L at T3'   = c(1, -1, 3, -3),
     
123.             'H vs. L at T4'   = c(1, -1, 4, -4)
     
124.             )
     
125. estimable(model.3.3.1, c.mat)

复制代码

1. ### R code to accompany Singer & Willett's ALDA book (12/27/06)
   
2. #
   
3. ################################################################################
   
4. ###
   
5. ###                         Chapter 4
   
6. ###
   
7. ################################################################################
   
8. 
   
9. alcohol1 <- read.table(file = file.choose(), header = TRUE, sep=",")
   
10. names(alcohol1) <- tolower(names(alcohol1))
    
11. alcohol1$id <- factor(alcohol1$id)
    
12. alcohol1$coa <- factor(alcohol1$coa, 0:1, c("Not CoA", "CoA"))
    
13. 
    
14. #    Fig. 4.1, p. 77.
    
15. #    Empirical growth plots with superimposed OLS trajectories.
    
16. 
    
17. ### how to subset a data.frame
    
18. mysub <- sample(levels(alcohol1$id), 5, replace = FALSE)
    
19. sample(1:20, 5)
    
20. ### subset the whole dataset by groups in "mysub"
    
21. mysub.df <- subset(alcohol1, id %in% mysub)
    
22. mysub.df <- subset(alcohol1, coa == "CoA")
    
23. 
    
24. library(lattice)
    
25. trellis.device(theme = col.whitebg)
    
26. xyplot(alcuse ~ age | id,
    
27.         data = alcohol1[alcohol1$id %in% c(4, 14, 23, 32, 41, 56, 65, 82), ],
    
28.         type = c("p","r"),
    
29.         ylim=c(-1, 4), as.table=T)
    
30. 
    
31. xyplot(alcuse ~ age | id,
    
32.         data = alcohol1,
    
33.         type = c("p","r"),
    
34.         ylim=c(-1, 4), as.table=T,
    
35.         subset = alcohol1$id %in% mysub)
    
36. 
    
37. 
    
38. #    Fig. 4.2, p.79.
    
39. #    Fitted OLS trajectories displayed separately by coa status and peer levels.
    
40. 
    
41. ### alternative
    
42. library(nlme)
    
43. alcuse.lis <- lmList(alcuse ~ age_14 | id, data = alcohol1)
    
44. alcuse.lis.aug <- coef(alcuse.lis, augFrame = TRUE)
    
45. names(alcuse.lis.aug)[1] <- "intercept"
    
46. 
    
47. trellis.device(theme = col.whitebg)
    
48. xyplot(alcuse ~ age | coa,
    
49.       data = alcohol1,
    
50.       groups = id,
    
51.       type = "r", lwd = 2)
    
52. 
    
53. ### alternate for peer plot - coplot
    
54. coplot(alcuse ~ jitter(age) | peer, data = alcohol1,
    
55.       panel = panel.smooth, lwd=2, col = "red", span=1)
    
56.       
    
57. #    Table 4.1, p. 94-95.
    
58. 
    
59. # Model A
    
60. model.a <- lme(alcuse ~ 1, data = alcohol1, random = ~1 |id,
    
61.               method = "ML")
    
62. summary(model.a)
    
63. VarCorr(model.a)
    
64. 
    
65. #Model B
    
66. model.b <- lme(alcuse ~ age_14 , data = alcohol1, random = ~ age_14 | id,
    
67.            method = "ML")
    
68. summary(model.b)
    
69. 
    
70. #Model C
    
71. model.c <- update(model.b, alcuse ~ coa*age_14)
    
72. summary(model.c)
    
73. 
    
74. #Model D
    
75. model.d <- update(model.b, alcuse ~ age_14*(coa + peer))
    
76. summary(model.d)
    
77. 
    
78. #Model E
    
79. model.e <- update(model.b, alcuse ~ coa + peer*age_14)
    
80. summary(model.e)
    
81. 
    
82. #Model F
    
83. model.f <- update(model.b, alcuse ~ coa + cpeer*age_14)
    
84. summary(model.f)
    
85. 
    
86. #Model G
    
87. model.g <- update(model.b, alcuse ~ ccoa + cpeer*age_14)
    
88. summary(model.g)
    
89. 
    
90. ### compare models
    
91. anova(model.a, model.b, model.c, model.d, model.e, model.f, model.g)
    
92. 
    
93. #    Fig. 4.3, p. 99.
    
94. #    Model B
    
95. 
    
96. ### alternative
    
97. #
    
98. ### could use predict for all these...
    
99. #
    
100. ### model b
     
101. pred.b <- data.frame(age_14 = 0:2)
     
102. pred.b$pred <- predict(model.b, pred.b, level=0)
     
103. plot(pred ~ I(age_14 + 14), data = pred.b, type = "l", lwd = 2,
     
104.     ylim = c(0, 2), xlim = c(13, 17),
     
105.     ylab = "Alcuse", xlab = "Age")
     
106. 
     
107. ### model c
     
108. pred.c <- expand.grid(age_14 = 0:2, coa = c("Not CoA","CoA"))
     
109. pred.c$pred <- predict(model.c, pred.c, level=0)
     
110. xyplot(pred ~ I(age_14 + 14), data = pred.c, type = "l", lwd = 2,
     
111.       groups = coa, lty = 1:2, auto.key = list(points = FALSE, lines = TRUE),
     
112.       ylim = c(0, 2), xlim = c(13, 17),
     
113.       ylab = "Alcuse", xlab = "Age")
     
114. 
     
115. ### model e
     
116. #
     
117. ### peer is continuous, so let's plot regression surface
     
118. quantile(alcohol1$peer, c(0.05, 0.95))
     
119. pred.e <- expand.grid(age_14 = 0:2, coa = c("Not CoA","CoA"),
     
120.                       peer = seq(0, 2, 0.1))
     
121. pred.e$pred <- predict(model.e, pred.e, level=0)
     
122. trellis.device(theme = col.whitebg)
     
123. wireframe(pred ~ I(age_14 + 14) + peer, data = pred.e, type = "l", lwd = 2,
     
124.       groups = coa, lty = 1:2, auto.key = list(points = FALSE, lines = TRUE),
     
125.       ylim = c(0, 2), xlim = c(13, 17),
     
126.       zlab = "Alcuse", ylab = "Peer", xlab = "Age")
     
127. 
     
128. ### NOTE: would be easier to see with dynamic 3D graph, eg, rgl or rggobi
     
129. 
     
130. 
     
131. #    Fig. 4.4
     
132. par(mfrow=c(2,2))
     
133. plot(intercept ~ coa, data = alcuse.lis.aug, col = "yellow")
     
134. plot(intercept ~ peer, data = alcuse.lis.aug, pch=16)
     
135. abline(h=0)
     
136. plot(age_14 ~ coa, data = alcuse.lis.aug, col = "yellow")
     
137. plot(age_14 ~ peer, data = alcuse.lis.aug, pch=16)
     
138. abline(h=0)
     
139. 
     
140. #    Fig. 4.5, p. 131.
     
141. #    Normality assumption plots.
     
142. 
     
143. ### alternative
     
144. qqnorm(model.f, ~ resid(., type = "r"))
     
145. 
     
146. #    Upper right panel
     
147. 
     
148. ### alternative
     
149. plot(model.f, resid(., type = "p") ~ as.numeric(id), abline = 0, pch = 16)
     
150. plot(model.f, resid(., type = "p") ~ fitted(.), abline = 0, pch = 16)
     
151. 
     
152. ### NOTE: id was converted to a factor, and plot() wants a numeric var
     
153. #    Middle left panel
     
154. 
     
155. ### alternative
     
156. qqnorm(model.f, ~ranef(.))
     
157. 
     
158. #    Middle right panel
     
159. 
     
160. ### alternative
     
161. qqnorm(model.f, ~ranef(.), standard = TRUE, abline = c(0,1))
     
162. apply(ranef(model.f, standard = TRUE), 2, mean)
     
163. apply(ranef(model.f, standard = TRUE), 2, sd)
     
164. 
     
165. ### NOTE: standardized REs don't have SD of 1...?
     
166. 
     
167. plot(model.f, ranef(., standard = TRUE)[[1]] ~ as.numeric(id), abline = 0, pch = 16)
     
168. plot(model.f, ranef(., standard = TRUE)[[2]] ~ as.numeric(id), abline = 0, pch = 16)
     
169. 
     
170. #    Fig. 4.6, p. 133.
     
171. #    Homoscedasticity plots.
     
172. #
     
173. 
     
174. ### alternative
     
175. plot(model.f, resid(.) ~ age, xlim = c(13,17), abline = 0,
     
176.     type = c("p","smooth"), lwd = 2)
     
177. ### NOTE: smoother suggests some nonlinearity
     
178. 
     
179. plot(model.f, ranef(.)[[1]] ~ coa, abline = 0, pch=16)
     
180. plot(model.f, ranef(.)[[1]] ~ peer, abline = 0, pch=16,
     
181.     type = c("p","smooth"), lwd = 2)
     
182. plot(model.f, ranef(.)[[2]] ~ coa, abline = 0, pch=16)
     
183. plot(model.f, ranef(.)[[2]] ~ peer, abline = 0, pch=16,
     
184.     type = c("p","smooth"), lwd = 2)
     
185. 
     
186. ###   Fig. 4.7, p. 136
     
187. plot(comparePred(alcuse.lis, model.b, primary = ~ alcohol1$age_14, length.out = 2))
     
188. ### NOTE: above not working quite right for some reason...
     
189. 
     
190. plot(augPred(model.b, level = 0:1, primary = ~ age_14, length.out = 2))

复制代码

1. ################################################################################
   
2. ###
   
3. ###                         Chapter 5
   
4. ###
   
5. ################################################################################
   
6. 
   
7. reading <- read.table(file = file.choose(), header=T, sep=",")
   
8. names(reading) <- tolower(names(reading))
   
9. reading$id <- factor(reading$id)
   
10. summary(reading)
    
11. 
    
12. #    Table 5.1, p. 141.
    
13. 
    
14. reading[reading$id %in% c(4, 27, 31, 33, 41, 49, 69, 77, 87), ]
    
15. 
    
16. #    Fig. 5.1, p. 143.
    
17. #    Empirical change plots with superimposed OLS trajectories.
    
18. #    The +'s and solid lines correspond to time using the child's target
    
19. #    age at data collection, whereas the dots and dashed lines correspond
    
20. #    to time using the child's observed age.
    
21. 
    
22. xyplot(piat ~ age | id,
    
23.       data = reading[reading$id %in% c(4, 27, 31, 33, 41, 49, 69, 77, 87), ],
    
24.         panel = function(x, y, subscripts){
    
25.                  panel.xyplot(x, y, pch=16)
    
26.                  panel.lmline(x,y, lty=4)
    
27.            panel.xyplot(reading$agegrp[subscripts], y, pch=3)
    
28.                  panel.lmline(reading$agegrp[subscripts],y)
    
29.     },
    
30.     ylim=c(0, 80), as.table=T, subscripts=T)
    
31. 
    
32. #    Creating the centered variables called agegrp.c and age.c, p. 144.
    
33. 
    
34. mat2 <- reading[ ,3:4] - 6.5
    
35. dimnames(mat2)[[2]] <- c("agegrp.c", "age.c")
    
36. reading <- cbind(reading, mat2)
    
37. 
    
38. ### or, a bit more transparently...
    
39. reading$age.c <- reading$age - 6.5
    
40. reading$agegrp.c <- reading$agegrp - 6.5
    
41. 
    
42. #    Table 5.2, p. 145.
    
43. #    Note: The degres of freedom used to in the calculations of the intercept
    
44. #    is different from the results in the book and this difference in partitioning
    
45. #    results also in slight differences in the standard error of the estimates.
    
46. 
    
47. #Using the agegrp variable.
    
48. lme.agegrp <- lme(piat ~ agegrp.c, reading, random = ~ agegrp.c | id, method="ML")
    
49. summary(lme.agegrp)
    
50. 
    
51. #Using the age variable.
    
52. lme.age <- lme(piat ~ age.c, reading, random= ~ age.c | id, method="ML")
    
53. summary(lme.age)
    
54. 
    
55. #    Inputting the wages data set.
    
56. 
    
57. wages <- read.table(file = file.choose(), header=T, sep=",")
    
58. names(wages) <- tolower(names(wages))
    
59. ### NOTE: for some reason, names have underscores and not periods
    
60. names(wages) <- gsub("_",".", names(wages))
    
61. 
    
62. wages$id <- factor(wages$id)
    
63. wages$black <- factor(wages$black, 0:1, c("White/Latino","Black"))
    
64. 
    
65. #    Table 5.3, p. 147.
    
66. 
    
67. wages[wages$id %in% c(206, 332, 1028), c(1, 3, 2, 6, 8, 10)]
    
68. 
    
69. #    Table 5.4, p. 149.
    
70. 
    
71. #Model A
    
72. model.5a <- lme(lnw ~ exper, wages, random= ~ exper | id, method="ML",
    
73.                 control = list(msVerbose = TRUE, niterEM = 200))
    
74. ### NOTE: boosting EM iterations speeds convergence
    
75. summary(model.5a)
    
76. 
    
77. #Model B
    
78. model.5b <- update(model.5a, lnw ~ exper*hgc.9 + exper*black)
    
79. summary(model.5b)
    
80. 
    
81. #Model C
    
82. model.5c <- update(model.5b, lnw ~ exper + exper:black + hgc.9)
    
83. summary(model.5c)
    
84. 
    
85. #    Fig. 5.2, p. 150.
    
86. #    Log wage trajectories for four prototypical dropouts from model C.
    
87. 
    
88. pred.5c <- expand.grid(exper = 0:11, black = c("White/Latino", "Black"),
    
89.                       hgc.9 = c(0, 3))
    
90. pred.5c$pred <- predict(model.5c, pred.5c, level=0)
    
91. pred.5c$hgc.9 <- factor(pred.5c$hgc.9, c(0,3), c("9th Grade","12th Grade"))
    
92. 
    
93. xyplot(pred ~ exper |hgc.9, data = pred.5c, groups = black, type = "l", lwd = 2,
    
94.       ylab = "LNW", xlab = "EXPER")
    
95. 
    
96. #    Inputting the small wages data set.
    
97. 
    
98. wages.small <- read.table(file.choose(), header=T, sep=",")
    
99. names(wages.small) <- tolower(names(wages.small))
    
100. ### NOTE: for some reason, names have underscores and not periods
     
101. names(wages.small) <- gsub("_",".", names(wages.small))
     
102. 
     
103. #    Table 5.5, p. 154.
     
104. 
     
105. #Model A
     
106. model.5.5a <- lme(lnw ~ exper + hgc.9 + exper:black, wages.small,
     
107.                   random = ~ exper | id, method = "ML",
     
108.                   control = list(msVerbose = TRUE, niterEM = 1000))
     
109. ### NOTE: jack up EM iterations to 1000
     
110. summary(model.5.5a)
     
111. 
     
112. #Model C
     
113. model.5.5c <- update(model.5.5a, random = ~ 1 | id)
     
114. anova(model.5.5c, model.5.5a)
     
115. 
     
116. ### Simpler model preferred
     
117. 
     
118. #    Inputting the unemployment data set.
     
119. 
     
120. unemployment <- read.table(file = file.choose(), header=T, sep=",")
     
121. names(unemployment) <- tolower(names(unemployment))
     
122. unemployment$id <- factor(unemployment$id)
     
123. 
     
124. #    Table 5.6, p. 161.
     
125. 
     
126. unemployment[unemployment$id %in% c(7589, 55697, 67641, 65441, 53782),]
     
127. 
     
128. #    Table 5.7, p. 163.
     
129. 
     
130. #Model A
     
131. model.5.7a <- lme(cesd ~ months, data = unemployment, random = ~ months | id,
     
132.               control = list(msVerbose = TRUE, niterEM = 100), method="ML")
     
133. summary(model.5.7a)
     
134. qqnorm(model.5.7a, ~ resid(., type = "r"))
     
135. hist(resid(model.5.7a), col = "yellow"))
     
136. ### NOTE: bit of skew in residuals
     
137. 
     
138. #Model B
     
139. model.5.7b <- update(model.5.7a, cesd ~ months + unemp)
     
140. summary(model.5.7b)
     
141. 
     
142. #Model C
     
143. model.5.7c <- update(model.5.7b, cesd ~ months*unemp)
     
144. summary(model.5.7c)
     
145. 
     
146. #Model D
     
147. model.5.7d <- update(model.5.7b, cesd ~ unemp + months:unemp,
     
148.                     random = ~ unemp + months:unemp | id,
     
149.                     control = list(msVerbose = TRUE, niterEM = 100,
     
150.                                     opt = "optim"))
     
151. summary(model.5.7d)
     
152. ### NOTE: interesting that switching the optimizer to "optim" vs. "nlminb"
     
153. ###       made the difference in convergence...
     
154. 
     
155. #    Fig. 5.3, p. 165
     
156. pred.5.3 <- data.frame(months = rep(0:15, 4),
     
157.                         unemp = c(rep(1, 16), rep(c(1,0), c(6,10)),
     
158.                                   rep(c(1,0), c(11,5)), rep(c(1,0,1), c(6, 5, 5))))
     
159. pred.5.3$quad <- factor(rep(1:4, each = nrow(pred.5.3)/4))
     
160. pred.5.3$pred <- predict(model.5.7b, pred.5.3, level=0)
     
161. 
     
162. trellis.device(theme = col.whitebg)
     
163. xyplot(pred ~ months | quad, data = pred.5.3, type = "l", lwd = 2,
     
164.       ylab = "CES-D", xlab = "Months since job loss")
     
165. ### NOTE: with significant futzing, we could make it identical...
     
166. 
     
167. #    Fig. 5.4, p. 167.
     
168. 
     
169. pred.5.7 <- expand.grid(months = 0:15, unemp = c(0,1))
     
170. pred.5.7$pred.b <- predict(model.5.7b, pred.5.7, level = 0)
     
171. pred.5.7$pred.c <- predict(model.5.7c, pred.5.7, level = 0)
     
172. pred.5.7$pred.d <- predict(model.5.7d, pred.5.7, level = 0)
     
173. 
     
174. xyplot(pred.b ~ months, data = pred.5.7, groups = unemp,
     
175.       type = "l", lwd = 2, lty = 1:2, auto.key = list(lines = TRUE),
     
176.       ylab = "CES-D", xlab = "Months since job loss")
     
177. xyplot(pred.c ~ months, data = pred.5.7, groups = unemp,
     
178.       type = "l", lwd = 2, lty = 1:2, auto.key = list(lines = TRUE),
     
179.       ylab = "CES-D", xlab = "Months since job loss")
     
180. xyplot(pred.d ~ months, data = pred.5.7, groups = unemp,
     
181.       type = "l", lwd = 2, lty = 1:2, auto.key = list(lines = TRUE),
     
182.       ylab = "CES-D", xlab = "Months since job loss")
     
183. 
     
184. #    Table 5.8, p. 175.
     
185. 
     
186. #Model A
     
187. model.5.8a <- lme(lnw ~ hgc.9 + ue.7 + exper + exper:black, data = wages,
     
188.                   random = ~ exper | id, method = "ML",
     
189.                   control = list(msVerbose = TRUE, niterEM = 100))
     
190. summary(model.5.8a)
     
191. 
     
192. #Model B
     
193. ### NOTE: don't see where
     
194. model.5.8b <- update(model.5.8a, lnw ~ hgc.9 + ue.mean + ue.person.centered +
     
195.                                   exper + exper:black)
     
196. summary(model.5.8b)
     
197. 
     
198. #Model C
     
199. model.5.8c <- update(model.5.8b, lnw ~ hgc.9 + ue1 + ue.centert1 + exper + exper:black)
     
200. summary(model.5.8c)
     
201. 
     
202. #    Inputting the medication data set.
     
203. 
     
204. medication <- read.table(file = file.choose(), header=T, sep=",")
     
205. names(medication) <- tolower(names(medication))
     
206. medication$id <- factor(medication$id)
     
207. 
     
208. #    Table 5.9, p. 182.
     
209. 
     
210. medication[c(1:6, 11, 16:21), c(3:8)]
     
211. 
     
212. #    Table 5.10, p. 184.
     
213. 
     
214. #Using time (Model A)
     
215. model.5.10a <- lme(pos ~ treat*time, data = medication,
     
216.                     random = ~ time | id, method = "ML")
     
217. summary(model.5.10a)
     
218. 
     
219. #Using time - 3.33 (Model B)
     
220. model.5.10b <- update(model.5.10a, pos ~ treat*time333)
     
221. summary(model.5.10b)
     
222. 
     
223. #Using time - 6.67 (Model C)
     
224. model.5.10c <- update(model.5.10b, pos ~ treat*time667)
     
225. summary(model.5.10c)
     
226. 
     
227. #    Fig. 5.5, p. 185.
     
228. #    Note: The vertical lines reflect the magnitude of the effect of treatment
     
229. #    when time is centered at different values.
     
230. 
     
231. days.seq <- seq(0, 7)
     
232. fixef.a <- fixef(model.5.10a)
     
233. 
     
234. trt <- fixef.a[[1]] + fixef.a[[2]] + (fixef.a[[3]]+fixef.a[[4]])*days.seq
     
235. cnt <- fixef.a[[1]] + fixef.a[[3]]*days.seq
     
236. 
     
237. ### NOTE: most of these plots look better with PDF or postscript
     
238. pdf("Figure 5.5.pdf", height = 8, width = 10.5)
     
239. plot(days.seq, trt, ylim=c(140, 190), xlim=c(0, 7), type="l", lwd = 2,
     
240.      xlab="Days", ylab=expression(widehat(POS)))
     
241. lines(days.seq, cnt, lty=4, lwd = 2)
     
242. legend(0, 190, c("treatment", "control"), lty=c(1, 4), lwd = 2)
     
243. segments(0, fixef.a[[1]] + fixef.a[[3]]*0, 0,
     
244.         fixef.a[[1]] + fixef.a[[2]] + (fixef.a[[3]]+fixef.a[[4]])*0, lty = 2)
     
245. segments(3.33, fixef.a[[1]] + fixef.a[[3]]*3.33, 3.33,
     
246.         fixef.a[[1]] + fixef.a[[2]] + (fixef.a[[3]]+fixef.a[[4]])*3.33, lty = 2)
     
247. segments(6.670, fixef.a[[1]] + fixef.a[[3]]*6.670, 6.670,
     
248.         fixef.a[[1]] + fixef.a[[2]] + (fixef.a[[3]]+fixef.a[[4]])*6.670, lty = 2)
     
249. dev.off()

复制代码

1. ################################################################################
   
2. ###
   
3. ###                         Chapter 6
   
4. ###
   
5. ################################################################################
   
6. 
   
7. #    Table 6.1, p. 192.
   
8. #    Creating the ged*exper variable in the wages data set.
   
9. 
   
10. wages$ged.exper <- wages$ged*wages$exper
    
11. print(wages[wages$id %in% c(206,2365,4384),c(1:5, 16)])
    
12. 
    
13. #    Table 6.2, p. 203.
    
14. 
    
15. model.6.2a <- lme(lnw ~ exper + hgc.9 + exper:black + ue.7, data = wages,
    
16.                   random = ~ exper | id, method = "ML",
    
17.                   control = list(msVerbose = TRUE, opt = "optim", trace = 1))
    
18. 2*model.6.2a$logLik
    
19. 
    
20. model.6.2b <- update(model.6.2a, lnw ~ exper + hgc.9 + exper:black + ue.7 + ged,
    
21.                     random = ~ exper + ged | id)
    
22. 2*model.6.2b$logLik
    
23. 
    
24. model.6.2c <- update(model.6.2b, random = ~ exper | id)
    
25. 2*model.6.2c$logLik
    
26. 
    
27. anova(model.6.2a, model.6.2c, model.6.2b)
    
28. 
    
29. model.6.2d <- update(model.6.2a, lnw ~ exper + hgc.9 + exper:black + ue.7 + postexp,
    
30.                       random = ~ exper + postexp | id,
    
31.                       control = list(msVerbose = TRUE, niterEM = 100, opt = "optim"))
    
32. 2*model.6.2d$logLik
    
33. 
    
34. anova(model.6.2a, model.6.2d)
    
35. 
    
36. model.6.2e <- update(model.6.2d, random = ~ exper| id)
    
37. 2*model.6.2e$logLik
    
38. 
    
39. model.6.2f <- update(model.6.2a, lnw ~ exper + hgc.9 + exper:black + ue.7 + postexp + ged,
    
40.                       random = ~ exper + postexp + ged | id,
    
41.                       control = list(msVerbose = TRUE, niterEM = 500, opt = "optim"))
    
42. ### NOTE: with big var-cov models, a bigger EM burn-in with optim boosts convergence
    
43. ###       notably
    
44. 2*model.6.2f$logLik
    
45. 
    
46. model.6.2g <- update(model.6.2f, random = ~ exper + ged | id)
    
47. 2*model.6.2g$logLik
    
48. 
    
49. model.6.2h <- update(model.6.2f, random = ~ exper + postexp | id)
    
50. 2*model.6.2h$logLik
    
51. 
    
52. model.6.2i <- update(model.6.2f, lnw ~ exper + hgc.9 + exper:black + ue.7 + ged + exper:ged,
    
53.                     random = ~ exper + ged + exper:ged | id)
    
54. 2*model.6.2i$logLik
    
55. 
    
56. model.6.2j <- update(model.6.2i, random = ~ ged + exper | id)
    
57. 2*model.6.2j$logLik
    
58. 
    
59. ### NOTE: examine all models; LR tests won't (all) be correct, but AIC/BIC can be
    
60. ###       compared
    
61. anova(model.6.2a, model.6.2b, model.6.2c, model.6.2d, model.6.2e, model.6.2f,
    
62.       model.6.2g, model.6.2h, model.6.2i, model.6.2j)
    
63. ### NOTE: AIC prefers model F, but BIC prefers several ahead of F with E its favorite
    
64. 
    
65. #    Table 6.3, p. 205.
    
66. #    Summary of model F.
    
67. 
    
68. summary(model.6.2f)
    
69. pairs(model.6.2f, pch = ".", cex = 2)
    
70. ### NOTE: notable correlations among random-effects
    
71. 
    
72. #    Fig. 6.4, p. 209.
    
73. #    We must first read in the alcohol data set and then run model e from table 4.1, p. 94.
    
74. 
    
75. ### we can do these using predict
    
76. pred.6.4e <- expand.grid(age_14 = seq(0, 2, 0.1), coa = c("Not CoA","CoA"),
    
77.                       peer = c(0.655, 1.381))
    
78. pred.6.4e$pred <- predict(model.e, pred.6.4e, level=0)
    
79. pred.6.4e$pred <- pred.6.4e$pred^2  # square it to go back to original metric
    
80. 
    
81. trellis.device(theme = col.whitebg)
    
82. xyplot(pred ~ I(age_14 + 14) | peer, data = pred.6.4e, type = "l", lwd = 2,
    
83.       groups = coa, lty = 1:2, auto.key = list(points = FALSE, lines = TRUE, lwd=2),
    
84.       xlim = c(13, 17),
    
85.       ylab = "Alcuse", xlab = "Age")
    
86. 
    
87. #reading in the alcohol data
    
88. alcohol <- read.table("d:/alcohol1_pp.txt", header=T, sep=",")
    
89. 
    
90. #table 4.1, model e
    
91. model.e <- lme(alcuse~coa+peer*age_14 , data=alcohol, random= ~ age_14 | id,
    
92.            method="ML")
    
93. #obtaining the fixed effects parameters
    
94. fixef.e <- fixef(model.e)
    
95. #obtaining the predicted values and squaring them
    
96. fit2.ec0p0 <- (fixef.e[[1]] + .655*fixef.e[[3]] +
    
97.              alcohol$age_14[1:3]*fixef.e[[4]] +
    
98.              .655*alcohol$age_14[1:3]*fixef.e[[5]])^2
    
99. fit2.ec0p1 <- (fixef.e[[1]] + 1.381*fixef.e[[3]] +
    
100.              alcohol$age_14[1:3]*fixef.e[[4]] +
     
101.              1.381*alcohol$age_14[1:3]*fixef.e[[5]] )^2
     
102. fit2.ec1p0 <- (fixef.e[[1]] + fixef.e[[2]] + .655*fixef.e[[3]] +
     
103.              alcohol$age_14[1:3]*fixef.e[[4]] +
     
104.              .655*alcohol$age_14[1:3]*fixef.e[[5]] )^2
     
105. fit2.ec1p1 <- (fixef.e[[1]] + fixef.e[[2]] + 1.381*fixef.e[[3]] +
     
106.              alcohol$age_14[1:3]*fixef.e[[4]] +
     
107.              1.381*alcohol$age_14[1:3]*fixef.e[[5]])^2
     
108. 
     
109. plot(alcohol1$age[1:3], fit2.ec0p0, ylim=c(0, 3), type="n",
     
110.      ylab="predicted alcuse squared", xlab="age")
     
111. lines(spline(alcohol1$age[1:3], fit2.ec0p0), pch=2, type="b")
     
112. lines(spline(alcohol1$age[1:3], fit2.ec0p1), type="b", pch=0)
     
113. lines(spline(alcohol1$age[1:3], fit2.ec1p0), type="b", pch=17)
     
114. lines(spline(alcohol1$age[1:3], fit2.ec1p1), type="b", pch=15)
     
115. 
     
116. title("Non-linear Change")
     
117. legend(14, 3, marks=c(2, 0, 17, 15), c("COA=0, low peer", "COA=0, high peer",
     
118.        "COA=1, low peer", "COA=1, high peer"))
     
119. 
     
120. #    Reading in the berkeley data set and then creating the transformed variables for iq and age.
     
121. 
     
122. berkeley <- read.table(file.choose(), header=T, sep=",")
     
123. names(berkeley) <- tolower(names(berkeley))
     
124. berkeley$age2.3 <- berkeley$age^(1/2.3)
     
125. berkeley$iq2.3 <- berkeley$iq^2.3
     
126. 
     
127. #    Fig. 6.6, p. 212.
     
128. 
     
129. par(mfrow=c(1,3))
     
130. plot(berkeley$age, berkeley$iq, type="p", ylim=c(0, 250), xlim=c(0, 60),
     
131.      ylab="IQ", xlab="TIME", pch=16)
     
132. plot(berkeley$age, berkeley$iq2.3, type="p", ylim=c(0, 300000), xlim=c(0, 60),
     
133.      ylab="IQ^(2.3)", xlab="TIME", pch=16)
     
134. plot(berkeley$age2.3, berkeley$iq, type="p", ylim=c(0, 250), xlim=c(0, 6),
     
135.      ylab="IQ", xlab="TIME^(1/2.3)", pch=16)
     
136. 
     
137. ### Table 6.4?
     
138. 
     
139. #    Reading in the external data set.
     
140. 
     
141. external <- read.table(file.choose(), header = TRUE, sep=",")
     
142. names(external) <- tolower(names(external))
     
143. 
     
144. #    Creating the higher order variables for grade.
     
145. #
     
146. ### NOTE: not needed...
     
147. external$grade2 <- external$grade^2
     
148. external$grade3 <- external$grade^3
     
149. external$grade4 <- external$grade^4
     
150. 
     
151. #    Fig. 6.7, p. 218.
     
152. 
     
153. par(mfrow=c(2,4))
     
154. model.6.7a <- lm(external ~ poly(grade, 2), data = external,
     
155.                   subset = id == "1")
     
156. plot(1:6, fitted(model.6.7a), type = "l", lwd = 2, col = "red",
     
157.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
158. points(1:6, subset(external, id == "1", external)[[1]], pch = 16)
     
159. lines(spline(1:6, subset(external, id == "1", external)[[1]]), lty = 3)
     
160. 
     
161. model.6.7b <- lm(external ~ poly(grade, 2), data = external,
     
162.                   subset = id == "6")
     
163. plot(1:6, fitted(model.6.7b), type = "l", lwd = 2, col = "red",
     
164.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
165. points(1:6, subset(external, id == "6", external)[[1]], pch = 16)
     
166. lines(spline(1:6, subset(external, id == "6", external)[[1]]), lty = 3)
     
167. 
     
168. model.6.7c <- lm(external ~ grade, data = external,
     
169.                   subset = id == "11")
     
170. plot(1:6, fitted(model.6.7c), type = "l", lwd = 2, col = "red",
     
171.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
172. points(1:6, subset(external, id == "11", external)[[1]], pch = 16)
     
173. lines(spline(1:6, subset(external, id == "11", external)[[1]]), lty = 3)
     
174. ### NOTE: looks like a quadratic would fit better...
     
175. 
     
176. model.6.7d <- lm(external ~ grade, data = external,
     
177.                   subset = id == "25")
     
178. plot(1:6, fitted(model.6.7d), type = "l", lwd = 2, col = "red",
     
179.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
180. points(1:6, subset(external, id == "25", external)[[1]], pch = 16)
     
181. lines(spline(1:6, subset(external, id == "25", external)[[1]]), lty = 3)
     
182. 
     
183. model.6.7e <- lm(external ~ poly(grade, 3), data = external,
     
184.                   subset = id == "34")
     
185. plot(1:6, fitted(model.6.7e), type = "l", lwd = 2, col = "red",
     
186.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
187. points(1:6, subset(external, id == "34", external)[[1]], pch = 16)
     
188. lines(spline(1:6, subset(external, id == "34", external)[[1]]), lty = 3)
     
189. 
     
190. model.6.7f <- lm(external ~ poly(grade, 4), data = external,
     
191.                   subset = id == "36")
     
192. plot(1:6, fitted(model.6.7f), type = "l", lwd = 2, col = "red",
     
193.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
194. points(1:6, subset(external, id == "36", external)[[1]], pch = 16)
     
195. lines(spline(1:6, subset(external, id == "36", external)[[1]]), lty = 3)
     
196. 
     
197. model.6.7g <- lm(external ~ poly(grade, 2), data = external,
     
198.                   subset = id == "40")
     
199. plot(1:6, fitted(model.6.7g), type = "l", lwd = 2, col = "red",
     
200.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
201. points(1:6, subset(external, id == "40", external)[[1]], pch = 16)
     
202. lines(spline(1:6, subset(external, id == "40", external)[[1]]), lty = 3)
     
203. 
     
204. model.6.7h <- lm(external ~ poly(grade, 2), data = external,
     
205.                   subset = id == "26")
     
206. plot(1:6, fitted(model.6.7h), type = "l", lwd = 2, col = "red",
     
207.     ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
208. points(1:6, subset(external, id == "26", external)[[1]], pch = 16)
     
209. lines(spline(1:6, subset(external, id == "26", external)[[1]]), lty = 3)
     
210. 
     
211. ### alternative in a for loop
     
212. ex.ids <- c("1","6","11","25","34","36","40","26")
     
213. 
     
214. par(mfrow=c(2,4))
     
215. for (i in 1:length(ex.ids)){
     
216.     fit <- lm(external ~ poly(grade, 4), data = external,
     
217.                   subset = id == ex.ids[i])
     
218.     pred <- predict(fit, data.frame(grade = seq(1, 6, 0.1)))
     
219.     plot(x = seq(1, 6, 0.1), y = pred, type = "l", lwd = 2, col = "red",
     
220.           ylim = c(0,60), ylab = "EXTERNAL", xlab = "GRADE")
     
221.     points(1:6, subset(external, id == ex.ids[i], external)[[1]], pch = 16)
     
222.     lines(spline(1:6, subset(external, id == ex.ids[i], external)[[1]]), lty = 3)
     
223.     }
     
224. ### NOTE: could use drop1() to find best fit if we wanted...
     
225. 
     
226. #    Table 6.5, p. 221.
     
227. #    Comparison of fitting alternative polynomial change trajectories to the external data set.
     
228. 
     
229. model.6.5a <- lme(external ~ 1, random =  ~ 1 | id, method = "ML", external)
     
230. summary(model.6.5a)
     
231. 
     
232. model.6.5b <- update(model.6.5a, external ~ time, random =  ~ time | id)
     
233. summary(model.6.5b)
     
234. ### NOTE: could use poly(time, degree) for orthogonal polynomials
     
235. 
     
236. model.6.5c <- update(model.6.5a, external ~ time + I(time^2),
     
237.                       random =  ~ time + I(time^2) | id,
     
238.                       control = list(msVerbose = TRUE, niterEM = 100, opt = "optim"))
     
239. summary(model.6.5c)
     
240. 
     
241. model.6.5d <- update(model.6.5c, external ~ time + I(time^2) + I(time^3),
     
242.                     random = ~ time + I(time^2) + I(time^3) | id)
     
243. summary(model.6.5d)
     
244. pairs(model.6.5d, pch = 16) # super high correlations
     
245. 
     
246. model.6.5d.1 <- update(model.6.5c, external ~ poly(time, 3),
     
247.                     random = ~ poly(time, 3) | id)
     
248. summary(model.6.5d.1)
     
249. ### NOTE: converges faster with orthogonal polynomials, which remove non-essential
     
250. ###       correlation
     
251. 
     
252. #    Reading in the fox and geese data.
     
253. 
     
254. fg.df <- read.table(file.choose(), header = TRUE, sep=",")
     
255. names(fg.df) <- tolower(names(fg.df))
     
256. fg.df$id <- factor(fg.df$id)
     
257. str(fg.df)
     
258. 
     
259. #    Fig. 6.8, p. 227.
     
260. #    Empirical growth plots for 8 children in the fox and geese data.
     
261. 
     
262. xyplot(nmoves ~ game | id,
     
263.         data=fg.df[fg.df$id %in% c(1, 4, 6, 7, 8, 11, 12, 15), ],
     
264.         ylim=c(0, 25), as.table = T, pch = 16)
     
265. 
     
266. #    Table 6.6, p. 231.
     
267. #    Fitting a logistic model to the fox and geese data.
     
268. 
     
269. model.6.6a <- nlme(nmoves ~ 1 + 19/(1 + xmid*exp( -scal*game + u)),
     
270.                     fixed = scal + xmid ~ 1,
     
271.                     random = scal + u ~ 1 |id,
     
272.                     start = c(scal = 0.2, xmid = 12),
     
273.                     data = fg.df,
     
274.                     control = list(msVerbose = TRUE, niterEM = 100))
     
275. summary(model.6.6a)
     
276. 
     
277. 
     
278. model.6.6b <- nlme(nmoves ~ 1 + 19/(1 + xmid*exp(-scal10*game - scal01*read.c - scal11*readc.game + u)),
     
279.                     fixed = scal10 + scal01 + scal11 + xmid ~ 1,
     
280.                     random = scal10 + u ~ 1 |id,
     
281.                     start = c(scal10 = 0.12, scal01 = -0.4, scal11 = 0.04, xmid = 12),
     
282.                     data = fg.df)
     
283. summary(model.b)
     
284. ### NOTE: need to check book for read.c and readc.game
     
285. 
     
286. #    Fig. 6.10, p. 232.
     
287. 
     
288. fixef.a <- fixef(model.a)
     
289. fit.a <- 1 + 19/(1 + fixef.a[[2]]*exp(-fixef.a[[1]]*fg.df$game[1:27]))
     
290. 
     
291. plot(fg.df$game[1:27], fit.a, ylim=c(0, 25), type="l",
     
292.      ylab="predicted nmoves", xlab="game")
     
293. title("Model A \n Unconditional logistic growth")
     
294. 
     
295. fixef.b <- fixef(model.b)
     
296. fit.b.high <- 1 + 19/(1+fixef.b[[4]]*exp(-fixef.b[[1]]*fg.df$game[1:27] -
     
297.                  fixef.b[[2]]*1.58 - fixef.b[[3]]*1.58*fg.df$game[1:27]))
     
298. 
     
299. fit.b.low <- 1 + 19/(1+fixef.b[[4]]*exp(-fixef.b[[1]]*fg.df$game[1:27] -
     
300.                  fixef.b[[2]]*(-1.58) - fixef.b[[3]]*(-1.58)*fg.df$game[1:27]))
     
301. 
     
302. plot(fg.df$game[1:27], fit.b.high, ylim=c(0, 25), type="l",
     
303.      ylab="predicted nmoves", xlab="game")
     
304. lines(fg.df$game[1:27], fit.b.low, lty=3)
     
305. 
     
306. title("Model B \n Fitted logistic growth by reading level")
     
307. legend(1, 25, c("High Reading level","Low reading level"), lty=c(1, 3))

复制代码

1. ################################################################################
   
2. ###
   
3. ###                         Chapter 7
   
4. ###
   
5. ################################################################################
   
6. 
   
7. opposites <- read.table(file.choose(), header = TRUE, sep=",")
   
8. names(opposites) <- tolower(names(opposites))
   
9. opposites$id <- factor(opposites$id)
   
10. 
    
11. #    Table 7.2, p. 246.
    
12. 
    
13. opp.reml <- lme(opp ~ time*ccog,
    
14.                 data = opposites,
    
15.                 random = ~ time | id,
    
16.                 control = list(msVerbose = TRUE, niterEM = 100, opt = "optim"))
    
17. summary(opp.reml)
    
18. 
    
19. #    Table 7.3, p. 258-259.
    
20. 
    
21. ### NOTE: why did S&W drop the random-effects in fitting these models?
    
22. #
    
23. #compound symmetry
    
24. fit.cs <- gls(opp ~ time * ccog, data = opposites,
    
25.               correlation = corCompSymm(form =  ~ time | id))
    
26. summary(fit.cs) ; 2 * logLik(fit.cs) ; fit.cs$sigma^2
    
27. 
    
28. #autoregressive
    
29. fit.ar1 <- gls(opp ~ time * ccog, data = opposites,
    
30.                 corr = corAR1(form =  ~ wave | id))
    
31. 2 * logLik(fit.ar1) ; fit.ar1$sigma^2
    
32. 
    
33. #    Table 7.4, p. 265.
    
34. 
    
35. # Standard error covariance structure
    
36. summary(opp.reml)
    
37. 
    
38. # Unstructured error covariance structure
    
39. opp.unstr <- update(opp.reml, correlation = corSymm(form =  ~ wave | id))
    
40. summary(opp.unstr)
    
41. 
    
42. ### Additional analyses
    
43. #
    
44. ### ACF plot from random intercept/slope model
    
45. plot(ACF(opp.reml), alpha = 0.01)
    
46. 
    
47. ### add AR1
    
48. opp.reml.ar1 <- update(opp.reml, corr = corAR1(form =  ~ wave | id))
    
49. anova(opp.reml, opp.reml.ar1)
    
50. 
    
51. ### not preferred
    
52. 
    
53. anova(fit.cs, fit.ar1, opp.reml, opp.reml.ar1, opp.unstr)
    
54. 
    
55. ### either the gls with AR1 or random intercepts/slopes
    
56. #
    
57. ### NOTE: possibly examine variogram and other correlation methods

复制代码

谢谢楼主分享！

好书，学习啦