For each individual i = 1, 2, …, n, T _i = min(D _i , C _i ) denotes the observed survival time, where D _i is the time from beginning to the event and C _i is the censoring time. δ _i = I(D _i ≤ C _i ) is the failure indicator, and I(⋅) is the indicator function. When D _i > C _i , the participant is said to be right-censored, which we consider more in this article. Censoring rate, representing the rate of participants whose information is not available due to loss to follow-up or nonoccurrence of the interested event within the trial duration, is significant to survival analysis (Prinja et al., 2010).

Figure 1 is a direct acyclic graph showing the relationship between exposure, outcome, covariates, and high-dimensional mediators. X is the exposure. M = { M 1 , M 2 , … , M p } T denotes the high-dimensional mediators, and p ≫ n. Y is the survival outcome. Z represents covariates. The additive hazards model with mediators is: λ i t | X i , M i , Z i = λ 0 i t + γ X i t + θ T Z i t + ∑ k = 1 p β k M k i i = 1,2 , … , n , (4) M k i = c k + α k X i t + ϑ T Z i t + e k i   k = 1,2 , … , p . (5)

Eq. 4 is an additive hazards model showing individual’s hazard rates. λ _i is associated with exposure, covariates, and high-dimensional mediators. λ _0i (t) indicates the time-varying intercept. Eq. 5 describes the way how exposure and covariates linearly influence mediators. c _k is the intercept and e _ki is random error.

To obtain a causal inference conclusion from the mediation analysis, we make some assumptions about mediators and confounders. Here, T(x, M ₁, M ₂, …, M _p ) denotes that the survival time depends on X and M _k (k = 1, 2, …, p). M _k (x*) represents the mediators with different exposure values. The consistency assumption matters to the proposed procedure requiring to hold the outcome once the exposure and mediators were set (VanderWeele and Vansteelandt, 2009; Rehkopf et al., 2016). Based on Luo et al. (2020) and Huang and Yang (2017), the assumptions for the proposed procedure are as follows:

1) X ⊥ T(x, m ₁, m ₂, …, m _p )|Z; there is no unmeasured confounding effect between X and T conditional on Z.

Referring to the counting process notation, N _i (t) = I(T _i ≤ t, δ _i = 1) represents the observed failure counting process, where δ _i = I(D _i ≤ C _i ). Y _i (t) = I(T _i ≥ t) is the at-risk indicator. And M i t = N i t − ∫ 0 t Y i s λ 0 i s + γ X i s + θ T Z i s + β T M i s d s is the additive martingale process. Let P = (γ, θ, β) and Q _i = (X _i , Z _i , M _i ). Then the martingale could be simplified as M i ( t ) = N i ( t ) − ∫ 0 t Y i ( s ) { λ 0 i ( s ) + P T Q } d s .

According to Lin and Ying (1994), the pseudo-likelihood score function of the proposed model is: U P = ∑ i = 1 n ∫ 0 ∞ Q i t − Q ̄ t d N i t − Y i t P T Q i t d t , where Q ̄ ( t ) = ∑ j = 1 n Y j ( t ) Q j ( t ) / ∑ j = 1 n Y j ( t ) . Referring to Lin and Lv (2013), we can write the score function into U P = b − V P , where b = 1 n ∑ i = 1 n ∫ 0 ∞ Q i t − Q ̄ t d N i t , V = 1 n ∑ i = 1 n ∫ 0 ∞ Y i t Q i t − Q ̄ t ⊗ 2 d t , and a ^⊗2 = aa ^T . Then the least-squares type loss function of the proposed model is: L P = 1 2 P T V P − b T P . (6)

However, the maximum likelihood estimation is not feasible when p ≫ n. To identify the true-positive mediators, we consider the “two-step” method for dimension reduction. First, we apply SIS to reduce dimension from an ultrahigh level to a moderate one (Fan and Lv, 2008). Then we perform the regularization method with SCAD penalty for the SIS-selected subset. The Sobel test is applied to identify true mediators in SCAD-selected subset. Figure 2 shows the overall workflow of the proposed procedure. We will introduce details below.

Step 1. (Screening) Using SIS to reduce candidate mediators from p dimension to d dimension, we identify a subset S ₁ = {M _k : 1 ≤ k ≤ p}. Here we select d = [2n/ log(n)] mediators instead of [n/ log(n)] recommended by Fan and Lv (2008) to contain more positive mediators in subset S ₁, because mediators are related to exposure and outcome simultaneously.

Step 2. (SCAD-penalized selection) Further selection with SCAD penalty for the subset S 2 = { M k : β ̂ k ≠ 0 } based on M _k ∈ S ₁ is applied by minimizing the following objective function with penalty: Q β = L β + ∑ j = 1 p p λ β j , where L( β ) has been shown in Eq. 6, and p λ ′ | β | = λ I | β | ≤ λ + a λ − | β | + a − 1 I | β | > λ a > 2 λ > 0 .

Here we choose the regularization parameters by 5-fold cross-validation. Gorst-Rasmussen and Scheike (2012) implemented the SCAD penalized method for additive hazards model in R package ahaz.

Step 3. (Effect decomposition and IE test) Referring to the single mediator (Lange and Hansen, 2011) and two mediators based on the additive hazards model (Huang and Yang, 2017), we use the counterfactual hazard difference to measure the effect difference when exposure changes from x to x*. Counterfactual hazard difference, also named total effect (TE), includes two parts: direct effect (DE) and IE. DE represents the exposure directly caused effect. And IE expresses the effect caused by exposure through mediators indirectly.

