Tutorials: Causal Models for Mediation Analyses

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

Zheng, C., Atkins, D. C., & Zhou, A. (2012). Causal models for mediation analyses: An introduction to structural mean models. Manuscript submitted for publication.

本帖隐藏的内容

• [Draft]
• [Data]
• [R code]
• [Simulation R code]
  

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：mediation Tutorials Tutorial Analyses analyse

1. ### R code to accompany (Updated, 29 August 2012):
   
2. #
   
3. # Zheng, C., Atkins, D. C., & Zhou, A. (2012). Causal models for mediation analyses: An introduction to structural mean models.  Manuscript submitted for publication.
   
4. #
   
5. ### Analyses demonstrating how to fit the RPM model described in
   
6. ### the paper are shown below.  Data originally from:
   
7. #
   
8. # Whiteside, U., Atkins, D. C., Kleiber, B. V., Neighbors, C., Witkiewitz, K., & Larimer, M. E. (2011).  DBT skills plus personalized normative feedback: Results of a randomized clinical trial.  Manuscript submitted for publication.
   
9. 
   
10. ### Import data
    
11. newdata <- read.csv("http://depts.washington.edu/cshrb/newweb/stats%20documents/Jan2012Mediation.csv",
    
12.                     header = TRUE)
    
13. head(newdata)
    
14. summary(newdata)
    
15. 
    
16. ### NOTE: variables with prefix "b" were assessed at baseline,
    
17. ###                  and variables with prefix "o" were assessed at one month
    
18. ###                  post intervention
    
19. #
    
20. ### NOTE: original study had 3 treatment groups, but our mediation
    
21. ###                  analyses focus only on DBT-BASICS and Control
    
22. #
    
23. ### NOTE: some individuals had missing data at one month and are not
    
24. ###                  included here
    
25. 
    
26. # Original variables will be re-named to generic names using the
    
27. # following conventions:
    
28. #   
    
29. # m = mediator
    
30. # y = outcome
    
31. # r = intervention indicator
    
32. # x.y = covariate for outcome
    
33. # x.m = covariate for mediator
    
34. #
    
35. ### Hypothesized mediator: One month DERS
    
36. m <- newdata$o_ders
    
37. 
    
38. ### Convert treatment to binary indicator
    
39. r <- 1 - as.numeric(newdata$condition == "Control")
    
40. 
    
41. ### Outcome: One month BDI
    
42. y <- newdata$o_bdi  
    
43. 
    
44. ### Baseline covariate (centered around mean)
    
45. x.m <- newdata$b_bai - mean(newdata$b_bai)
    
46. x.y <- newdata$b_bai - mean(newdata$b_bai)
    
47. 
    
48. ### NOTE: In the current usage, x.m and x.y (covariates for mediator model
    
49. ###             and outcome model) are identical; however, they do not have to be,
    
50. ###                    and so we show them as two separate terms here
    
51. #
    
52. ### Pull data together into data.frame
    
53. df <- data.frame(y = y, m = m, r = r, x.m = x.m, x.y = x.y)
    
54. head(df)
    
55. 
    
56. ### The rank preserving model (RPM), which is a particular way of
    
57. ### estimating a structural mean model, is fit using estimating
    
58. ### equations.  The following R code gives an example of how to
    
59. ### fit the RPM to the college student drinking data and directly
    
60. ### follows the equations in the technical appendix of the ms.
    
61. #
    
62. ### Fit model for m to obtain weight W_i(X) in w1, w2
    
63. #
    
64. ### First, fit linear model of mediator predicted from treatment by
    
65. ### baseline interaction
    
66. modelm <- lm(m ~ r*x.m, data = df)
    
67. summary(modelm)
    
68. 
    
69. ### NOTE: Interaction of treatment and baseline covariate (here, BAI),
    
70. ###       is not very strong (p = 0.19), which will impact the variance
    
71. ###       of the RPM estimate.
    
72. 
    
73. ### The following code generates the predicted value of the mediator
    
74. ### given baseline covariates for each level of the treatment, based
    
75. ### on Equation 8 in the Technical Appendix
    
76. tmp.data <- expand.grid(x.m = df$x.m,
    
77.                         r = c(0,1))
    
