[Jared E. Knowles]Getting Started with Multilevel Modeling in R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

1. Getting Started with Multilevel Modeling in R
   
2. Jared E. Knowles
   
3. Introduction
   
4. Analysts dealing with grouped data and complex hierarchical structures in their data ranging from measurements nested within participants, to counties nested within states or students nested within classrooms often find themselves in need of modeling tools to reflect this structure of their data. In R there are two predominant ways to fit multilevel models that account for such structure in the data. These tutorials will show the user how to use both the lme4 package in R to fit linear and nonlinear mixed effect models, and to use rstan to fit fully Bayesian multilevel models. The focus here will be on how to fit the models in R and not the theory behind the models. For background on multilevel modeling, see the references. [1]
   
5. 
   
6. This tutorial will cover getting set up and running a few basic models using lme4 in R.Future tutorials will cover:
   
7. 
   
8. constructing varying intercept, varying slope, and varying slope and intercept models in R
   
9. generating predictions and interpreting parameters from mixed-effect models
   
10. generalized and non-linear multilevel models
    
11. fully Bayesian multilevel models fit with rstan or other MCMC methods
    
12. Setting up your enviRonment
    
13. Getting started with multilevel modeling in R is simple. lme4 is the canonical package for implementing multilevel models in R, though there are a number of packages that depend on and enhance its feature set, including Bayesian extensions. lme4 has been recently rewritten to improve speed and to incorporate a C++ codebase, and as such the features of the package are somewhat in flux. Be sure to update the package frequently.
    
14. 
    
15. To install lme4, we just run:
    
16. 
    
17. # Main version
    
18. install.packages("lme4")
    
19. 
    
20. # Or to install the dev version
    
21. library(devtools)
    
22. install_github("lme4", user = "lme4")

复制代码

扫码加我 拉你入群

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

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

分享0 收藏2 回帖

关键词：Multilevel Modeling getting Started Level themselves structure account complex nested

1. Read in the data
   
2. Multilevel models are appropriate for a particular kind of data structure where units are nested within groups (generally 5+ groups) and where we want to model the group structure of the data. For our introductory example we will start with a simple example from the lme4 documentation and explain what the model is doing. We will use data from Jon Starkweather at the University of North Texas. Visit the excellent tutorial available here for more.
   
3. 
   
4. library(lme4)  # load library
   
5. library(arm)  # convenience functions for regression in R
   
6. lmm.data <- read.table("http://www.unt.edu/rss/class/Jon/R_SC/Module9/lmm.data.txt",
   
7.     header = TRUE, sep = ",", na.strings = "NA", dec = ".", strip.white = TRUE)
   
8. # summary(lmm.data)
   
9. head(lmm.data)
   
10. ##   id extro  open agree social class school
    
11. ## 1  1 63.69 43.43 38.03  75.06     d     IV
    
12. ## 2  2 69.48 46.87 31.49  98.13     a     VI
    
13. ## 3  3 79.74 32.27 40.21 116.34     d     VI
    
14. ## 4  4 62.97 44.41 30.51  90.47     c     IV
    
15. ## 5  5 64.25 36.86 37.44  98.52     d     IV
    
16. ## 6  6 50.97 46.26 38.83  75.22     d      I
    
17. Here we have data on the extroversion of subjects nested within classes and within schools.

复制代码

1. Fit the Non-Multilevel Models
   
2. Let's start by fitting a simple OLS regression of measures of openness, agreeableness, and socialability on extroversion. Here we use the display function in the excellent arm package for abbreviated output. Other options include stargazer for LaTeX typeset tables, xtable, or the ascii package for more flexible plain text output options.
   
3. 
   
4. OLSexamp <- lm(extro ~ open + agree + social, data = lmm.data)
   
5. display(OLSexamp)
   
6. ## lm(formula = extro ~ open + agree + social, data = lmm.data)
   
7. ##             coef.est coef.se
   
8. ## (Intercept) 57.84     3.15  
   
9. ## open         0.02     0.05  
   
10. ## agree        0.03     0.05  
    
11. ## social       0.01     0.02  
    
12. ## ---
    
13. ## n = 1200, k = 4
    
14. ## residual sd = 9.34, R-Squared = 0.00
    
15. So far this model does not fit very well at all. The R model interface is quite a simple one with the dependent variable being specified first, followed by the ~ symbol. The righ hand side, predictor variables, are each named. Addition signs indicate that these are modeled as additive effects. Finally, we specify that datframe on which to calculate the model. Here we use the lm function to perform OLS regression, but there are many other options in R.
    