Defining M k x and M k x * as the mediator value with different exposure value x and x* separately, we have the following decomposition of TE (for more details see Supplementary Appendix): TE = λ T x * , M 1 x * , … , M p x * ; t | Z − λ T x , M 1 x , … , M p x ; t | Z = λ T x * , M 1 x * , … , M p x * ; t | Z − λ T x * , M 1 x , … , M p x ; t | Z + λ T x * , M 1 x , … , M p x ; t | Z − λ T x , M 1 x , … , M p x ; t | Z = γ x * − x + α 1 β 1 + ⋯ + α p β p x * − x = DE + IE

Then we apply the Sobel mediation significance test to subset S ₂ to pick out true-positive mediators from candidates by significant IE. According to Sobel (1987), we have the null hypothesis H ₀: α _k β _k = 0 and following p value calculating formula: P r a w , k = 2 1 − ϕ | α k ̂ β k ̂ | σ ̂ α k β k , (7) where σ ̂ α k β k = α ̂ k 2 σ ̂ β k 2 + β ̂ k 2 σ ̂ α k 2 is the estimated standard error, and α ̂ k is the estimator of α _k , β ̂ k is the estimator of β _k , σ ̂ α k 2 is the estimated variance of α _k , and σ ̂ β k 2 is the estimated variance of β _k .

This section demonstrates the simulation results of the proposed procedure with high-dimensional mediator’s selection and IE estimation in a series of simulation studies.

Lung cancer is still the most fatal cancer worldwide with many pathogenic factors such as tobacco smoking and air pollution; 80% to 85% of lung cancers were caused by smoking (Wild et al., 2020). Nicotine in tobacco may result in genetic mutations. To find out whether smoking leads to lung cancer by affecting the DNA methylation, we applied the proposed procedure to the TCGA lung cancer cohort study involving DNA methylation data (907 samples measured by Illumina Infinium HumanMethylation450 platform), phenotype data (1,299 samples), and survival data (1,145 samples) for lung squamous cell carcinoma and lung adenocarcinoma. DNA methylation values recorded via BeadStudio software were continuous from 0 to 1 representing the intensity ratio. Thus, a higher value represents a higher degree of methylation, and so does the lower one.

After sample matching and data cleaning among the above data sets, we obtained 833 patients; 41.2% (343) were female, and 68.4% (570) were smokers. The patients’ ages ranged from 33 to 90 years with a median of 67 years. The overall survival time represented the days from first diagnosed to death or the last follow-up date. The median survival time was 652 days (1.79 years).

SIS based on the marginal correlation between tobacco smoking status and DNA methylation was first applied to reduce DNA methylation sites from 365,306 to 2n/log(n) (=248). Then we applied the SCAD penalty for a further dimension reduction and get a 25 sites subset. We applied the Sobel test and BH method to that subset for IE hypothesis testing. cg19757631, cg08636115, cg05147638, cg24720672, and cg08530838 are significant DNA methylation sites with adjusted p value < 0.05. We are interested in mediating DNA methylation markers, which increase lung cancer patients’ survival hazards. Therefore, we focus on the CpG sites with positive IEs ( α k ̂ β k ̂ > 0): cg19757631, cg08636115, cg05147638, and cg24720672.

Table 3 shows mediated CpG sites with positive IE. The estimated IE was represented by α ̂ β ̂ . The TE (effect between exposure and outcome with covariates) of tobacco smoking on lung cancer patients’ survival equaled 0.0137 (95% CI = −0.0252–0.0526), and its DE (effect between exposure and outcome adjusting for mediators and covariates) equals 0.0171 (95% CI = −0.0244–0.0585). The IEs of four significant mediated CpG sites cg19757631, cg08636115, cg05147638, and cg24720672 are equal 0.0296 (95% CI = 0.0129–0.0464), 0.0263 (95% CI = 0.0093–0.0433), 0.0185 (95% CI = 0.0047–0.0323), and 0.0269 (95% CI = 0.0100–0.0438), respectively.

Bakulski et al. (2019) studied DNA methylation sites associated with smoking exposure in TCGA lung adenocarcinoma tissue samples and found cg19757631 is significant (FDR-adjusted p value < 0.05). In their study, the estimated methylation change of smokers versus never smokers is −12.28% (adjusted p value = 4.81E-06), which is consistent with ours ( β ̂ = − 0.2806 ). The experiment results from Fei et al. (2019) about gene PRDM16 (cg08636115 located) suggest that PRDM16 is a metastasis suppressor and potential therapeutic target for lung adenocarcinomas, which has the same conclusion as ours ( β ̂ = − 0.3811 ) . Shtutman et al. (2011) explained the operational mechanism of COPZ1 (cg05147638 located) in the tumor cell: the function-based genomic screening identified COPZ1 gene is essential in different tumor cell types instead of normal cells. Gene COPZ1 methylation is harmful. This conclusion approves our results: β ̂ = 0.4918. As for CpG site cg24720672, we find some researches about leukemia—a kind of cancer, and we infer it has the similar mechanism in tumor tissue as lung cancer (Nair, 2016; Zhang et al., 2018; Jiang et al., 2020).

The real data application identifies four significant mediated DNA methylation sites with positive IEs between tobacco smoking and lung cancer patients’ survival. CpG site cg19757631 is a mediator having a known relationship with tobacco smoking (Bakulski et al., 2019). CpG sites cg05147638, cg08636115, and cg24720672 are newly addressed mediators. Besides, we also apply the naive method to the TCGA lung cancer data set, but nothing has been identified.