78. tmp.data$pred <- predict(modelm, tmp.data, type = "response")
    
79. head(tmp.data)
    
80. 
    
81. ### Now, get treatment difference on mediator (conditional on x.m)
    
82. w <- tmp.data[tmp.data$r == 1, "pred"] -
    
83.      tmp.data[tmp.data$r == 0, "pred"]
    
84. 
    
85. ### "weight 1" is the difference between treatment and average treatment,
    
86. ### that is, r^tilde - r^bar in equation 7:
    
87. w1 <- df$r - mean(df$r)
    
88. 
    
89. ### "weight 2" is this difference multiplied by w (ie, treatment difference in
    
90. ### mediator conditional on x.m; see equation 8)
    
91. w2 <- (df$r - mean(df$r))*w
    
92. 
    
93. ### For Equation 7 in the Appendix, we do not include additional
    
94. ### baseline covariates.  See below for estimation using Equation 9.1
    
95. ### that does include baseline covariates.
    
96. #
    
97. ### To more clearly connect the R code with the Technical Appendix
    
98. ### we include a representation of Equation 7 here (as best we can
    
99. ### given unformatted text); plugging in Y^00(b, c') = Y - c'r - bm:
    
100. #
     
101. ### Sigma_i (r-Er) W(x) (Y-c'r-bm) = 0      [7]
     
102. 
     
103. ### We have already estimated the first two parts, which we combine:
     
104. WW <- cbind(w1,w2)
     
105. head(WW)
     
106. 
     
107. ### Relating this to Equation 7 in Technical Appendix:
     
108. #
     
109. ### WW = (r-Er) W(x)
     
110. 
     
111. ### Re-arranging Eq. 7 above:
     
112. #
     
113. ### (r-Er) W(x) Y = (r-Er) W(x) (r,m) %*% (c',b)
     
114. #
     
115. ### Combine observed treatment and mediator, where (r,m) = MM
     
116. MM <- cbind(df$r, df$m)
     
117. 
     
118. ### To solve for (c',b), we take final steps:
     
119. #
     
120. ### YY = MM %*% (c',b) and
     
121. #
     
122. ### theta = (c',b) = MM^{-1} YY.
     
123. XX <- crossprod(WW, MM)
     
124. YY <- crossprod(WW, df$y)
     
125. theta <- tcrossprod(solve(XX), t(YY))
     
126. theta
     
127. 
     
128. ### Calculate y00 and then get sandwich variance estimation
     
129. y0 <- df$y - tcrossprod(MM, t(theta))
     
130. V <- tcrossprod(tcrossprod(solve(XX), WW * as.vector(y0)))
     
131. B3 <- solve(solve(crossprod(WW * as.vector(y0)))[1:2,1:2])
     
132. 
     
133. ### "Sandwich" variance estimates
     
134. V3 <- solve(XX) %*% B3 %*% solve(t(XX))
     
135. V3
     
136. 
     
137. ### Calculate 95% CI
     
138. #
     
139. ### Extract SE from variance-covariance matrix
     
140. se <- sqrt(diag(V3))
     
141. res3 <- cbind(b = theta[1:2],
     
142.                           lo = theta[1:2] + qnorm(0.025)*se,
     
143.                           hi = theta[1:2] + qnorm(0.975)*se)
     
144. round(res3, 2)
     
145. 
     
146. ### Note: The results reported in the ms use baseline covariates; thus,
     
147. ### see the following section.
     
148. 
     
149. ############# Augmented with Baseline Covariates #########
     
150. #
     
151. ### As noted in the Technical Appendix, it is possible to reduce
     
152. ### the variance of the estimated parameters (i.e., increase
     
153. ### efficiency), by including baseline covariates as shown in
     
154. ### Equation 9.1
     
155. #
     
156. ### Based on equations 9.1 and 10.1, we augment the estimating
     
157. ### equations with an intercept and baseline covariate:
     
158. w4 <- 1
     
159. w5 <- df$x.y
     
160. WW.2 <- cbind(w1,w2,w4,w5)
     
161. head(WW.2)
     
162. 
     
163. ### Note that WW.2 is the matrix shown in the explicit form of the solution,
     