16. 
    
17. If we want to extract measures such as the AIC, we may prefer to fit a generalized linear model with glm which produces a model fit through maximum likelihood estimation. Note that the model formula specification is the same.
    
18. 
    
19. MLexamp <- glm(extro ~ open + agree + social, data = lmm.data)
    
20. display(MLexamp)
    
21. ## glm(formula = extro ~ open + agree + social, data = lmm.data)
    
22. ##             coef.est coef.se
    
23. ## (Intercept) 57.84     3.15  
    
24. ## open         0.02     0.05  
    
25. ## agree        0.03     0.05  
    
26. ## social       0.01     0.02  
    
27. ## ---
    
28. ##   n = 1200, k = 4
    
29. ##   residual deviance = 104378.2, null deviance = 104432.7 (difference = 54.5)
    
30. ##   overdispersion parameter = 87.3
    
31. ##   residual sd is sqrt(overdispersion) = 9.34
    
32. AIC(MLexamp)
    
33. ## [1] 8774
    
34. This results in a poor model fit. Let's look at a simple varying intercept model now.

复制代码

1. Fit a varying intercept model
   
2. Depending on disciplinary norms, our next step might be to fit a varying intercept model using a grouping variable such as school or classes. Using the glm function and the familiar formula interface, such a fit is easy:
   
3. 
   
4. MLexamp.2 <- glm(extro ~ open + agree + social + class, data = lmm.data)
   
5. display(MLexamp.2)
   
6. ## glm(formula = extro ~ open + agree + social + class, data = lmm.data)
   
7. ##             coef.est coef.se
   
8. ## (Intercept) 56.05     3.09  
   
9. ## open         0.03     0.05  
   
10. ## agree       -0.01     0.05  
    
11. ## social       0.01     0.02  
    
12. ## classb       2.06     0.75  
    
13. ## classc       3.70     0.75  
    
14. ## classd       5.67     0.75  
    
15. ## ---
    
16. ##   n = 1200, k = 7
    
17. ##   residual deviance = 99187.7, null deviance = 104432.7 (difference = 5245.0)
    
18. ##   overdispersion parameter = 83.1
    
19. ##   residual sd is sqrt(overdispersion) = 9.12
    
20. AIC(MLexamp.2)
    
21. ## [1] 8719
    
22. anova(MLexamp, MLexamp.2, test = "F")
    
23. ## Analysis of Deviance Table
    
24. ##
    
25. ## Model 1: extro ~ open + agree + social
    
26. ## Model 2: extro ~ open + agree + social + class
    
27. ##   Resid. Df Resid. Dev Df Deviance    F  Pr(>F)   
    
28. ## 1      1196     104378                             
    
29. ## 2      1193      99188  3     5190 20.8 3.8e-13 ***
    
30. ## ---
    
31. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
32. This is called a fixed-effects specification often. This is simply the case of fitting a separate dummy variable as a predictor for each class. We can see this does not provide much additional model fit. Let's see if school performs any better.
    
33. 
    
34. MLexamp.3 <- glm(extro ~ open + agree + social + school, data = lmm.data)
    
35. display(MLexamp.3)
    
36. ## glm(formula = extro ~ open + agree + social + school, data = lmm.data)
    
37. ##             coef.est coef.se
    
38. ## (Intercept) 45.02     0.92  
    
39. ## open         0.01     0.01  
    
40. ## agree        0.03     0.02  
    
41. ## social       0.00     0.00  
    
42. ## schoolII     7.91     0.27  
    
43. ## schoolIII   12.12     0.27  
    
44. ## schoolIV    16.06     0.27  
    
45. ## schoolV     20.43     0.27  
    
46. ## schoolVI    28.05     0.27  
    
47. ## ---
    
48. ##   n = 1200, k = 9
    
49. ##   residual deviance = 8496.2, null deviance = 104432.7 (difference = 95936.5)
    
50. ##   overdispersion parameter = 7.1
    
51. ##   residual sd is sqrt(overdispersion) = 2.67
    
52. AIC(MLexamp.3)
    
53. ## [1] 5774
    
54. anova(MLexamp, MLexamp.3, test = "F")
    
55. ## Analysis of Deviance Table
    
56. ##
    
57. ## Model 1: extro ~ open + agree + social
    
58. ## Model 2: extro ~ open + agree + social + school
    
