[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]

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]]