164. ### directly following presentation of equation 10.1 in the ms.  Moreover,
     
165. ### MM.2 (below) is the middle matrix:
     
166. MM.2 <- cbind(df$r, df$m, 1, df$x.y)
     
167. 
     
168. ### We can then estimate coefficients via:
     
169. XX.2 <- crossprod(WW.2, MM.2)
     
170. YY.2 <- crossprod(WW.2, df$y)
     
171. theta.2 <- tcrossprod(solve(XX.2), t(YY.2))
     
172. theta.2
     
173. 
     
174. ### Calculate y00 and then get sandwich variance estimation
     
175. y0.2 <- df$y - tcrossprod(MM.2, t(theta.2))
     
176. V.2 <- tcrossprod(tcrossprod(solve(XX.2), WW.2 * as.vector(y0.2)))
     
177. A.2 <- XX.2[1:2, 1:2]
     
178. B3.2 <- solve(solve(crossprod(WW.2 * as.vector(y0.2)))[1:2,1:2])
     
179. 
     
180. ### "Sandwich" variance estimates
     
181. V3.2 <- solve(A.2) %*% B3.2 %*% solve(t(A.2))
     
182. V3.2
     
183. 
     
184. ### Calculate 95% CI
     
185. #
     
186. ### Extract SE from variance-covariance matrix
     
187. se2 <- sqrt(diag(V3.2))
     
188. res3.2 <- cbind(b = theta.2[1:2],
     
189.               lo = theta.2[1:2] + qnorm(0.025)*se2,
     
190.               hi = theta.2[1:2] + qnorm(0.975)*se2)
     
191. round(res3.2, 2)
     
192. round(res3, 2)
     
193. 
     
194. ### Note that we also see that including the baseline covariate increased the efficiency of the estimates.

复制代码

1. ### Import original data and use the treatment as well as
   
2. ### baseline BAI as z and x
   
3. 
   
4. newdata <- read.csv("Mediation.csv", header = TRUE)
   
5. 
   
6. ### Generate simulation data, alphau and betau measure
   
7. ### violation of sequential ignorability and gamma measures
   
8. ### interaction
   
9. 
   
10. gendata <- function(thetaz, thetam, betax, betau, beta0,
    
11.                                         alphaz, alphax, alpha0, gamma, alphau) {
    
12. 
    
13. z <- 1 - as.numeric(newdata$condition=="Control")
    
14. x.m <- newdata$b_bai
    
15. x.y <- newdata$b_bai
    
16. n <- length(z)
    
17. u <- rnorm(n)
    
18. m <- alphaz*z + alphax*x.m + gamma*x.m*z + alphau*u +
    
19.           alpha0 + rnorm(n, sd = 0.2)
    
20. 
    
21. y <- thetaz*z + thetam*m + betax*x.y + betau*u + beta0 +
    
22.                 rnorm(n, sd = 0.5)
    
23. 
    
24. list(z, m, y, x.m, x.y)
    
25. }
    
26. 
    
27. library(sandwich) # for robust SE
    
28. 
    
29. ### RPM Estimation function, return with estimates and CI
    
30. ### from RPM and OLS
    
31. 
    
