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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
