# Causal Mediation Analysis of Survival Outcome with Multiple Mediators

Background: Mediation analyses have been a popular approach to investigate the effect of an exposure on an outcome through a mediator. Mediation models with multiple mediators have been proposed for continuous and dichotomous outcomes. However, development of multimediator models for survival outcomes is still limited.

Methods: We present methods for multimediator analyses using three survival models: Aalen additive hazard models, Cox proportional hazard models, and semiparametric probit models. Effects through mediators can be characterized by path-specific effects, for which definitions and identifiability assumptions are provided. We derive closed-form expressions for path-specific effects for the three models, which are intuitively interpreted using a causal diagram.

Results: Mediation analyses using Cox models under the rare-outcome assumption and Aalen additive hazard models consider effects on log hazard ratio and hazard difference, respectively; analyses using semiparametric probit models consider effects on difference in transformed survival time and survival probability. The three models were applied to a hepatitis study where we investigated effects of hepatitis C on liver cancer incidence mediated through baseline and/or follow-up hepatitis B viral load. The three methods show consistent results on respective effect scales, which suggest an adverse estimated effect of hepatitis C on liver cancer not mediated through hepatitis B, and a protective estimated effect mediated through the baseline (and possibly follow-up) of hepatitis B viral load.

Conclusions: Causal mediation analyses of survival outcome with multiple mediators are developed for additive hazard and proportional hazard and probit models with utility demonstrated in a hepatitis study.

From the ^{a}Institute of Statistical Science, Academia Sinica, Taipei, Taiwan; ^{b}Departments of Epidemiology and Biostatistics, Brown University, Providence, RI; and ^{c}Genomics Research Center, Academia Sinica, Taipei, Taiwan.

Submitted May 4 2016; accepted March 13 2017.

Financially supported by NIH/NCI 5R03CA 182937-02 and NIH/NIA 1R01AG048825-01 and Taiwan MOST 105-2118-M-001-014-MY3.

The authors report no conflicts of interest.

Computation code and data are available at http://www.stat.sinica.edu.tw/ythuang/M3S-Cox%20and%20M3S-Aalen.zip and http://www.stat.sinica.edu.tw/ythuang/M3S-Probit.zip.

Supplemental digital content is available through direct URL citations in the HTML and PDF versions of this article (www.epidem.com).

Correspondence: Yen-Tsung Huang, Institute of Statistical Science, Academia Sinica, 128 Academia Road, Section 2, Nankang, Taipei 11529, Taiwan. E-mail: ythuang@stat.sinica.edu.tw.

Mediation analyses first proposed in psychology literature^{1},^{2} have become a very popular approach and been applied to a wide range of epidemiologic studies.^{3–7} Mediation models characterize relationships of an exposure, a mediator, and an outcome. Specifically, in mediation analyses, the effect of the exposure on the outcome can be decomposed into natural direct effect, the effect of the exposure on the outcome not mediated by the mediator, and natural indirect effect, the effect of the exposure on the outcome mediated through the mediator. By employing the counterfactual framework,^{8} causal mediation models can be formulated as a graphic model illustrated using a directed acyclic graph (DAG),^{9} and causal assumptions for effect identifiability have been carefully studied.^{10} Built upon the framework of causal inference, the methodology of mediation analyses has been generalized from continuous outcomes to dichotomous outcomes^{11},^{12} and time-to-event survival outcomes.^{13},^{14} The identification of natural direct and indirect effect in mediation analyses with an exposure *S*, a mediator *M*, and a survival outcome *Y* requires a set of assumptions.^{14} The assumptions can be stated as that, conditional on measured confounders ** X**, there is (1) no confounding for the relationship of

*S*and

*Y*, (2) no confounding for the relationship of

*M*and

*Y*, conditional on

*S*, (3) no confounding for the

*S*–

*M*association, and (4) no confounder for the

*M*–

*Y*relationship caused by