59. ##   Resid. Df Resid. Dev Df Deviance    F Pr(>F)   
    
60. ## 1      1196     104378                           
    
61. ## 2      1191       8496  5    95882 2688 <2e-16 ***
    
62. ## ---
    
63. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
64. The school effect greatly improves our model fit. However, how do we interpret these effects?
    
65. 
    
66. table(lmm.data$school, lmm.data$class)
    
67. ##      
    
68. ##        a  b  c  d
    
69. ##   I   50 50 50 50
    
70. ##   II  50 50 50 50
    
71. ##   III 50 50 50 50
    
72. ##   IV  50 50 50 50
    
73. ##   V   50 50 50 50
    
74. ##   VI  50 50 50 50
    
75. Here we can see we have a perfectly balanced design with fifty observations in each combination of class and school (if only data were always so nice!).
    
76. 
    
77. Let's try to model each of these unique cells. To do this, we fit a model and use the : operator to specify the interaction between school and class.
    
78. 
    
79. MLexamp.4 <- glm(extro ~ open + agree + social + school:class, data = lmm.data)
    
80. display(MLexamp.4)
    
81. ## glm(formula = extro ~ open + agree + social + school:class, data = lmm.data)
    
82. ##                  coef.est coef.se
    
83. ## (Intercept)       80.36     0.37
    
84. ## open               0.01     0.00
    
85. ## agree             -0.01     0.01
    
86. ## social             0.00     0.00
    
87. ## schoolI:classa   -40.39     0.20
    
88. ## schoolII:classa  -28.15     0.20
    
89. ## schoolIII:classa -23.58     0.20
    
90. ## schoolIV:classa  -19.76     0.20
    
91. ## schoolV:classa   -15.50     0.20
    
92. ## schoolVI:classa  -10.46     0.20
    
93. ## schoolI:classb   -34.60     0.20
    
94. ## schoolII:classb  -26.76     0.20
    
95. ## schoolIII:classb -22.59     0.20
    
96. ## schoolIV:classb  -18.71     0.20
    
97. ## schoolV:classb   -14.31     0.20
    
98. ## schoolVI:classb   -8.54     0.20
    
99. ## schoolI:classc   -31.86     0.20
    
100. ## schoolII:classc  -25.64     0.20
     
101. ## schoolIII:classc -21.58     0.20
     
102. ## schoolIV:classc  -17.58     0.20
     
103. ## schoolV:classc   -13.38     0.20
     
104. ## schoolVI:classc   -5.58     0.20
     
105. ## schoolI:classd   -30.00     0.20
     
106. ## schoolII:classd  -24.57     0.20
     
107. ## schoolIII:classd -20.64     0.20
     
108. ## schoolIV:classd  -16.60     0.20
     
109. ## schoolV:classd   -12.04     0.20
     
110. ## ---
     
111. ##   n = 1200, k = 27
     
112. ##   residual deviance = 1135.9, null deviance = 104432.7 (difference = 103296.8)
     
113. ##   overdispersion parameter = 1.0
     
114. ##   residual sd is sqrt(overdispersion) = 0.98
     
115. AIC(MLexamp.4)
     
116. ## [1] 3396
     
117. This is very useful, but what if we want to understand both the effect of the school and the effect of the class, as well as the effect of the schools and classes? Unfortunately, this is not easily done with the standard glm.
     
118. 
     
119. MLexamp.5 <- glm(extro ~ open + agree + social + school * class - 1, data = lmm.data)
     
120. display(MLexamp.5)
     
121. ## glm(formula = extro ~ open + agree + social + school * class -
     
122. ##     1, data = lmm.data)
     
123. ##                  coef.est coef.se
     
124. ## open              0.01     0.00  
     
125. ## agree            -0.01     0.01  
     
126. ## social            0.00     0.00  
     
127. ## schoolI          39.96     0.36  
     
128. ## schoolII         52.21     0.36  
     
129. ## schoolIII        56.78     0.36  
     
130. ## schoolIV         60.60     0.36  
     
131. ## schoolV          64.86     0.36  
     
132. ## schoolVI         69.90     0.36  
     
133. ## classb            5.79     0.20  
     
134. ## classc            8.53     0.20  
     
135. ## classd           10.39     0.20  
     
136. ## schoolII:classb  -4.40     0.28  
     
137. ## schoolIII:classb -4.80     0.28  
     
138. ## schoolIV:classb  -4.74     0.28  
     
139. ## schoolV:classb   -4.60     0.28  
     