32. mymed <- function(z, m, y, x.m, x.y){
    
33. 
    
34. data <- na.omit(data.frame(z = z, m = m, y = y,
    
35.                                 x.m = x.m, x.y = x.y))
    
36. 
    
37. modelm <- lm(m ~ z*x.m, data = data)
    
38. data1 <- data0 <- data
    
39. data1$z <- 1
    
40. data0$z <- 0
    
41. 
    
42. w <- as.matrix(predict(modelm, data1, type = "response")-predict(modelm, data0, type = "response"))
    
43. w1 <- (data$z-mean(data$z))
    
44. w2 <- (data$z-mean(data$z))*w
    
45. 
    
46. # w4,  w5 is for the augmented estimating equations
    
47. 
    
48. w4 <- 1
    
49. w5 <- data$x.y
    
50. WW <- cbind(w1, w2, w4, w5)
    
51. 
    
52. # solve the equation by its explicit expression
    
53. 
    
54. MM <- cbind(data$z, data$m, 1, data$x.y)
    
55. XX <- crossprod(WW, MM)
    
56. YY <- crossprod(WW, data$y)
    
57. theta <- tcrossprod(solve(XX), t(YY))
    
58. 
    
59. # calculate y00 and then get sandwich variance estimation
    
60. 
    
61. y0 <- data$y-tcrossprod(MM, t(theta))
    
62. V <- tcrossprod(tcrossprod(solve(XX), WW*as.vector(y0)))
    
63. B3 <- solve(solve(crossprod(WW*as.vector(y0))))
    
64. A <- XX
    
65. V3 <- (solve(A)%*%B3%*%solve(t(A)))[1:2, 1:2]
    
66. 
    
67. # Calculate 95% CI
    
68. res3 <- cbind(theta[1:2], theta[1:2] + qnorm(0.025)*sqrt(diag(V3)), theta[1:2] + qnorm(0.975)*sqrt(diag(V3)))
    
69. 
    
70. #Indirect effect and Mediated effect via OLS
    
71. res <- lm(y~z + m + x.y, data = data)
    
72. 
    
73. #Compute 95% CI for indirect and mediated effect
    
74. rbind(res3, cbind(coef(res), coef(res)-1.96*sqrt(diag(vcovHC(res))), coef(res) + 1.96*sqrt(diag(vcovHC(res))))[2:3, ])
    
75. 
    
76. }
    
77. 
    
78. # one time simulation
    
79. onesim <- function(betau = betau, gamma = gamma){
    
80. 
    
81. thetaz <- (-1.55)
    
82. thetam <- 0.37
    
83. betax <- 0.1
    
84. beta0 <- 0
    
85. alphaz <- -2.2
    
86. alphax <- 0.9
    
87. alpha0 <- 0
    
88. alphau <- 1
    
89. 
    
90. simdata <- gendata(thetaz, thetam, betax, betau, beta0, alphaz, alphax, alpha0, gamma, alphau)
    
91. 
    
92. z <- simdata[[1]]
    
93. m <- simdata[[2]]
    
94. y <- simdata[[3]]
    
95. x.m <- simdata[[4]]
    
96. x.y <- simdata[[5]]
    
97. 
    
98. result <- mymed(z, m, y, x.m, x.y)
    
99. 
    
100. bias <- result[,1]-rep(c(thetaz, thetam), 2)
     
101. 
     
102. cr <- as.numeric(((result[,2]-rep(c(thetaz, thetam), 2))*(result[,3]-rep(c(thetaz, thetam), 2)))<0)
     
103. 
     
104. c(bias, cr)
     
105. 
     
106. }
     
107. 
     
108. # Simulation
     
109. M <- 1000
     
110. 
     
111. myfun <- function(betau = betau, gamma = gamma) {
     
112. 
     
113.         res <- replicate(M, onesim(betau, gamma), simplify = TRUE)
     
114. 
     
115.         # calculate bias and coverage rate
     
116.         a1 <- apply(res, 1, mean)
     
117. 
     
118.         # calculate sd
     
119.         a2 <- apply(res, 1, sd)
     
120. 
     
121.         out <- cbind(a1[1:4], a2[1:4], a1[5:8])
     
122.         colnames(out) <- c("Bias","SD","CR")
     
123.         rownames(out) <- c("RPM.Z","RPM.M","OLS.Z","OLS.M")
     
124.         round(out, 2)
     
125. }
     
126. 
     
127. # Sequential ignorability fails and weak interaction
     
128. set.seed(111)
     
129. myfun(2, -0.1)
     
130. ### NOTE: ~17 sec for 1,000 simulations on a current MacBook Pro
     
131. 
     
132. # Sequential ignorability fail and strong interaction
     
133. set.seed(111)
     
134. myfun(2, -0.5)
     
135. 
     
136. # Sequential ignorability hold and weak interaction
     
137. set.seed(111)
     
138. myfun(0, -0.1)
     
139. 
     
140. # Sequential ignorability hold and strong interaction
     
141. set.seed(111)
     
142. myfun(0, -0.5)

复制代码

已经阅读论坛新手操作内容已经阅读论坛新手操作内容已经阅读论坛新手操作内容已经阅读论坛新手操作内容

感谢楼主提供好资料！

thanks