*S*; additional standard assumptions include positivity and consistency which are discussed in eAppendix A1 (http://links.lww.com/EDE/B185).

This article is motivated by a hepatitis study where mediation by hepatitis B and C in relation to liver cancer was investigated.^{15} The study was conducted in a community-based prospective cohort study in Taiwan, in which viral load of hepatitis B (HBV) and C (HCV) viruses was measured at the baseline, and the incidence of liver cancer was recorded prospectively.^{16} Based on the existing scientific evidence,^{17–23} a mediation model was proposed with the exposure being hepatitis C viral load, the mediator being hepatitis B viral load, and the outcome being the incidence of liver cancer. In addition, hepatitis B viral load was also measured during the follow-up. We extend the above mediation model to include exposure, HCV viral load (*S*), the outcome, the liver cancer incidence (*Y*) and two mediators, one being the baseline HBV viral load (*M*_{1}) and the other being the follow-up HBV viral load (*M*_{2}; Figure 1). Note that if one is interested in the mediator model with only one mediator *M*_{2}, then the existence of *M*_{1} violates the aforementioned assumption (4). To address the issue, we study three effects: (1) the effect of HCV on the liver cancer incidence mediated through the baseline HBV viral load and possibly through the follow-up HBV viral load (the black path), (2) the effect of HCV on the liver cancer incidence mediated through the follow-up HBV viral load, but not through the baseline HBV viral load (the dark gray path), and (3) the effect of HCV on the liver cancer incidence not mediated by HBV regardless through the baseline or follow-up viral load (the light gray path).

To identify the above effects, one has to conduct mediation analyses for a model with a survival outcome and multiple mediators. The three different effects have been termed as path-specific effects, which characterize mediation effects for various pathways through different mediators.^{24} Regression- and weighting-based approaches have been proposed to estimate path-specific effects.^{25} The path-specific effect approach has also been proposed as a method to adjust for exposure-induced confounding for the mediator–outcome association,^{26} similar to our problem in the hepatitis study. Causal mediation models have been generalized to incorporate mixed variable types such as a combination of continuous and dichotomous mediators^{27} or a set of high-dimensional continuous mediators.^{28} However, these methods focus on noncensored outcomes. There has been some work on mediation analyses for survival outcomes. Lange and Hansen^{13} introduced mediation analyses using Aalen additive hazard model^{29}; VanderWeele^{14} presented mediation analyses using Cox proportional hazard models^{30} and accelerated failure time models^{31}; Tchetgen Tchetgen^{32} proposed studying mediation in Cox models by a doubly robust estimator. Except for a resampling-based method presented by Lange et al.,^{33} literature on survival models that incorporate multiple mediators is still very limited.

Because the numerical approach may have its limitation in computation costs, this article aims to present analytic methods, which also provide better mechanistic understanding of the parameters involved in mediation effects. Huang and Cai^{34} have recently developed a new mediation analysis for survival outcome using semiparametric probit models. Instead of focusing on a specific effect scale, for example, hazard difference, hazard ratio, or ratio of survival time, the probit model is able to quantify the effect on the entire survival probability. Due to the conjugacy property of normal distributions, this model can be easily used to study multiple mediators. Moreover, we develop two additional multimediator models that use Aalen additive hazard and Cox proportional hazard models. We illustrate the utility of the three multimediator models in the motivating hepatitis study where we examine path-specific effects on different effect scales, including the difference in transformed survival time, the difference in survival probability across follow-up time, the hazard difference, and the hazard ratio.

## DEFINITIONS AND ASSUMPTIONS

We first define a general outcome *Y* to be a function of time to liver cancer development *T*,

. Let

be the counterfactual outcome that would have been observed had *S* (HCV baseline viral load), *M*_{1} (HBV baseline viral load), and *M*_{2} (HBV follow-up viral load) been set to *s*, *m*_{1}, and *m*_{2}, respectively;

be the counterfactual value of *M*_{2} (HBV follow-up viral load) had *S* and *M*_{1} been set to *s* and *m*_{1}, respectively; and

be the counterfactual value of *M*_{1} (HBV baseline viral load) had *S* been set to *s*. By extension of natural direct and indirect to two-mediator models, we define three path-specific effects on *Y*:

A series of no unmeasured confounding assumptions needs to be made to identify the three path-specific effects.^{26},^{34} With

for the condition that *A* is independent of *B* conditional on *C*, we list six sufficient conditions for identifying the path-specific effect (PSE) in the context of the motivating hepatitis study. To simplify notation, the following assumptions are presented as marginal exchangeability, and they can be generalized to be conditionally exchangeable on covariates ** X**, in which known confounders such age and gender can be adjusted to achieve exchangeability.

1. : there is no confounding for the joint effect of the baseline and follow-up hepatitis B viral load on the time to liver cancer development, conditional on the baseline hepatitis C viral load (S).

2. : there is no confounding for the effect of the baseline hepatitis C viral load (S) on the time to liver cancer incidence.

3. : there is no confounding for the joint effect of , the baseline HBV and HCV viral load on the follow-up HBV viral load (M2).

4. : there is no confounding for the effect of the baseline HCV viral load (S) on the baseline HBV viral load (M1).

5. : there is no baseline HCV viral load (S)-induced factor that can confound the baseline HBV viral load–survival time (M1–T) and the follow-up HBV viral load-survival time (M2–T joint relation, where and are interventions for the baseline HCV viral load with different values than s and each other.

6. : there is no baseline HCV viral load-induced factor that confounds the M1–M2 (baseline vs. follow-up HBV viral load) association.

Besides no unmeasured confounding assumptions, standard assumptions of positivity and consistency are also required (see eAppendix A1; http://links.lww.com/EDE/B185).

## MULTIMEDIATOR MODELS OF SURVIVAL OUTCOME

In this section, we present three methods for mediation analyses with one exposure *S*, two mediators *M*_{1} and *M*_{2}, and a survival outcome *Y*, as shown in Figure 1. The methods can easily be extended to more than two mediators.

### Aalen Additive Hazard Model

We propose two linear regression models for *M*_{1} and *M*_{2}, respectively:

where ** X** is covariates with the first element being 1 for the intercept; the error terms

and

are independent and normally distributed with mean zero and respective variances,

and

. We propose the following additive hazard model for the outcome

:

where

is the hazard of developing liver cancer for subject *i*;

is the baseline hazard; and

,

,

, and

are regression coefficients for the covariates

(** X** without the first element), the baseline HCV viral load

*S*, the baseline HBV viral load

*M*

_{1}and the follow-up HBV viral load

*M*

_{2}, respectively. Based on the six assumptions in the “Definitions and Assumptions” section, one can express the counterfactual hazard as follows:

where

. Derivation of the above expression is in the eAppendix (http://links.lww.com/EDE/B185). The model for the second mediator (3) and the survival model (4) focus on the main effects, and they can incorporate *S*-by-*M*_{1} and *S*-by-*M*_{2} cross-product interaction terms by replacing

,

, and

with

,

, and

, respectively. We then use the definition in (1) to re-expressed path-specific effects on the scale of hazard difference by using the above result:

The results in (5) have intuitive interpretations and can be easily visualized using the DAG in Figure 1. Each arrow in Figure 1 has its corresponding effect parameter in models (2)–(4):

represents the effect of the baseline HCV viral load *S* on the baseline HBV viral load *M*_{1};

and

, respectively, represent the effects of the baseline HCV viral load *S* and the baseline HBV viral load *M*_{1} on the follow-up HBV viral load *M*_{2};

,

, and

, respectively, represent the effects of *S*, *M*_{1}, and *M*_{2} on the hazard of developing liver cancer. The path of

has only one arrow with effect parameter

.

has two arrows:

and

with respective effect parameters being

and

, and the path-specific effect is the product of the two parameters

.

contains two paths:

and

. The path

consists of three arrows:

,

, and

with effect parameters being

,

, and

, respectively; the path

consists of two arrows:

and

with effect parameters being

and

, respectively. The effect of

is the sum of

and

, which are the products of effect parameters along the two paths.

Total effect, defined as

, can be expressed as the sum of the three path-specific effects:

. The additive hazard model (4) can incorporate time- dependent effects by generalizing

,

,

, and

to

,

,

, and

, respectively.

### Cox Proportional Hazard Model

Unlike additive hazard models, which assume hazard linearly determined by the predictors, Cox models assume that the hazard is determined by the predictors exponentially (or linearly on the log hazard scale):

where

,

,

, and

are regression coefficients for

, *S*, *M*_{1}, and *M*_{2}, respectively. Provided that the assumptions in the “Definitions and Assumptions” section hold and the outcome event is rare, one can approximate the counterfactual log hazard as follows:

where

. Derivation of the above expression is provided in the eAppendix (http://links.lww.com/EDE/B185; Section A5). Again, it can be extended to incorporate *S*-by-*M*_{1} and *S*-by-*M*_{2} cross-product interaction terms by replacing

,

, and

in (3) and (6) with

,

, and

, respectively. We let

, and re-express the three path-specific effects in (1) as effects on log hazard ratio by using the above expression:

Similar to mediation analyses using Aalen additive hazard models, one can visualize the three path-specific effects

as products of parameters in Figure 1. Total effect is sum of the three path-specific effects:

Variances and confidence intervals for

and

can be calculated using a resampling method that takes random draws from multivariate normal distribution of estimates for

or

^{13},^{28} with detail provided in the eAppendix (http://links.lww.com/EDE/B185; Sections A4 and A5).

### Probit Model

For completeness and comparison with the above two models, we present another multimediator model that has been published.^{34} This model proposes a semiparametric probit model for the survival time:

where the outcome *H*(*T*) is a nonparametric transformation

of the survival time *T*;

is a standard normal random variable, independent of

and

; and

,

,

, and

are regression coefficients for the covariates ** X**, the baseline HCV viral load

*S*, the baseline HBV viral load

*M*

_{1}, and the follow-up HBV viral load

*M*

_{2}, respectively. By the assumptions in the section of “Definitions and Assumptions”, one can show that

The derivation is in the eAppendix (http://links.lww.com/EDE/B185; Section A3). Similar to the additive hazard and proportional hazard models, it can be extended to incorporate *S*-by-*M*_{1} and *S*-by-*M*_{2} cross-product interaction terms by replacing

,

, and

in (3) and (8) with

,

, and

, respectively. By letting

and using the above result, the path-specific effects defined in (1) can be expressed as follows:

Note that the three path-specific effects

are effects on the difference in the transformed survival time. Similar to the other two models, the three path-specific effects

can also be visualized as products of parameters in Figure 1.

Next, we present another set of path-specific effects, the difference in survival probability. Note that the survival probability in the motivating hepatitis study represents the probability that a subject has not yet developed liver cancer. Because

,

, and

are all normally distributed, one can further show that the distribution function for the survival time *T*,

is conjugate with the distributions for *M*_{2} and *M*_{1}. The counterfactual outcome defined as a cumulative distribution function of counterfactual survival time

can be expressed as follows:

where

and

is the baseline cumulative hazard (see eAppendix; http://links.lww.com/EDE/B185 for the derivation). Therefore, by letting

, we obtain another set of path-specific effects

that characterizes the effects on entire survival probability as a function of the follow-up time:

The difference in survival probability is a very general effect because other effect scales such as the average survival time and the hazard are functions of the survival probability.

,

,

’s,

’s,

’s all can be estimated consistently with least square estimators and a nonparametric maximum likelihood estimator. It follows that path-specific effects can be estimated with the closed-form formula in (9) and (10) by plugging in the estimates. Total effect is defined as

and

can be expressed as

and

respectively. Variances and confidence intervals for

and

can be calculated using the nonparametric maximum likelihood estimator and the functional delta method.^{34},^{35}

The proposed method for mediation analyses using semiparametric probit models has been implemented in Matlab code, and those using Aalen additive hazard models and Cox proportional hazard models have been implemented in R code. We term these methods M3S (multimediation models for survival data). Relevant code is available on our website (http://www.stat.sinica.edu.tw/ythuang/M3S-Probit.zip; http://www.stat.sinica.edu.tw/ythuang/M3S-Cox%20and%20M3S-Aalen.zip). R scripts for the Aalen and Cox models are also in the eAppendix (http://links.lww.com/EDE/B185; Section A7).

## DATA APPLICATION

Details of the study design are provided in eAppendix A6 (http://links.lww.com/EDE/B185). The ethical review is not required for the secondary analyses of the existing and deidentified data. We first conducted one-mediator analyses, similar to Huang et al.,^{15} to examine the effect of the HCV viral load (*S*) on the incidence of liver cancer (*Y*) mediated through the follow-up HBV viral load (*M*). Natural direct effect, natural indirect effect, and total effect for the HCV viral load increasing from the minimum to the first quartile of the detectable log-transformed viral level are presented in Table 1. Mediation analyses using probit models estimated that an increase in the HCV viral load decreased the transformed survival time *H*(*T*) using the natural direct effect (−0.32 [95% confidence interval {CI}: −0.57, −0.076]), which is not through HBV, but increased *H*(*T*) by the natural indirect effect (0.19 [95% CI: 0.12, 0.26]); the opposite natural direct and natural indirect effects masked each other in the total effect (−0.13 [95% CI: −0.38, 0.12]). We noted that the nonparametric transformation

made it less intuitive to interpret the effects in absolute units such as months. Results using Aalen additive hazard models and Cox proportional hazard models revealed a similar pattern (Table 1). Although the three models provided similar results under different effect scales, probit models can further characterize effects on the survival probability as a function of the follow-up time. The natural direct and indirect effects, respectively, decreased and increased the survival probability from getting liver cancer monotonically across the follow-up time: for an increase from the minimum to the first quartile of the detectable HCV viral level, the effect independent of the follow-up HBV DNA (natural direct effect) decreased the 15-year probability of not having an hepatocellular carcinoma (HCC) event by 0.3%, and the effect mediated by the follow-up HBV DNA (natural indirect effect) increased the probability by 0.18% (Figure 2A). Again, the opposite effects masked each other in total and marginal effects (Figure 2B) where the marginal effect was estimated from the probit model that simply regressed *H*(*T*) on HCV viral load *S*, adjusting for covariates ** X**. Besides the total effect, the marginal effect was another estimate for the overall effect, which did not decompose the effect into the natural direct and indirect effects; but similar to the total effect, the marginal effect was also conditional on

**.**

*X*According to the epidemiology of hepatitis in Taiwan, most dual infected participants were chronic hepatitis B carriers superinfected by HCV.^{17},^{20} The literature supports the assumption that HBV is suppressed by the superimposed HCV.^{18},^{19},^{21–23} Therefore, the HCV (*S*)-induced baseline HBV viral load (*M*_{1}) may be a common cause of the follow-up HBV viral load (*M*_{2}) and the outcome (*Y*) (Figure 1), which violates the identifiability assumption

in the single-mediator model. To address the issue, we resorted to the two-mediator model investigating path-specific effects^{26} (Table 2). In analyses using Cox models, an increase in the HCV viral load increased the risk of liver cancer through the mechanism of

with a hazard ratio of 1.17 (95% CI: 1.06, 1.29), but decreased the risk through the mechanism of

with a hazard ratio of 0.91 (95% CI: 0.88, 0.94); the HCV viral load did not change the risk through the mechanism only mediated by the follow-up HBV viral load (

: 1.00 [95% CI: 0.99, 1.01]); and the total effect was 1.06 (95% CI: 0.96, 1.18). Probit models characterized path-specific effects on the entire survival function, which suggests that

and

, respectively, decreased and increased the survival probability (from getting liver cancer) monotonically across the follow-up time;

covered the probability difference of zero; and the three effects masked each other in total and marginal effects (Figure 3).

The nonparametric transformation function

depicted in Figure 4 is a monotonically increasing function with slight nonlinearity and was estimated as a step function using the nonparametric maximum likelihood estimator.^{34} The normality assumption was also evaluated for probit models using a bootstrap procedure,^{34} which revealed no departure from the normal distribution.

## DISCUSSION

We present multimediator analyses under semiparametric probit models, additive hazard models, and Cox proportional hazard models and provide closed-form expressions for path-specific effects. The closed-form expression for Cox models can only be obtained under the rare-outcome assumption, which may be a limitation for its application. Mediation analyses under additive hazard models and Cox models characterize effects on specific scales; additive hazard models focus on the hazard difference and Cox models focus on the log hazard ratio. Probit models examine the effect on the transformed survival time difference and the difference in survival probability across the follow-up time. The closed-form expressions under the three models have intuitive interpretations using a causal diagram labeled with effect parameters on each arrow. Effects estimated from the three models convey complementary information. The effect on survival probability by probit models is very general and may be used to estimate any summary outcome such as average survival time or 5-year survival risk. The effects estimated using probit models is in the opposite direction compared with the other two models because probit models characterize the difference in survival time and survival probability, which is qualitatively the opposite of the hazard difference or the hazard ratio. Because of the lack of conjugacy between product of normal and normal distributions, we may not be able to provide closed-form expressions for a mediator-by-mediator interaction effect (e.g.,

) or a nonlinear mediator effect (e.g.,

). One may have to utilize the simulation-based approach to approximate these effects.^{33} In the hepatitis study, there was no evidence for such effects.

Causal multimediator models have several advantages. First, models with multiple mediators are able to characterize more comprehensive mechanistic structure. It is often not plausible to assume that the effect of an exposure on an outcome is mediated only by a single mediator. Second, path-specific effects in multimediator models constitute an approach to address identifiability issue due to the undue influence by the exposure-induced mediator–outcome confounder,^{26} for example, the baseline HBV viral load in the hepatitis study. Finally, being able to incorporate more potential mediators renders the identifiability assumptions (5) and (6) more plausible^{25} because the inclusion of additional mediators makes it less likely to miss alternative pathways from the exposure to the outcome.

The analyses of the hepatitis study were limited to participants with no missing data on the hepatitis B and C viral load and covariates, and assumed that the viral load was accurately measured. Mediation analyses accounting for missing data and measurement error have been developed,^{36},^{37} but how to incorporate in the proposed multimediator analyses for survival outcomes may require additional development.

The finding that

had stronger effects than

may not exclude the possibility that the follow-up HBV DNA is a mediator because

contains the effect mediated by both baseline and follow-up HBV DNA levels. However, our analyses does not support that there existed an effect mediated only through the follow-up HBV. The finding has an important public health implication: Taiwanese patients with dual infection of hepatitis B and C may not require more intensive follow-up for HBV activity than those with HBV single infection if their baseline HBV viral load is low. Analytically, one may propose a mediation model with the baseline and follow-up HBV viral load as a single set of mediators, and interrogates the longitudinal mediation effect through HBV.^{38} Such analyses evaluate an overall mediation effect but would not be able to provide the respective effect measures of different mechanisms/paths.

## REFERENCES

Equation (Uncited)