A cornerstone of epidemiologic research is the understanding of causal pathways from an exposure to an outcome. For example, how much of the observed increase in the risk of breast cancer among postmenopausal women from alcohol consumption or obesity is mediated through an increase in women's estrogen levels? Cohort studies could provide measurements of risk factors and estrogen levels at baseline and data on time to diagnosis of breast cancer or study end (whichever comes first). A second example, which will be used throughout this paper, is the effect of socioeconomic status (SES) on the risk of experiencing long-term absence due to sickness. It is well documented that SES strongly affects the risk of experiencing long-term sickness absence. This effect seems to be at least partly mediated through differences in the physical work environment among SES groups,1 but how much?
A common causal structure underlies many epidemiologic problems. In the case of an outcome that is either binary or normal, a number of techniques have been developed to assess the relative importance of the various paths from the exposure to the outcome.2–6 However, as illustrated by our 2 examples, the outcome of interest will often be time-to-event (or, in other words, survival time). In the SES example, this would be number of weeks from baseline to onset of first long-term sickness absence or study end. As is well known, survival times are non-normal and will typically be right censored.
In this paper, we propose a simple measure of mediation in a survival setting. The measure is based on the counterfactual framework7 and has direct causal interpretation. The method allows the total effect (adjusted for confounders) of changing an exposure, eg, SES, on a time-to-event outcome to be measured as the number of additional events per unit of time. Furthermore, the method can decompose this number of additional events into a part attributable to the direct pathway and a part mediated through a given mediator. In the SES example, the method enables us to compute how many long-term sickness cases could be prevented per week if the SES groups were changed from the lowest to the highest. This number of prevented cases can then be decomposed into a part mediated by physical work environment and a direct part.
To assess the magnitude of the pathway from SES through physical work environment to time to onset of long-term sickness absence, the traditional approach estimates hazard ratios of SES from Cox models both with and without adjusting for physical work environment (the potential mediator). A change in hazard ratios is taken as evidence of mediation through physical work environment. However, as previously pointed out in the literature,8,9 such an analysis of mediation has severe shortcomings. Most importantly, the observed changes in hazard ratios cannot be given a causal interpretation. In addition, the important assumption of proportional hazards can never be satisfied for both models with and without the mediator. In other words, it is not mathematically consistent to use a Cox model both with and without a potential mediator (mathematically, this is due to the fact that the class of proportional hazard models is not closed under marginalization).
The paper proceeds as follows. We introduce the model framework and describe the counterfactual variables and the required assumptions. Next, we present our proposed measures of mediation and illustrate the technique on the problem of mediation in the relationship between SES and long-term sickness absence. The subsequent section discusses generalizations and limitations of our proposed method. All technical details are deferred to eAppendix 1 (https://links.lww.com/EDE/A476), eAppendix 2 (https://links.lww.com/EDE/A476) contains a detailed guide to implementation in the software package R, and eAppendix 3 (https://links.lww.com/EDE/A476) presents a small simulation study.
THE MODEL FRAMEWORK
Assume that measurements have been made of the time to an event (T) or censoring, whichever comes first, an exposure (X), a potential mediator (M), and other baseline covariates (Z). In the SES example, T is time to onset of long-term sickness absence; X is SES group; M is physical work environment; and Z is sex, age, etc. Recall that the rate at time t measures the probability of experiencing an event within the next unit of time, given that a person has not experienced an event before time t. The usual Cox model applied to the SES example yields an estimate of how many times greater the rate is for a given SES group relative to a reference group, ie, the hazard ratio. However, because this ratio cannot be related to an absolute number of events and the mathematical structure of the Cox model fits poorly with mediation analysis, we suggest modeling the rate by the Aalen additive hazard model10 instead. Applied to the SES example, this model yields an estimate of the absolute change in the rate when comparing a given SES group to a reference group. This estimate (multiplied by, say, 10,000) can be directly interpreted as the number of additional cases of long-term sickness absence per week per 10,000 persons at risk, when compared with the reference group. The Aalen additive hazard model is a flexible semiparametric model for survival outcomes and is at least as flexible as the more widely used Cox model. In the Aalen model, the rate as a function of exposure (x), mediator (m), and other baseline covariates (z) can be written as
;)
where λj(t) are potentially time-dependent coefficient functions. We have allowed for all coefficients to be time-dependent (eg, the effect of SES might change over time), but in many applications non–time-dependent coefficients will suffice. Standard software for estimating additive hazard models (eg, the “timereg” package in the software package R) have built-in tools for determining which coefficients need to be time-dependent.
In addition, assume that the mediator is a normal variable that can be modeled by a simple linear regression. Thus, given exposure (x) and other covariates (z), the mediator can be written as
where e is a mean zero normally distributed error with variance σ2. Given a set of independent observations, standard techniques (eg, the “timereg” package) can be used to estimate the parameters α0, α1, α2, σ2, and the collection of functions λ0(t), ..., λ3(t).
COUNTERFACTUALS AND ASSUMPTIONS
To model the causal effect of the exposure, we also define the following counterfactual variables, which are variables describing what would have happened if, perhaps contrary to fact, exposure and mediator were set to specific values. For an introduction to counterfactuals see the book by Pearl.11
indicates the time to event when the exposure is set to x and the mediator is set to m.
indicates the value of the mediator when the exposure is set to x.
In the SES example,
describes the time to onset of long-term sickness absence person i would experience if SES group was set to x and physical work environment set to m. To facilitate the discussion of mediation, note that the aforementioned counterfactuals can be combined—as for instance
which denotes the event time when the exposure is set to x, but the mediator is set to the value it would have had if the exposure had been set to x*. As it is impossible to measure individual causal effect, we will drop the subscript i on all counterfactuals and discuss only population-wide causal effects (also called average causal effects).
Following the tradition in mediation analysis,4 we will compare
with
to obtain a measure of the natural direct effect of changing the exposure from x to x*. Likewise, we will compare
with
to obtain a measure of the natural indirect effect. In this context, the word “natural” refers to the fact that we let the mediator take the value it would take naturally when the exposure is set to x. This is in contrast to controlled effects, where the mediator is kept fixed at a controlled level m. As controlled effects do not allow for an obvious definition of indirect effects, we will focus on natural effects in this paper. Natural and controlled effects have been discussed extensively.7,3,12
For future use, let γ(t; x, m) denote the counterfactual rate for the event when the exposure is set to x and the mediator to m; that is, the rate for the counterfactual variable Tx,m.
ASSUMPTIONS
Whenever one wishes to draw causal conclusion from an analysis, a number of assumptions about absence of unmeasured confounders (sometimes referred to as exchangeability assumptions) must be satisfied. In particular, we must assume that there are no unmeasured confounders for the exposure-outcome (A.1), the mediator-outcome (A.2), and the exposure-mediator relations (A.3).
Furthermore, to identify natural effects we adopt a modified version of Pearl7 identifiability condition: (A.4)
. On an intuitive level, this assumption ensures that the effect of the exposure has its effect through 2 distinct and nonintertwined causal pathways. The assumption would, for instance, be violated if there is a variable that is affected by the exposure and itself affects both the mediator and the outcome. Indeed, it has been shown13 that if such an exposure-dependent confounder (for the mediator-outcome relation) exists, then direct and indirect effects will not be nonparametrically identified. This applies to any mediation analysis, not just in a survival setting. The literature contains alternative (and less restrictive) versions of this identification assumption.3,14,15 Common to these, however, is the fact that they only work in the simpler setup of a nonsurvival outcome.
Finally, we must impose the consistency assumptions,16 denoted (A.5). This ensures that we do not alter the outcome if we set the exposure and mediator to the values they would naturally take.
THE PROPOSED MEASURE OF MEDIATION
Unlike traditional mediation analysis, where the outcome is measured at a single time point, a causal measure of mediation in survival analysis must take into account (1) how much of the effect of the exposure is mediated through the mediator and (2) how these proportions change over time. We therefore suggest using the counterfactual rate difference as the effect measure of the exposure changing from x to x*. In the SES example, the counterfactual rate difference measures the added number of long-term sickness absence cases per week caused by changing SES group from x to x*. In many respects, the counterfactual rate difference can be interpreted as a counterfactual risk difference.17 As will be demonstrated later in the text, this measure can naturally be separated into 2 parts; one measuring the natural direct effect and the other measuring the natural indirect effect. The proof of the following theorem can be found in eAppendix 1 (https://links.lww.com/EDE/A476).
Theorem 1.
Under assumptions (A.1) to (A.5), it holds that the total causal effect of changing the exposure from x* to x, measured on the rate difference scale at time t, can be expressed as
with TE, DE, and IE denoting total effect, natural direct effect, and natural indirect effect, respectively.
To interpret the counterfactual rate differences (ie, TE(t), DE(t), and IE(t) in Theorem 1), imagine 3 versions of the world at baseline. In the first version, X is set to the value x* for the whole population and the mediator to the value it naturally takes when X is set to x*. In the second version of the world, X is kept at x*, but the mediator is set to the value it naturally takes when X is set to x. In the third version of the world, X is set to the value x and the mediator to the value it naturally takes when X is set to x. In all other regards, these 3 versions of the world are identical at the time when X is set. At time t, we take 3 large samples (say 10,000 persons) from the persons still alive in each of the 3 versions of the world. The total effect from Theorem 1 (multiplied by 10,000) is exactly the difference in the number of deaths in the sample from version 1 and the sample from version 3 per unit of time. Likewise, the direct effect is the difference in the number of deaths per unit of time between the samples from version 2 and 3. Finally, the indirect effect is the difference in the number of deaths per unit of time between the samples from version 1 and 2.
In summary, the indirect effect is the number of deaths that can be attributed to mediation through the mediator, whereas the direct effect is the number of deaths that can be attributed to a direct path (or to other mediators not included in the analysis). The total effect is the number of deaths caused by changing the exposure—ie, equals the sum of the direct and indirect effects. In the next section, the use of the theorem is illustrated by an application to the SES example.
If neither the exposure nor the mediator has time-dependent effects in the Aalen model (ie, that λ1(t) and λ3(t) are both constant), Theorem 1 simplifies to
Hence the direct as well as the indirect effects can be expressed by a simple number instead of a function of time (t).
COMPUTING CONFIDENCE INTERVALS
This section explains how to compute 95% confidence intervals (CIs) for total, direct, and indirect effects when the effects are non–time-dependent, as in equation (3). The slightly more involved case of time-dependent effects is addressed in eAppendix 1 (https://links.lww.com/EDE/A476).
Replacing λ1, λ3, and α1 in equation (3) by their estimates 1, 3, and 1, which are available from standard statistical software, yields our estimate of the natural direct and indirect effects. Under mild regularity conditions (eg, condition 5.1 of Martinussen and Scheike17), it holds that the 3 estimators are asymptotically normally distributed and that (1, 3) is uncorrelated with 1. The covariance matrices for both 1 and 3, and 1 are available as output from any software package capable of performing the estimation.
From equation (3), it is evident that the direct effect is asymptotically normal, the indirect effect is asymptotically distributed as the product of 2 uncorrelated normal random variables, and the total effect is the sum of the direct and indirect effects. Confidence intervals and tests involving either the indirect effect or the total effect can be computed using the delta rule, or better yet by simulation. In eAppendix 2 (https://links.lww.com/EDE/A476) we provide detailed guidelines on implementing the techniques in R. Additional R-code is available from the authors upon request.
ILLUSTRATION
To illustrate the new approach we reexamine the relation between socioeconomic status (SES, grouped in 5 categories with “I” being the highest), work environment, and onset of long term sickness absence that is analyzed in the study by Christensen et al.1 This example has been referred to as the SES example earlier in the text. In the paper, the authors analyzed Danish data obtained from a random sample of 11,437 people collected as part of the Danish Work Environment Cohort Study. Long-term sickness absence was defined as 8 weeks of consecutive sickness absence, and follow-up covered 18 months; for further details of the study see Christensen et al.1 To assess the proportion of the effect of SES on long-term sickness absence that is mediated through work environment, Christensen et al1 estimated a number of Cox models, each including additional covariates. If inclusion of a potential mediator in the Cox model reduced the hazard ratios associated with SES, this was interpreted as mediation through that factor. We reanalyze the problem using the approach developed in this paper with a newer data set collected in December 2005. The respondents were followed for 18 months from 1 January 2006 to 30 June 2007. Note that, in this new data set, long-term sickness absence is defined as 3 weeks of consecutive sickness absence. Descriptive statistics are presented in Table 1. All analyses were conducted in R (version 2.8.1).
TABLE 1: Descriptive Statistics by Sex and Socioeconomic Status
For ease of presentation, we will focus only on the physical work environment as mediator. The physical work environment is quantified on a score between 1 and 100 (with 100 corresponding to a physically very demanding work environment), obtained as the average of 5 questions. We will use the log transform of this score, log(phys), which by inspection appears normally distributed.
Table 2 presents the results of a Cox-based mediation analysis similar to the one conducted by Christensen et al.1Table 2 shows that the updated data set yields qualitatively the same results. For example, men in the lowest SES group (V) have a 3.25 times greater hazard of entering a long-term sickness absence period than men in SES group I, when adjusted for age and family status. When log(phys) is included, this hazard ratio drops to 2.04, indicating mediation through the physical work environment.
TABLE 2: Hazard Ratios (HRs) and 95% Confidence Intervals for Onset of Long-term Sickness Absence During 18 Months of Follow-up Using Cox Models, by SES
The first step of our approach is to run a regression of log(phys) on SES adjusted for age and family status. Diagnostic tests indicate that the model is well specified; Table 3 presents the estimates. The next step is to fit the Aalen additive model to the onset of long-term sickness absence with age, family status, and log(phys) as covariates. Standard techniques17 find no indication of time-dependent effects. Successively including interaction terms shows that none is significant after accounting for the number of tests conducted. Diagnostic plots (not reported) indicate that the model is well specified. Non–time-dependent parameter estimates are presented in Table 4. A detailed step-by-step description, including required computer code, of the estimation and subsequent calculation of the direct and indirect effects can be found in eAppendix 2 (https://links.lww.com/EDE/A476).
TABLE 3: Parameter Estimates and Standard Errors (SEs) for the Regression (
Equation 2) of Log(phys) on SES Adjusting for Age and Family Status
Parameter Estimates and Standard Errors From the Aalen Additive Hazard Model (
Equation 1) Adjusting for Log(phys), SES, Age, and Family Status
Considering again the men in SES group V, it can be observed from Table 3 that these on average have a log(phys) score of 1.83 units higher than the men in SES group I when adjusted for age and family status.
As the coefficients of SES and log(phys) in the Aalen model are non–time-dependent, Theorem 1 and equation (3) can be used to separate the total causal effect of moving from one SES group to another into a direct-effect part and a part mediated through log(phys). Confidence intervals for the direct effects are available directly from the Aalen additive model, whereas confidence intervals for the indirect and total effects are constructed by combining the covariance matrices for the parameter estimates of the regression and the Aalen model using simulation. Further details of this procedure are mentioned in eAppendix 2 (https://links.lww.com/EDE/A476).
The direct effects can be found in Table 4 and the total and indirect effects in Table 5. As expected, high SES causes a lower incidence of long-term sickness absence. For instance, among men, a change from SES group I to V would increase the number of long-term sickness absence cases by 10.3 (95% CI = 6.5–14.1) persons per week per 10,000 men at risk. Of these additional cases 4.3 (2.4–6.2) or 43% (22%–74%) can be attributed to the pathway through physical work environment. In other words, if an intervention could improve the physical work environment of the SES V group to the level of the SES I group without affecting other aspects of social deprivation, 43% of the socioeconomic effect on long-term sickness absence could be eliminated.
TABLE 5: Indirect Effects (IE) and Total Effects (TE) According to Theorem 1 and the Ratio of Indirect to Total Effect
Likewise, Table 4 shows that men in SES group V have a rate that is 6.01 × 10−4 units higher than for the men in SES group I, adjusted for age, family status, and log(phys). Thus, the direct effect of changing SES group from I to V is 6.01 additional long-term sickness cases per week per 10,000 men at risk.
Qualitatively, this is in accordance with the results of Christensen et al1 and the results in Table 2, but with the befit of allowing a causal interpretation of the results and with confidence intervals for the ratio of indirect to total effect. For women, the confidence intervals for this fraction corresponding to SES groups III and IV include values above 1. The explanation for this apparent paradox is that the direct and indirect effects could have different signs, such that the total effect is smaller than the indirect effect.
DISCUSSION
To the best of our knowledge, the only method in the literature to quantify mediation in a survival context is dynamic path analysis.18–20 However, this procedure does not provide a causal interpretation and cannot be implemented through standard software. In fairness, it should be said that the dynamic path-analysis framework allows both the exposure and the mediator to be time dependent, which we do not consider in this paper. Recent work21,22 has discussed how to extend formal mediation analysis to nonlinear models (eg, logistic models). Our work can be viewed as a continuation of these principles to a survival context.
In this paper, we focus on the somewhat restrictive case of a conditional normal mediator and no interaction involving either the exposure or the mediator. The simplicity of the expressions for both the natural direct and indirect measures hinges on these assumptions. If they were relaxed, the counterfactual framework could indeed be adapted in an obvious way, but deriving formulas for direct and indirect effects would be challenging as the exact probability distribution of the mediator must be employed to obtain an expression for the causal rate. This generalization is therefore left for future research.
By adapting Aalen's additive hazard model, we obtain a simple and intuitively understandable measure of mediation, namely the number of additional cases per unit of time. In addition, the separation of effects can be performed based on output from standard statistical software. It should be noted that similar results are not, as far as we can tell, attainable with Cox models. However, it is in fact possible to infer a controlled direct effect (on the hazard-ratio scale) from a Cox model by including both exposure and mediator, but neither natural nor indirect effects can be addressed.
We consider a single event type (which covers most epidemiologic applications) for ease of presentation. However, the results also hold with multiple event types. Indeed the mathematical proofs are given in the more general form of multiple event types.
In the analysis of the SES example, the choice of covariates was similar to those in the study of Christensen et al.1 A potential confounder for the relationship between log(phys) and long-term sickness absence could be attitude toward health, which itself can be affected by SES, but no such measurement was available. Thus, we cannot rule out that some of the observed mediated effect of physical work environment is due to difference in attitude toward health.
CONCLUSION
This paper has presented a novel way of quantifying mediation in a survival context. The method provides a simple and straightforwardly interpretable measure of both the natural direct and indirect effects along with their confidence intervals. Effects are calculated on the additive hazard scale and can therefore be directly translated to expected number of additional cases per unit of time.
The method is illustrated by analyzing the relationship between SES, physical work environment, and long-term sickness absence previously examined in the study by Christensen et al.1 As in the original analysis, we found that a substantial part of the effect of SES on long-term sickness absence is mediated through the physical work environment. In addition, our analysis provides confidence intervals for the mediated proportion, and the effect measure can be directly interpreted as the number of additional cases of long-term sickness absence due to differences in the physical work environment and the number due to direct effects. eAppendix 2 (https://links.lww.com/EDE/A476) of this paper provides a detailed step-by-step description of the analysis (including computer code). A major issue for future research is to extend the approach to more general types of mediators and to incorporate interactions.
ACKNOWLEDGMENTS
We gratefully acknowledge valuable discussions with Naja H. Rod.
REFERENCES
1.Christensen KB, Labriola M, Lund T, Kivimäki M. Explaining the social gradient in long-term sickness absence: a prospective study of Danish employees.
J Epidemiol Community Health. 2008;62:181–183.
2.Rubin DB. Direct and indirect causal effects via potential outcomes.
Scand Stat Theory Appl. 2004;31:161–170.
3.Petersen ML, Sinisi SE, Laan MJ. Estimation of direct causal effects.
Epidemiology. 2006;17:276–284.
4.Hafeman DM, Schwartz S. Opening the black box: a motivation for the assessment of mediation.
Int J Epidemiol. 2009;38:838–845.
5.Vansteelandt S. Estimating direct effects in cohort and case-control studies.
Epidemiology. 2009;20:851–860.
6.VanderWeele TJ. Marginal structural models for the estimation of direct and indirect effects.
Epidemiology. 2009;20:18–26.
7.Pearl J. Direct and indirect effects. In: Proceedings of the Seventeenth Conference on Uncertainty and Artificial Intelligence. San Francisco, CA: Morgan Kaufmann; 2001;411–420.
8.Cole SR, Hernán MA. Fallibility in estimating direct effects.
Int J Epidemiol. 2002;31:163–165.
9.Kaufman JS, Maclehose RF, Kaufman S. A further critique of the analytic strategy of adjusting for covariates to identify biologic mediation.
Epidemiol Perspect Innov. 2004;1:4.
10.Aalen O. A model for non-parametric regression analysis of counting processes.In: Klonecki W, Kozek A, Rosinski J, eds.
Lecture Notes in Statistics, Mathematical Statistics and Probability Theory. New York: Springer-Verlag; 1980:1–25.
11.Pearl J.
Causality: Models, Reasoning, and Inference. New York: Cambridge University Press; 2000.
12.Goetgeluk S, Vansteelandt S, Goetghebeur E. Estimation of controlled direct effects.
J R Stat Soc Ser B Stat Methodol. 2009;70:1049–1066.
13.Avin C, Shpitser I, Pearl J. Identifiability of path-specific effects. In: Proceedings of the International Joint Conferences on Artificial Intelligence; 2005:357:363.
14.Robins JM, Greenland S. Identifiability and exchangeability for direct and indirect effects.
Epidemiology. 1992;3:143–155.
15.Hafeman DM, VanderWeele TJ. Alternative assumptions for the identification of direct and indirect effects.
Epidemiology. 2011;22:
16.VanderWeele TJ, Vansteelandt S. Conceptual issues concerning mediation, interventions and composition.
Stat Interface. 2009;2:457–468.
17.Martinussen T, Scheike TH.
Dynamic Regression Models for Survival Data. New York: Springer; 2006.
18.Fosen J, Ferkingstad E, Borgan Ø, Aalen O. Dynamic path analysis—a new approach to analyzing time dependent covariates.
Lifetime Data Anal. 2006;12:143–167.
19.Martinussen T. Dynamic path analysis for event time data: large sample properties and inference.
Lifetime Data Anal. 2010;6:85–101.
20.Commenges D, Gégout-Petit A. A general dynamical statistical model with causal interpretation.
J R Stat Soc Ser B Stat Methodol. 2009;71:719–736.
21.Pearl J.
The Mediation Formula: A Guide to the Assessment of Causal Pathways in Non-Linear Models [technical report]. Los Angeles, CA: University of California; 2010.
22.VanderWeele TJ, Vansteelandt S. Odds ratios for mediation analysis for a dichotomous outcome.
Am J Epidemiol. 2010;172:1339–1348.