140. ## schoolVI:classb  -3.87     0.28  
     
141. ## schoolII:classc  -6.02     0.28  
     
142. ## schoolIII:classc -6.54     0.28  
     
143. ## schoolIV:classc  -6.36     0.28  
     
144. ## schoolV:classc   -6.41     0.28  
     
145. ## schoolVI:classc  -3.65     0.28  
     
146. ## schoolII:classd  -6.81     0.28  
     
147. ## schoolIII:classd -7.45     0.28  
     
148. ## schoolIV:classd  -7.24     0.28  
     
149. ## schoolV:classd   -6.93     0.28  
     
150. ## schoolVI:classd   0.06     0.28  
     
151. ## ---
     
152. ##   n = 1200, k = 27
     
153. ##   residual deviance = 1135.9, null deviance = 4463029.9 (difference = 4461894.0)
     
154. ##   overdispersion parameter = 1.0
     
155. ##   residual sd is sqrt(overdispersion) = 0.98
     
156. AIC(MLexamp.5)
     
157. ## [1] 3396

复制代码

1. Exploring Random Slopes
   
2. Another alternative is to fit a separate model for each of the school and class combinations. If we believe the relationsihp between our variables may be highly dependent on the school and class combination, we can simply fit a series of models and explore the parameter variation among them:
   
3. 
   
4. require(plyr)
   
5. 
   
6. modellist <- dlply(lmm.data, .(school, class), function(x) glm(extro ~ open +
   
7.     agree + social, data = x))
   
8. display(modellist[[1]])
   
9. ## glm(formula = extro ~ open + agree + social, data = x)
   
10. ##             coef.est coef.se
    
11. ## (Intercept) 35.87     5.90  
    
12. ## open         0.05     0.09  
    
13. ## agree        0.02     0.10  
    
14. ## social       0.01     0.03  
    
15. ## ---
    
16. ##   n = 50, k = 4
    
17. ##   residual deviance = 500.2, null deviance = 506.2 (difference = 5.9)
    
18. ##   overdispersion parameter = 10.9
    
19. ##   residual sd is sqrt(overdispersion) = 3.30
    
20. display(modellist[[2]])
    
21. ## glm(formula = extro ~ open + agree + social, data = x)
    
22. ##             coef.est coef.se
    
23. ## (Intercept) 47.96     2.16  
    
24. ## open        -0.01     0.03  
    
25. ## agree       -0.03     0.03  
    
26. ## social      -0.01     0.01  
    
27. ## ---
    
28. ##   n = 50, k = 4
    
29. ##   residual deviance = 47.9, null deviance = 49.1 (difference = 1.2)
    
30. ##   overdispersion parameter = 1.0
    
31. ##   residual sd is sqrt(overdispersion) = 1.02
    
32. We will discuss this strategy in more depth in future tutorials including how to performan inference on the list of models produced in this command.

复制代码

1. Fit a varying intercept model with lmer
   
2. Enter lme4. While all of the above techniques are valid approaches to this problem, they are not necessarily the best approach when we are interested explicitly in variation among and by groups. This is where a mixed-effect modeling framework is useful. Now we use the lmer function with the familiar formula interface, but now group level variables are specified using a special syntax: (1|school) tells lmer to fit a linear model with a varying-intercept group effect using the variable school.
   
3. 
   
4. MLexamp.6 <- lmer(extro ~ open + agree + social + (1 | school), data = lmm.data)
   
5. display(MLexamp.6)
   
6. ## lmer(formula = extro ~ open + agree + social + (1 | school),
   
7. ##     data = lmm.data)
   
8. ##             coef.est coef.se
   
9. ## (Intercept) 59.12     4.10  
   
10. ## open         0.01     0.01  
    
11. ## agree        0.03     0.02  
    
12. ## social       0.00     0.00  
    
13. ##
    
14. ## Error terms:
    
15. ##  Groups   Name        Std.Dev.
    
16. ##  school   (Intercept) 9.79   
    
17. ##  Residual             2.67   
    
18. ## ---
    
19. ## number of obs: 1200, groups: school, 6
    
20. ## AIC = 5836.1, DIC = 5789
    
21. ## deviance = 5806.5
    
22. We can fit multiple group effects with multiple group effect terms.
    
23. 
    
24. MLexamp.7 <- lmer(extro ~ open + agree + social + (1 | school) + (1 | class),
    
25.     data = lmm.data)
    
26. display(MLexamp.7)
    
27. ## lmer(formula = extro ~ open + agree + social + (1 | school) +
    
