[R] mediation analysis in R using the method described by Lange T., et al. in the paper "A simple unified approach for estimating natural direct and direct effects
Ruzan Udumyan
ruzudumyan at gmail.com
Sun Dec 21 09:03:04 CET 2014
Dear All,
I need to do a mediation analysis with survival data (using Cox model).
- The independent variable is categorical (with 3 levels, coded as 1, 2,
3).
- Mediator variable is a 1-10 score variable (derived from a continuous
variable).
I must confess, I am a Stata user. However, could not find an appropriate
command in Stata for mediation analysis with Cox model, and was very happy
to find a way of doing it in R.
I am running into problems though, which I do not know how to solve. More
specifically, I get an error message from RStudio after running the syntax
mentioned below.
I am kindly asking you to help me fix the problem.
With many thanks and my best regards,
Ruzan
Error in model.frame.default(formula = Surv(time, event) ~ factor(stress) + :
variable lengths differ (found for 'factor(stress)')
In addition: Warning messages:1: In checkwz(wz, M = M, trace = trace,
wzepsilon = control$wzepsilon) :
31 elements replaced by 1.819e-122: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
187 elements replaced by 1.819e-123: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
1492 elements replaced by 1.819e-124: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
2081 elements replaced by 1.819e-125: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
2081 elements replaced by 1.819e-126: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
4783 elements replaced by 1.819e-127: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
8309 elements replaced by 1.819e-128: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
8642 elements replaced by 1.819e-129: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
8642 elements replaced by 1.819e-1210: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
8642 elements replaced by 1.819e-12
>
Here is the syntax:
doEffectDecomp = function(myData)
{
#step 1
myData$stresstemp <- myData$stress
library(VGAM)
fitM <-vglm(phys1 ~ factor(stresstemp)+ C, data = myData, family =
multinomial()) # C stands for confounders
#summary(fitM)
#Step 2. Next we construct new variables id and star
N <- nrow(myData)
myData$id <- 1:N # construct id variable
levelsOfstress <- unique(myData$stress)
myData1 <- myData
myData2 <- myData
myData3 <- myData
myData1$star <- levelsOfstress[1]
myData2$star <- levelsOfstress[2]
myData3$star <- levelsOfstress[3]
newMyData <- rbind(myData1, myData2, myData3)
#Step 3. compute weights
newMyData$stresstemp <- newMyData$stress
tempDir <- as.matrix(predict(fitM,type = "response",
newdata=newMyData))[cbind(1:(3*N),newMyData$phys1)]
newMyData$stresstemp <- newMyData$star
tempIndir <- as.matrix(predict(fitM,type = "response",
newdata=newMyData))[cbind(1:(3*N),newMyData$phys1)]
newMyData$weightM <- tempIndir/tempDir
#hist(newMyData$weightM)
#Step 4. Finally the MSM model for direct and indirect effects can be
estimated.
# weighted Cox model
library(survival)
cox <- coxph(Surv(time,event) ~factor(stress) + factor(star) + C,
method="efron", data=newMyData, weights=newMyData$weightM)
summary(cox)
# return value: TE; DE; IE
# Return value: Estimates for total, direct, indirect effect
TE = exp(sum(coef(cox)[c('stressTRUE', 'starTRUE')])) # I am not
sure this is correct for the categorical exposure with 3 levels
DE = exp(unname(coef(cox)['stressTRUE']))
IE = exp(sum(coef(cox)['starTRUE']))
PM = log(IE) / log(TE)
return(c(exp(coef(cox)), TE=TE, DE=DE, IE=IE, PM=PM))
}
doEffectDecomp(myData)
[[alternative HTML version deleted]]
More information about the R-help
mailing list