28. ##     (1 | class), data = lmm.data)
    
29. ##             coef.est coef.se
    
30. ## (Intercept) 60.20     4.21  
    
31. ## open         0.01     0.01  
    
32. ## agree       -0.01     0.01  
    
33. ## social       0.00     0.00  
    
34. ##
    
35. ## Error terms:
    
36. ##  Groups   Name        Std.Dev.
    
37. ##  school   (Intercept) 9.79   
    
38. ##  class    (Intercept) 2.41   
    
39. ##  Residual             1.67   
    
40. ## ---
    
41. ## number of obs: 1200, groups: school, 6; class, 4
    
42. ## AIC = 4737.9, DIC = 4683
    
43. ## deviance = 4703.6
    
44. And finally, we can fit nested group effect terms through the following syntax:
    
45. 
    
46. MLexamp.8 <- lmer(extro ~ open + agree + social + (1 | school/class), data = lmm.data)
    
47. display(MLexamp.8)
    
48. ## lmer(formula = extro ~ open + agree + social + (1 | school/class),
    
49. ##     data = lmm.data)
    
50. ##             coef.est coef.se
    
51. ## (Intercept) 60.24     4.01  
    
52. ## open         0.01     0.00  
    
53. ## agree       -0.01     0.01  
    
54. ## social       0.00     0.00  
    
55. ##
    
56. ## Error terms:
    
57. ##  Groups       Name        Std.Dev.
    
58. ##  class:school (Intercept) 2.86   
    
59. ##  school       (Intercept) 9.69   
    
60. ##  Residual                 0.98   
    
61. ## ---
    
62. ## number of obs: 1200, groups: class:school, 24; school, 6
    
63. ## AIC = 3568.6, DIC = 3508
    
64. ## deviance = 3531.1
    
65. 
    
66. Here the (1|school/class) says that we want to fit a mixed effect term for varying intercepts 1| by schools, and for classes that are nested within schools.

复制代码

1. Fit a varying slope model with lmer
   
2. But, what if we want to explore the effect of different student level indicators as they vary across classrooms. Instead of fitting unique models by school (or school/class) we can fit a varying slope model. Here we modify our random effect term to include variables before the grouping terms: (1 + open|school/class) tells R to fit a varying slope and varying intercept model for schools and classes nested within schools, and to allow the slope of the open variable to vary by school.
   
3. 
   
4. MLexamp.9 <- lmer(extro ~ open + agree + social + (1 + open | school/class),
   
5.     data = lmm.data)
   
6. display(MLexamp.9)
   
7. ## lmer(formula = extro ~ open + agree + social + (1 + open | school/class),
   
8. ##     data = lmm.data)
   
9. ##             coef.est coef.se
   
10. ## (Intercept) 60.26     3.93  
    
11. ## open         0.01     0.01  
    
12. ## agree       -0.01     0.01  
    
13. ## social       0.00     0.00  
    
14. ##
    
15. ## Error terms:
    
16. ##  Groups       Name        Std.Dev. Corr
    
17. ##  class:school (Intercept) 2.62         
    
18. ##               open        0.01     1.00
    
19. ##  school       (Intercept) 9.51         
    
20. ##               open        0.00     1.00
    
21. ##  Residual                 0.98         
    
22. ## ---
    
23. ## number of obs: 1200, groups: class:school, 24; school, 6
    
24. ## AIC = 3574.7, DIC = 3506
    
25. ## deviance = 3529.3

复制代码

[code]Conclusion
Fitting mixed effect models and exploring group level variation is very easy within the R language and ecosystem. In future tutorials we will explore comparing across models, doing inference with mixed-effect models, and creating graphical representations of mixed effect models to understand their effects.

Appendix
print(sessionInfo(), locale = FALSE)
## R version 3.0.1 (2013-05-16)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
##
## other attached packages:
## [1] plyr_1.8        arm_1.6-10      MASS_7.3-29     lme4_1.0-5     
## [5] Matrix_1.1-0    lattice_0.20-24 knitr_1.5      
##
## loaded via a namespace (and not attached):
##  [1] abind_1.4-0    coda_0.16-1    evaluate_0.5.1 formatR_0.10  
##  [5] grid_3.0.1     minqa_1.2.1    nlme_3.1-113   splines_3.0.1
##  [9] stringr_0.6.2  tools_3.0.1
[1] Examples include Gelman and Hill, Gelman et al. 2013, etc.

thank you