key: cord-0638335-holm86rd
authors: Benkeser, David; D'iaz, Iv'an; Ran, Jialu
title: Inference for natural mediation effects under case-cohort sampling with applications in identifying COVID-19 vaccine correlates of protection
date: 2021-03-03
journal: nan
DOI: nan
sha: 5bb19b10abc148c40406616163ce9d8c3c4e2a40
doc_id: 638335
cord_uid: holm86rd

Combating the SARS-CoV2 pandemic will require the fast development of effective preventive vaccines. Regulatory agencies may open accelerated approval pathways for vaccines if an immunological marker can be established as a mediator of a vaccine's protection. A rich source of information for identifying such correlates are large-scale efficacy trials of COVID-19 vaccines, where immune responses are measured subject to a case-cohort sampling design. We propose two approaches to estimation of mediation parameters in the context of case-cohort sampling designs. We establish the theoretical large-sample efficiency of our proposed estimators and evaluate them in a realistic simulation to understand whether they can be employed in the analysis of COVID-19 vaccine efficacy trials.

The best hope for stemming the SARS-CoV2 pandemic is to develop a safe and effective preventive vaccine and to distribute the vaccine widely (Corey et al., 2020) . Several vaccines have already been demonstrated to exhibit high efficacy and have been granted emergency use authorizations around the world (Baden et al., 2021; Polack et al., 2020) . However, supply issues for COVID-19 vaccines remain a pressing concern. Thus, it is of critical importance to continue to bring new effective vaccines to the market.

An important means of accelerating the vaccine development process is to establish immune correlates of vaccine efficacy. Immune correlates are immunogenicity assays that are predictive of the vaccine's effect on infection with SARS-CoV2 and/or COVID-19 disease (Qin et al., 2007) . These immune responses are typically measured shortly after receipt of the final dose of a vaccine. For SARS-CoV2, scientists are keenly interested in antibody responses to the vaccine insert SARS-CoV-2 proteins, including binding, pseudovirus neutralizing, and live virus neutralizing antibodies. Discovery and validation of a strong immune correlate would have an immense impact on SARS-CoV2 vaccine development, as it could provide a valid surrogate endpoint and thereby open an expedited regulatory pathway for licensure of vaccine products. Currently, in order to adequately power a Phase III licensure trial in the U.S., studies are enrolling between 30,000 and 60,000 individuals and following for up to two years. If a licensure pathway were opened based on an immune correlate, licensure trials could be designed to yield results in weeks rather than years.

Because of the importance of immune correlates to the success of a vaccine development program, sub-studies of ongoing Phase III randomized trials are designed specifically to examine immune correlates in a harmonized fashion (Gilbert et al., 2021) . There is a wide body of literature available on statistical approaches for establishing immune correlates. These methods broadly fall into two classes: correlates of risk and correlates of protection (also known as, respectively, non-mechanistic and mechanistic correlates of protection). The former approaches seek to evaluate relative and absolute risks of infection and disease across levels of an immune response. These methods can be framed as a standard supervised learning problem, where one is predicting infection/disease status based on one or several sets of participant information. Mechanistic correlates of protection, on the other hand, seek to establish a causal relationship between vaccination, immune response and outcome. There is a rich literature available on statistical approaches based on principal stratification (Gilbert and Hudgens, 2008; Joffe and Greene, 2009; Li et al., 2010; Wolfson and Gilbert, 2010) . An alternative approach is to evaluate immune responses as mediators of the vaccine's effect (Cowling et al., 2019) .

Mediation is a very deep literature with contributions coming from across many diverse fields over several decades. VanderWeele (2016) and Nguyen et al. (2020) provide excellent overviews of current research in the area. In this work, we focus on natural mediation parameters (Robins and Greenland, 1992; Pearl, 2001) . These parameters decompose a vaccine's effect into direct pathways (i.e., pathways not involving the measured immune responses) and indirect pathways (i.e., pathways through the measured immune responses). Natural indirect effects involve a comparison of two counterfactual risks: (i) the risk of infection or disease under a hypothetical intervention wherein individuals receive the vaccine, but rather than allowing the vaccine to determine the level of immune response, we fix their immune responses to the values they would naturally assume under placebo; (ii) the risk of infection or disease when individuals receive the vaccine and their immune responses are allowed to naturally respond to the vaccine. Comparison of these risks yields one causal effect that quantifies the impact of the immune response on risk of infection and disease.

Identification and estimation of natural mediation effects has been extensively studied (among others, Petersen et al. (2006) ; van der Laan and Petersen (2008) ; Imai et al. (2010) ; Tchetgen Tchetgen and Shpitser (2012); Zheng and van der Laan (2012)). However, none of these approaches is well-suited for use in the sampling frame of COVID-19 vaccine efficacy trials, where immune responses are measured using case-cohort sampling (Breslow et al., 2013) . In this approach, all Phase III trial participants have blood drawn shortly after receipt of the final vaccination. However, immune responses are only measured in a stratified random sub-sample (i.e., cohort) of individuals, with inclusion probabilities dependent on participants' age, race/ethnicity, and SARS-CoV2 baseline serostatus. In order to adequately power correlates analysis, data are augmented to include immune responses from all individuals who eventually become infected with SARS-CoV2 (i.e., cases). Thus, the observed data represent a biased sub-sample of the trial population.

In this work, we describe methods for estimation and inference on natural mediation effects in this context. Our work builds on past work on multiply robust estimation (Tchetgen Tchetgen and Shpitser, 2012; Zheng and van der Laan, 2012) . We combine these approaches with theory of two-phase sampling (Rose and van der Laan, 2011) to establish the nonparametric efficiency bound for regular, asymptotically Normal estimators in this problem. We establish the efficient influence function of these parameters in a nonparamet-ric model, which naturally suggests one-step estimators of the effects of interest. Using these results we propose two multiply robust estimators. One is derived using a direct application of the strategy suggested by Rose and van der Laan (2011) . The second is novel in that it relies on a new representation of the efficient influence function. For both estimators, we establish regularity conditions under which they converge weakly to Normal random variables. Furthermore, we study and compare their performance in finite samples using simulations targeted towards the design of the randomized COVID-19 vaccine efficacy trials.

Suppose that we have access to a randomized trial that generates independent observations of the vector X = (W, A, S, C, CY ) ∼ P X , where W is a vector of baseline covariates, A is the randomized vaccine assignment (A = 0 denotes randomization to placebo, A = 1 to vaccine), S is a vector of immune responses, C is an indicator of full follow-up for a participant (i.e., that the participant is not lost-to-followup), and Y is the study endpoint (e.g., infection with SARS-CoV2 or clinical COVID-19 disease). Here, we represent the data in such a way that the outcome CY is set to zero is C = 0, though this is arbitrary and does not affect our subsequent developments. Note that this trial is hypothetical, because in reality we will not have observed the immune responses S for all trial participants. We assume a model for P X that makes certain positivity assumptions on the conditional distributions of A, S and C. Specifically, for a 0 = 0, 1, we define g A|W (a 0 | W ) := P X (A = a 0 | W ) and assume a positivity condition, P X {δ A < g A|W (a 0 | W ) < 1 − δ A } = 1 for some δ A > 0. With respect to censoring, we assume that C does not depend on S, which seems reasonable in our context since individuals are blinded to the value of their immune response measurements during the trial. The methods below could easily be adapted to handle the case where censoring depends on S. Moreover, we define g C (1 | a 0 , w) := P X (C = 1 | A = a 0 , W = w) and assume a positivity condition P X {g C (1 | a 0 , W ) > δ C } = 1 for some δ C > 0. Finally, we make a positivity assumption on the conditional immune response distribution. For a 0 = 0, 1, we denote by q S|a0,W (· | W ) the conditional density of S given A = a 0 , W . We assume for a 1 = 0, 1 and a 2 = 1 − a 1 , P X {q S|a2,W (S | W ) / q S|a1,W (S | W )] < ∞ = 1, which is an assumption of common support of the W -conditional immune response distributions between vaccinated and placebo-recipients. It is possible that this assumption could be violated in the context of COVID-19 vaccine trials for certain vaccines and particular immune responses. In particular, for a highly immunogenic vaccine, it may be that all participants have a positive immune response. In this case, there would be positive mass at zero in the immune response distributions amongst SARS-CoV2 seronegative placebo recipients, while amongst SARS-CoV2 seronegative vaccine recipients there would always be a positive immune response. In such cases, one possible path forward is to re-frame the mediation question, for example, by defining S as a binary fold-rise in antibody response.

Beyond these positivity conditions, our model makes no assumption on the distribution P X . However, our results apply to settings where additional assumptions are placed on g A|W and g C , including the possibility that these quantities are known exactly.

We define Y (a, s) as the value Y would assume under a hypothetical intervention that sets A = a, S = s, and C = 1, though our notation suppresses this latter intervention for simplicity. Thus, we can define a counterfactual observation X * = (W, S(a 1 ), Y (a 2 , s) : (a 1 , a 2 ) ∈ {0, 1} 2 , s ∈ S)) ∼ P * X for some set S.

We use E * X {f (X * )} = f (x)dP * X (x) to denote the expectation of any P * X -measurable function f under P * X . The effect of the vaccine can be quantified in terms of a ratio of counterfactual risks, E * X {Y (1, S(1))}/E * X {Y (0, S(0))}. Notice that this is a contrast of average counterfactual outcomes generated under two distinct interventions. The first intervention (in the numerator) randomizes a participant to receive vaccine and then sets their immune response to S(1), that is, the value their immune response would naturally assume under vaccine. The second (in the denominator), similarly randomizes to placebo and sets immune responses to their natural value under placebo. This effect can be decomposed as follows:

where the first term is a natural indirect effect and the latter a natural direct effect. The former compares counterfactual risks under interventions that only differ in how they assign S and thus describes the impact of modulating S on the protective efficacy of the vaccine. The latter compares interventions that only differ in the vaccination status, thereby evaluating the impact of the vaccine that is not influenced by S. For brevity, we hence discuss identification and estimation of ψ * (a 1 , a 2 ) := E * X {Y (a 1 , S(a 2 ))}, for (a 1 , a 2 ) ∈ {0, 1} 2 , noting that each of these effects can be written in terms of these quantities.

Under certain causal assumptions, one can show that ψ * (a 1 , a 2 ) = ψ(a 1 , a 2 ), where, using E X to denote expectation under P X ,

(1)

The form of these assumptions has been the subject of debate; see e.g., Zheng and van der Laan (2012). We omit such a discussion here and focus instead on the statistical estimation problem of the parameter (1).

An estimatorψ(a 1 , a 2 ) of ψ(a 1 , a 2 ) is said to be asymptotically linear ifψ(a 1 , a 2 )−ψ(a 1 , a 2 ) = n −1 n i=1 D(P X )(X i ) + o p (n −1/2 ) for some function D(P X ) such that E X {D(P X )(X)} = 0 and E X {D 2 (P X )(X)} < ∞. The function D(P X ) is called the influence function ofψ(a 1 , a 2 ) and is a key object for (i) characterizing the asymptotic sampling distribution of the estimator; (ii) establishing an efficiency bound for regular estimators of ψ(a 1 , a 2 ); (iii) describing robustness properties of efficient estimators; (iv) constructing efficient estimators.

With respect to the sampling distribution of the estimator, we note that asymptotically linear estimators are particularly convenient for asymptotic analysis, as their large sample behavior is described by the weak law of large numbers and the central limit theorem. The latter implies that n 1/2 {ψ(a 1 , a 2 ) − ψ(a 1 , a 2 )} converges weakly to a random variable with a mean-zero Normal distribution with variance E X {D 2 (P X )(X)}. Influence functions are also useful for establishing local efficiency of estimators. Because the asymptotic variance of an asymptotically linear estimator equals the variance of the its influence function, we can often describe an efficient influence function, that is, the influence function that has the smallest variability amongst all influence functions of regular, asymptotically linear estimators. An estimator with influence function equal to the efficient influence function is said to be asymptotically efficient. The influence function can be used to describe which combinations of parameters of P X must be consistently estimated in order to achieve a consistent estimator of ψ(a 1 , a 2 ). For many parameters, one can obtain a consistent estimator of the parameter of interest with two or more distinct combinations. We describe such parameters and estimators as doubly-or multiply-robust.

Influence functions are also used to construct asymptotically linear and efficient estimators via the one-step estimation framework. Suppose we have an estimate of P X , sayP X , that we use to construct a plug-in estimate of ψ(a 1 , a 2 ), sayψ(a 1 , a 2 ) :=Ê X [Ê X {Ê X (CY | A = a 1 , C = 1, S, W ) | A = a 2 , W }], whereÊ X denotes expectation underP X . A one-step estimator can be constructed asψ + (a 1 , a 2 ) :=ψ(a 1 , a 2 ) + n −1 n i=1 D(P X )(X i ). Under regularity conditions onP X ,ψ + (a 1 , a 2 ) can be shown to be asymptotically linear and efficient. In practice, as we will soon see, it is often not necessary to estimate the entire distribution P X , but rather only key parameters of the distribution, which we refer to as nuisance parameters.

The efficient influence function for ψ(a 1 , a 2 ) is provided in Tchetgen Tchetgen and Shpitser (2012). To write down its form, we first require notation to describe key nuisance parameters. We define g A|W,S (a 0 | W, S) := P X (A = a 0 | W, S), the conditional probability of randomization status given W and S, as well asQ Y (W, S) := E X (Y | A = a 1 , W, S), the conditional probability of the outcome. Subsequently, we can defineQQ Y (W,S) (W ) := E X {Q Y (W, S) | A = a 2 , W } as the conditional mean ofQ Y (W, S) given A = a 2 , W . Finally, we use Q W (w) := P X (W ≤ w) to denote the cumulative distribution function of W . Using these definitions, we can write the nonparametric efficient influence function for ψ(a 1 , a 2 ) evaluated on an observation x as

The efficient influence function is multiply robust in the sense that it is a valid estimating equation for ψ(a 1 , a 2 ) if any of the following conditions holds: (i)Q Y andQQ Y (W,S) are correct; (ii) g A|W,S , g A|W , and g C (1 | a 1 , ·) are correct; or (iii)Q Y , g A|W , and g C (1 | a 1 , ·) are correct. Here we use "correct" to mean equal to the true regression implied by P X . When it comes to estimation, this robustness will imply that we only need to consistently estimate certain combinations of these nuisance parameters to obtain a consistent estimate of ψ(a 1 , a 2 ).

In the context of COVID-19 vaccines, S is measured subject to case-cohort sampling. The observed data can thus be represented as n independent copies of O = (W, A, R, RS, C, CY ) ∼ P , where R is an indicator of having S measured. As with Y , we arbitrarily set S = 0 if immune responses are not measured, though this does not impact subsequent developments. We denote by g R (1 | W, A, C, CY ) = P (R = 1 | W, A, C, CY ) the probability of having immune responses measured. Recall that, because all cases are sampled by design, g R (1 | w, a, 1, 1) = 1 for all w and a, but otherwise will equal to the probability of being randomly selected given covariates and vaccine assignment. The model for P can be defined in terms of the model described above for P X and a model for the sampling mechanism. The only assumption we make on the sampling mechanism is that P {g R (1 | W, A, C, CY ) > 0} = 1, which is satisfied by design in COVID-19 vaccine trials. Rose and van der Laan (2011) provide a convenient characterization of the efficient influence function for pathwise differentiable parameters in two-phase sampling settings. Since ψ(a 1 , a 2 ) is indeed pathwise differentiable and case-cohort sampling is a special case of general two-phase sampling, we can apply these results directly. Their results imply that the efficient influence function of ψ(a 1 , a 2 ) in a nonparametric observed data model for P , when evaluated on a typical observation o can be written

(3) Above, we only need evaluate D(P X ) on observations O such that R = 1 and when R = 1 O = X and thus we are able to evaluate D(P X )(O). Because D(P X ) depends on parameters of the full data distribution, in the next section we establish identifiability of these quantities based on the observed data distribution.

The influence function (3) is doubly robust in the sense that it is a valid estimating equation for ψ(a 1 , a 2 ) if D(P X ) is a valid estimating function in the full data model and either g R or E{D(P X )(O) | R = 1, W = w, A = a, C = c, CY = cy} corresponds to the true value of the regressions implied by P .

A one-step estimator can be constructed directly from the results of Rose and van der Laan (2011). Recall that one-step estimators involve two ingredients: (i) a plug-in estimator of the target parameter and (ii) the empirical average of its efficient influence function at the estimated nuisance parameters evaluated on the observed data.

With respect to (i), we note that a plug-in estimator of ψ(a 1 , a 2 ) can be obtained via inverse probability weighted sequential regression. First, we require an estimate g n,R of g R , if it is unknown. This could be obtained using regression of the binary outcome R onto W, A, C, and CY . The estimated sampling probabilities are then used to estimateQ Y . Here, we regress the outcome Y onto covariates W and S amongst the subset of data with R = 1 and A = a 1 , while including inverse probability weights R/g n,R (1 | W, A, C, CY ) into the learning procedure. For example, logistic regression working models could be used with regression parameters estimated by maximizing an inverse weighted log-likelihood. More flexible learning approaches could also be adopted. Irrespective of the learning approach chosen, we letQ n,Y denote the estimated regression function. The second step in the procedure involves evaluatingQ n,Y (W i , S i ) for i such that R i = 1. This prediction serves as the outcome in a regression onto covariates W amongst the subset of data with R = 1, A = a 2 , again including inverse probability weights. LetQ n,Q Y (W,S) denote the estimated regression. A plug-in estimate is obtained as ψ n,1 (a 1 , a 2 ) = n −1 n i=1Qn,QY (W,S) (W i ). Inspection of (3) reveals that it involves both the full data efficient influence function (2) (in the first line), as well as an additional nuisance parameterQ D(P X ) (W, A, C, CY ) := E{D(P X )(O) | R = 1, W, A, C, CY } (second line). To evaluate the full data efficient influence function, we will need estimates of g A|W , g C , and g A|W,S . The randomization probability g A|W is known our example, but in general could be estimated. To maintain generality, we denote by g n,A|W the estimated regression function. An estimate g n,C of the censoring probability g C could be obtained via regression of the binary outcome C onto W amongst observations with A = a 1 . Importantly these two regression, do not require inverse probability weights, since the distribution of A, W, and C are observed for everyone. On the other hand, an estimate g n,A|W,S of g A|W,S would need to involve inverse probability weights, for example as in a weighted regression of the binary outcome A onto W and S.

Once all of the above regressions have been estimated, we can replace the true nuisance parameters in (2) with their estimated counterparts and evaluate the expression for all O i such that R i = 1,

Next, we obtain an estimateQ n,D(P X ) of the additional nuisance parameterQ D(P X ) by regressing the pseudo-outcomeD

The one-step estimator is ψ + n,1 (a 1 , a 2 ) := ψ n,1 (a 1 , a 2 ) + n −1 n i=1D 1 (O i ). We have the following theorem describing the behavior of the one-step estimator. We write the squared L 2 (P )-norm of a square-integrable function f as ||f || 2

Theorem 1. Assume the following conditions:

andD 1 falls in a P -Donsker class with probability tending to 1.

Under these assumptions n 1/2 {ψ + n,1 (a 1 , a 2 )−ψ(a 1 , a 2 )} converges in distribution to a random variable with a Normal(0, E{D 2 (P )(O)}) distribution.

We provide a proof in the web supplement. The uniform consistency assumptions in A1-A3 ensure estimates of treatment and censoring probabilities are bounded away from zero. The L 2 rate conditions in A1-A5 ensure negligibility of a second-order remainder term. The n −1/4 convergence rates are sufficient for the results of the theorem, but could also be weakened to conditions involving the product of rates for various combinations of nuisance parameters. We note that the n −1/4 rate is slower than the standard parametric rate of n −1/2 , implying that flexible regression techniques could be used. However, due to the curse of dimensionality, in settings with even moderately large dimensions of S or W , the n −1/4 rates may be difficult to achieve without further smoothness assumptions. For example, the highly adaptive lasso assumes that the nuisance functions have variation norm bounded by a constant and achieves requisite L 2 convergence irrespective of the dimension of the regression (Benkeser and van der Laan, 2016). However, this bounded variation condition can be rather stringent in higher dimensions. Assumption A6 is needed to ensure that an empirical process term is negligible (van der Vaart and Wellner, 1996) , but could be removed by using cross-fitting as in Zheng and van der Laan (2011) and Chernozhukov et al. (2017) .

The robustness of the full-and observed-data efficient influence functions imply that ψ + n,1 (a 1 , a 2 ) is consistent if: (i) either (i.a)Q n,Y andQ n,Q Y (W,S) are consistent for their respective targets or (i.b) g n,A|W,S , g n,A|W , and g n,C (1 | a 1 , ·) are consistent or (i.c)Q n,Y , g n,A|W , and g n,C (1 | a 1 , ·) are consistent and (ii) either (ii.a) g n,R is consistent for g R or (ii.b)Q n,D(P X ) is consistent forQ D(P X ) .

Under the assumptions of the theorem, we also have that σ 2 n,

is a consistent estimate of the asymptotic variance of n 1/2 ψ + n,1 (a 1 , a 2 ) that can be used to construct Wald-type confidence intervals.

Our second estimator differs both in terms of the construction of the plug-in estimator, as well as in the form of the efficient influence function used in the one-step estimator. Both changes are the result of a more explicit examination ofQ D(P X ) in (3). We definẽ

which is similar toQ D(P X ) , but only involves the first term in D(P X ). We also definẽ

In the appendix, we show thatQQ (W,A,C,CY ) = QQ Y (W,S) defined above.

Using the above-defined nuisance parameters, we can rewrite (3) as follows,

To construct a one-step estimator based on this representation, we describe construction of a plug-in estimate and an approach for evaluating the one-step correction term.

For the plug-in estimator, we begin as above by obtainingQ n,Y , an estimate ofQ Y using inverse weighted regression and evaluatingQ n,Y (W i , S i ) for all i such that R i = 1. We then use these estimates as the outcome in a regression onto W, A, C, CY amongst those with R = 1, which yields an estimateQ n,

. . , n and use this as the outcome in a regression onto W amongst those with A = a 2 , thereby providing an estimateQ n,Q(W,A,C,CY ) of QQ (W,A,C,CY ) . Finally, the plug-in estimator is ψ n,2 (a 1 , a 2 ) := n −1 n i=1Qn,Q(W,A,C,CY ) (W i ). With the estimates of the sequential outcome regression described above and the estimates of and g R , g A|W , g A|W,S described above, the only additional nuisance parameter estimate we require is one ofQ D . Such an estimate can be obtained by evaluatinĝ

for all i such that R i = 1. The quantityD X,1 (O i ) then serves as an outcome in a regression

We denote byD 2 the efficient influence function in (5), but with true nuisance parameters substituted with their estimated counterparts, just as we did in (4) for the classic one-step estimator. Our proposed one-step estimator is ψ + n,2 (a 1 , a 2 ) = ψ n,2 (a 1 , a 2 ) + n −1 n i=1D 2 (O i ). We have the following theorem describing its behavior. LetQ Under these assumptions n 1/2 {ψ n,2 (a 1 , a 2 )−ψ(a 1 , a 2 )} converges in distribution to a random variable with a Normal(0, E{D 2 (P )(O)}) distribution.

The conditions are similar to those of Theorem 1. The estimator ψ n,2 (a 1 , a 2 ) is also multiply robust in the sense of the classic one-step estimator, but where we replace (i.a) with condition (i.a')Q n,Y ,Q n,Q(W,A,C,CY ) are consistent for their respective targets and replace (ii.b) with condition (ii.b')Q n,D(P X ) is consistent forQ D(P X ) andQ n,Q Y (W,S) is consistent forQQ Y (W,S) . Under the assumptions of the theorem, σ 2 n,2 = n −1 n i=1 {D 2 (O i )− n −1 n j=1D 2 (O j )} 2 is a consistent estimate of the asymptotic variance of n 1/2 ψ n,2 (a 1 , a 2 ). We provide a detailed theoretical comparison of the two one-step estimators in the supplementary materials.

We turn to a discussion of the implications of our theoretical results on planning for the analysis of randomized trials of COVID vaccines. First, since we are in the context of a randomized trial, g A|W is known exactly. Thus, the conditions pertaining to consistent estimation of this quantity are easily satisfied. For example, we could use the known randomization probabilities or fit a low-dimensional parametric regression to adjust for chance imbalances in covariates. Similarly, the sampling probability g R is known by design. If we use the known sampling probability in the estimation procedure, then we can use any estimator ofQ D(P X ) (or, for the alternative one-step,Q D(P X ) ) and the results of the theorem hold irrespective of the consistency ofQ n,D (Q n,D ). All that we need is that these estimates converge to something at a reasonable rate. Second, we comment on estimation of regression quantities that involve conditioning on C and CY . Note that there are three categories into which a participant's data may fall, right-censored (C = 0, CY = 0), not right-censored and not a case (C = 1, CY = 0), and not right-censored and a case (C = 1, CY = 1). Thus, for regression quantities that require conditioning on these data, we can construct dummy variables for two of these cases to include in the regression.

Third, in the above developments, we have left unspecified which regression techniques may be used to estimate the nuisance parameters. We provide a software implementation that implements two specific approaches, but also that allows users to define their own regression approaches. Specifically, we adopt approaches based on working linear or logistic (depending on the parameter space of the particular nuisance parameter) regression models, as well as based on super learner (van der Laan et al., 2007) . In the latter approach, one defines a pre-specified library of candidate regression estimators. The resultant regression estimator is a convex combination of the library of candidate estimators where the weights are selected by minimizing a cross-validated risk criterion.

The goal of this simulation was to verify the statistical properties of our estimators as established by our theorems. Accordingly, we used a discrete data generating process so that it was possible to (i) numerically approximate the efficiency bound and (ii) use nonparametric maximum likelihood estimators for nuisance parameters thereby guaranteeing that requisite regularity conditions required by our theorems were satisfied. We generated data as follows. A bivariate variable W = (W 1 , W 2 ) was generated by independently drawing two Bernoulli(1/2) variables. Given W = w, a binary treatment A was generated according to a Bernoulli distribution with g A|W (1 | w) = expit(w 1 − w 2 ). Given W = w and A = a, S was drawn from a Binomial(2, p(a, w)) distribution with p(a, w) = expit(−1+w 1 /4−w 2 /3+a/2). Given W = w, A = a, S = s, a binary outcome Y was drawn from a Bernoulli distribution with P (Y = 1 | W = w, A = a, S = s) = expit(−2 + a/2 + w 1 /2 − s/2), while C was drawn from a Bernoulli distribution with g C (1 | a, w) = expit(2 + w 1 /2 − w 2 /3). Subsequently, case-cohort sampling was applied to S, by sampling a random 1/4 of the cohort, in addition to all observations with CY = 1. We evaluated our two proposed estimators of ψ(1, 0), the true value of which was approximately 0.187 with efficient variance bound equal to E{D 2 (P )(O)} = 0.509.

For each sample size n ∈ {500, 1000, 2000, 4000, 8000}, we simulated 1000 data sets from this data generating process and computed our two proposed estimators along with a confidence interval for ψ(1, 0). We evaluated the estimators in terms of their bias (scaled by n 1/2 ), their standard error (scaled by n 1/2 ), the coverage probability of a nominal 95% confidence interval, and the ratio of the scaled standard error to the square root of the efficient variance.

We first evaluated estimators under the conditions of the theorem where all nuisance parameters are consistently estimated at appropriate rates. To achieve this, we used parsimonious logistic regression models for g A|W , g C , g R ,Q Y and fully saturated logistic regression models (i.e., nonparametric maximum likelihood estimates) for the remaining nuisance parameters. In this setting, our theory predicts that bias of the estimators should converge to zero faster than n −1/2 , the scaled standard error of the estimators should converge to the square root of the efficiency bound, and confidence intervals should attain nominal coverage.

After confirming asymptotic properties in the setting with consistently estimated nuisance parameters, we subsequently studied six scenarios in which various combinations of nuisance parameters were misspecified. The six scenarios that we studied are ones in which our multiple robustness results imply that the one-step estimators should remain consistent. In these settings, certain nuisance parameters were inconsistently estimated by using a simple intercept-only regression model (i.e., the sample average of the outcome of the regression). In these settings, our results suggest that we should still see diminishing bias and stabilizing scaled standard error as sample size increases. However, neither estimators' influence function will equal the efficient influence function. Consequently, we expect non-nominal coverage probabilities for the confidence intervals.

The results of the simulation results confirmed our theory and estimators enjoyed expected behavior in all settings considered (Table 1 ). In the setting where all nuisance parameters were consistently estimated, the two one-step estimators had similar performance across the board. However, in settings where the regression of the full data efficient influence function was inconsistently estimated, the alternative one-step estimator ψ + n,2 (1, 0) tended to have higher variability.

In our second simulation, we wished to study the finite-sample performance of our estimators in a setting similar to what might be expected in a correlates study of a COVID-19 vaccine. To that end, we developed a data generating process meant to mimic the expected data outputs of a large (n = 30000), placebo-controlled trial. We simulated a threedimensional covariate W = (W 1 , W 2 , W 3 ) with three Bernoulli-distributed components with success probabilities of 2/5, 1/4, 1/4, meant to represent measurements of participants' age (≥ 65 vs. < 65 years old), racial/ethnic minority status (yes vs. no), and presence of COVID-19 risk factors (yes vs. no). For each simulated participant, we drew their randomization assignment A from a Bernoulli(1/2) distribution. A continuous immune marker was then generated as follows. Given W 1 = w 1 , we drew S * from a Normal(2 -0.5 w 1 , 1) distribution. Next, given A = a, we set S = a1(S * > 0)S * . That is, placebo recipients uniformly have S = 0, which is reasonable as the primary analysis of COVID-19 trials focuses on previously uninfected individuals. For individuals in the vaccine arm, there are some "vaccine failures," that is individuals in whom the vaccine fails to induce a positive immune response. Under our simulation setting, the probability of failure increases with age. Given W = w, S = s, A = a we set the probability of having an observed COVID-19 disease endpoint as expit(α − 0.5s − 1.8a + 0.2w 1 + 0.1w 2 + 0.7w 3 ). Note that the parameter α controls the total number of breakthrough endpoints. We considered performance of our estimators under several values of α that may be realistic for COVID vaccine efficacy studies. Under this data generating process, the vaccine has high efficacy, as was observed with the first reported COVID-19 vaccines. For simplicity, we assumed there was no censoring of endpoints. The corresponding true parameter values for each setting we considered is shown in Table 2 . Two-phase sampling of S was simulated according to the sampling plan developed for all Phase 3 trials run in the United States (Gilbert et al., 2021) . In this design, a stratified random sample is drawn based on the sixteen strata defined by W and A. Specifically, 113 vaccine recipients and 15 placebo recipients are sampled from each of the eight strata defined by W . All observed COVID-19 endpoints are additionally sampled.

We were interested primarily in comparing the two one-step estimators in terms of their performance at realistic sample sizes in terms of point estimation and confidence interval coverage for the indirect effect ψ(1, 1)/ψ(1, 0) and for the proportion mediated, defined as 1−log{ψ(1, 0)/ψ(0, 0)}/log{ψ(1, 1)/ψ(0, 0)}, as these are two key parameters for establishing vaccine correlates. For each setting described above, we simulated 1000 data sets and computed the two one-step point estimates and confidence intervals.

We repeated the simulation under several different modeling strategies for nuisance parameters. The first considered regression models that included all possible interactions between regressors, where S was modeled as a linear term. In this approach, we do not expect that the nuisance parameters are consistently estimated, though the bias from model misspecification is expected to be relatively small. Our second strategy for nuisance parameter estimation considered using the super learner. In our analysis, we used four different regressions that included a main-terms generalized linear model, a stepwise generalized linear model that included all two-way interaction terms, an intercept-only model, and polynomial multivariate adaptive regression splines. Due to the added flexibility of the ensemble model, we believe this approach should come closer to consistent estimation of the nuisance parameters. However, it is a highly data-driven and so we wished to evaluate whether it could be appropriately employed at these sample sizes. Finally, we considered using simple main-terms regression models that did not include any interactions between variables. These models were expected to have the worst performance in modeling the nuisance parameters and served as a "worst-case scenario" in our simulation. We report on the generalized linear models with interaction and super learner results in the body and describe the misspecified main terms models in the supplement.

For both one-step estimators and irrespective of the nuisance parameter estimation approach, we found that when events were very rare in the vaccine group, the estimators exhibited considerable bias. In these rare-event settings, confidence intervals for the indirect effect and proportion mediated had slightly less than nominal coverage. However, as the number of expected events amongst the vaccinated rose to around 50, the bias of the estimators decreased to a reasonable level and the confidence intervals attained nominal coverage probabilities. Comparing the two modeling strategies, we found comparable results between the GLM with interactions and the SuperLearner. Comparing the two one-step estimators, we saw similar performance in terms of bias and confidence interval coverage. However, we found that confidence intervals for the standard one-step estimator when combined with super learner tended to be more erratic (Figure 1 ). The same phenomenon was not observed when the GLM with interactions was used to estimate nuisance parameters (Figure 2 ). Figure 1 : Indirect effect confidence intervals using classic (left column) and alternative (right column) one-step estimators constructed using super learner to estimate nuisance quantities. The confidence intervals are ordered across the 1000 simulations from smallest to largest and are displayed on the log scale. In the rare event setting (α = −5), we include a subfigure that is zoomed in and has the extreme confidence intervals removed. Figure 2 : Indirect effect confidence intervals using classic (left column) and alternative (right column) one-step estimators constructed using GLMs with interactions to estimate nuisance quantities. The confidence intervals are ordered across the 1000 simulations from smallest to largest and are displayed on the log scale. In the rare event setting (α = −5), we include a subfigure that is zoomed in and has the extreme confidence intervals removed.

We found that basing inference on the misspecified main terms models tended to lead to estimators with non-negligible bias and conservative inference, with confidence interval coverage probabilities near to 1 (Supplemental material). However, we did find more robust performance of the alternative one-step estimator relative to its standard counterpart.

Overall, we found that the estimators performed reasonably well in all but the smallest sample sizes considered. We also found that using super learning to estimate nuisance parameters appears acceptable even in settings with a limited number of observed events.

The supplemental material contains additional simulation results for operating characteristics of the estimators of the counterfactual risk parameters ψ(a 1 , a 2 ). Of note, we show there that the classic one-step estimator ψ + n,1 (1, 0) < 0 for a non-negligible fraction of simulations. The results are particularly bad for the misspecified nuisance models. However, even with super learner ψ + n,1 (1, 0) was less than zero up to 5% of simulated data sets. On the other hand, the alternative one-step estimator was always positive.

A central contribution of this study is establishing that nonparametric estimators of mediation effects can reasonably be employed in the analysis of COVID-19 correlates data assuming a reasonable number (≈ 50) of vaccine breakthrough cases are observed. The fact that super learner can be used to estimate nuisance parameters while retaining reasonable operating characteristics is highly appealing given (i) the inherent complexity of some of the nuisance parameters that need to be estimated as part of our procedure (e.g.,Q D(P X ) ) and (ii) our limited understanding of the interaction between patient-level COVID-19 risk factors and vaccine-induced immune responses on risk of COVID-19.

Beyond this practical contribution, we believe that the alternative approach outlined to one-step estimation is also a valuable contribution to the literature on estimation of pathwise differentiable parameters in two-phase sampling designs. An attractive feature of this alternative approach is the decoupling of estimation of the sampling probabilities and estimation of full-data nuisance parameters. This should provide more practical avenues for achieving a robust estimator in settings where sampling probabilities are unknown. However, the results of the simulations were mixed in terms of the performance of this alternative onestep estimator vs. the standard one. In our COVID-19-oriented simulation, we saw more stable inference in small samples using the alternative one-step. The alternative estimator of ψ(1, 0) never slipped below 0, while this occurred with some frequency for the classic estimator. However, in the large sample simulation, we saw that the alternative one-step may not be as robust to certain forms of nuisance parameter misspecification. A careful comparison of these estimation approaches across a more diverse set of scenarios is warranted and is an important direction for future research.

Another potential area for expanding this work area is in developing a targeted minimum loss-based estimation (TMLE) framework for estimation of these mediation parameters. As substitution estimators, TMLEs enjoy the property of being guaranteed to fall in the parameter space. This property may be particularly important in the context of studies of highly effective COVID-19 vaccines where risks of disease amongst the vaccinated are extremely small. Very modest one-step correction terms can easily push estimators below zero, as we saw in the simulation study. TMLEs would likely enjoy more robust behavior in these settings.

The code needed to reproduce the both simulation studies is available on in a GitHub repository at https://github.com/benkeser/natmed2 sims. 

The R package natmed2 is available on Github at https://github.com/benkeser/natmed2. Web Appendices referenced appear below. Appendix A includes a proof of Theorems 1 and 2; Appendix B contains theoretical comparison of the classic and alternative one-step approaches; Appendix C provides an alternative identification result that does not rely at all on full-data nuisance parameters; and Appendix D contains additional simulation results.

For proving asymptotic linearity of the proposed estimators, it will first be useful to establish the following lemmas. We recall the following definitions:

The first line follows by definition ofQ Y , the second by the tower rule. The third follows from the fact that two-phase sampling status R is independent ofQ Y (W, S) given (W, A, C, CY ). The fourth follows again by definition, the fifth recognizing that the conditional distribution of (A, C, CY ) given A = a 2 , W implied by P X is the same as that implied by P . The final line follows by definition.

Lemma 2. Let M F denote the non-parametric full-data model described in the main document. Let Ψ F : M F → R denote the parameter mapping that maps any distribution P X in the model to a number

where E X denotes expectation under P X . The efficient influence function D(P X ) of Ψ in the full data model satisfies

Proof. This lemma follows from Theorem 4 in the supplementary materials of Díaz et al. (2020) .

Lemma 3. Let V = (W, A, C, CY ), and let P denote the observed data distribution implied by coarsening P X according to g R . The efficient influence function D(P ) of Ψ in the observed data model satisfies

where R 2 (P X , P X , P , P ) = R F 2 (P X , P X )

Proof. This follows from Lemma 2 and the definition of D(P ) after computing

using the assumption R ⊥ ⊥ X | V .

Corollary 1. The object R F 2,1 can equivalently be written as

Proof. The equality of R F 2,1 follows immediately from Lemma 1. The equality of R 2 comes from an algebraic simplification of (6) in Lemma 3.

We are now ready to prove the theorems. We start with the proof of asymptotic linearity for the classic one-step estimator. Below, we useP andP X denote any distribution in the observed and full data models, respectively, that is compatible with our estimated nuisance parameters. We define the shorthand notation g j (W ) := g A|W (a j | W ) and g n,j (W ) := g n,A|W (a j | W ). Similarly, we defineg j (W, S) := g A|W,S (a j | W, S), g n,j (W, S) := g n,A|W,S (a j | W, S), and g C,j (W ) := g C (1 | a j , W ), g n,C,j (W ) := g n,C (1 | a j , W ). We also make use of the shorthand notation P f to denote expectation of a P -measurable function f of O under a given probability distribution P ; thus, P f = f (o)dP (o). In particular, we denote by P n the empirical distribution based on O 1 , . . . , O n so that P n f = n −1 n i=1 f (O i ). Note that based on these definitions, the classic one-step estimator can be written as ψ + n,1 (a 1 , a 2 ) = Ψ F (P X ) + P n D(P X ). Thus, by Lemma 3, ψ + n,1 (a 1 , a 2 ) − ψ(a 1 , a 2 ) = (P n − P )D(P ) + R 2 (P X , P X ,P , P ) = (P n − P )D(P ) + R 2 (P X , P X ,P , P ) + o p (n −1/2 ) = P n D(P ) + R 2 (P X , P X ,P , P ) + o p (n −1/2 ) which follows from assumption A6 and the fact that D(P ) has mean zero. Thus, it remains to establish that R 2 (P X , P X ,P , P ) = o p (n −1/2 ). We note that

where the last line follows from assumptions A1 and A5. Next, we note that R F 2,2 (P X , P X ) can be written as P X g 2 g n,2

(g 1 −g n,1 ) + (g n,2 −g 2 ) +g 2 (g 1 −g n,1 ) +g n,1 (g n,2 −g 2 ) g n,1g1 (Q Y −Q n,Y (W, S)) .

Thus, R F 2,2 (P X , P X ) can be split into four terms that each may be analyzed separately in a similar fashion as with R F 2,1 . For example, we have P X g 2 g n,2

which follows from Assumptions A3 and A4. Likewise, we can bound R F 2,3 (P X , P X ),

as well as R 2 ,

Thus, ψ + n,1 (a 1 , a 2 ) − ψ(a 1 , a 2 ) = P n D(P ) + o p (n −1/2 ). Because O 1 , . . . , O n are independent, the central limit then implies the result.

The proof of Theorem 2 follows almost exactly as above. We start by noting that based on assumption B7, we have ψ + n,2 (a 1 , a 2 ) − ψ(a 1 , a 2 ) = P n D(P ) + R 2 (P X , P X ,P , P ) + o p (n −1/2 ) .

All that remains is to establish the negligibility of the remainder terms as written in Corollary 1. This proceeds similarly as above and we can establish that under the conditions of the theorem R F 2,1 (P X , P X ) = P X g −1 n,2 (Q n,Q(W,A,C,CY ) −QQ (W,A,C,CY ) ) (g n,2 − g 2 )

as well as that R 2,1 (P X , P X ,P , P ) = P g −1 n,

The main difference between the classic and alternative approach is that the latter does not involve inverse probability weights when integratingQ Y (W, S) over S. Instead, we first integrate out S over the conditional distribution of S given R = 1, W, A, C, CY . By assumption, this conditional distribution is is equal to the conditional distribution of S given W, A, C, CY in the full data distribution, so there is no need to bother with inverse probability weights in its estimation. A second regression is then applied to integrate out A, C, and CY to obtain the appropriate function of W . This strategy may have a theoretical advantage in situations where g R is unknown. To understand why, note that the multiple robustness of the classic approach allows for inconsistent estimation of g R , so long asQ n,D is consistent. However, estimation of the sequential regressionsQ Y andQQ Y relies on inverse weighting by the estimate of g R . Thus, if g R is inconsistently estimated, it is likely that these quantities will be inconsistently estimated as well. In this case, consistency of the estimator is fully reliant on consistent estimation of the other nuisance quantities. On the other hand, the alternative one-step estimator decreases the reliance of the estimator on g R by replacing one step of inverse weighted estimation. Thus, we hypothesize that the alternative one-step may enjoy better performance in situations where g R is inconsistently estimated. However, this is not fully satisfactory, as estimation of the sequential regression quantityQQ Y still relies on estimation ofQ Y , which itself relies on estimation of g R . Thus, the alternative estimator may also be subject to poor performance when g R is inconsistently estimated.

7.1 C. A one-step estimator that does not rely on inverse-weighted learning

In fact, we can outline an estimation approach that fully disentangles estimation of g R and estimation of the required full data nuisance parameters. However, the approach does not appear to scale well in settings where S is continuous valued. Let q S (s | W, A, Y ) be defined as the conditional density of S given R = 1, W, A, C = 1, Y evaluated at s. Note that this parameter of P is the same as the conditional density of S given W, A, C = 1, Y under P X . For each s ∈ S, we define the nuisance quantities

We now derive a relationship between E X (CY | C = 1, W, A = a 1 , S) and these two quantities. For simplicity and without loss of generality we show the derivation when data are discrete.

Thus, we have the following identification result,

This suggests the following estimation strategy for evaluating the plug-in estimator. First, obtain an estimate q n,S of the conditional mediator density q S . For all i such that R i = 1, evaluate q n,S (S i | W = W j , A = a 1 , Y = Y j ) for all j such that C j = 1 and A j = a 1 .

Next, obtain an estimateQ n,Y q S (S i | W ) ofQ Y q S (S i | W ) by using observations such that C j = 1 and A j = a 1 to regress the outcome Y j × q n,S (S i | W = W j , A = a 1 , Y = Y j ) onto W j . Similarly, obtain an estimateQ n,q S (S i | W ) ofQ q S (S i | W ) by regressing the outcome q n,S (S i | W = W j , A = a 1 , Y = Y j ) onto W j . Next, for all i such that R i = 1, evaluateQ n,Y q S (S i | W i )/Q n,q S (S i | W i ). This quantity serves as our estimate ofQ n,Y and we proceed as with the plug-in estimator for the alternative one-step.

The key element of this estimator's construction is that it does not rely on inverse weighting by an estimate of g R to estimate nuisance parameters. However, there is a price to be paid in terms of complexity. Whereas the previous estimators each allowed us to avoid estimation of the conditional mediator density, in this approach we require such an estimator. Moreover, we require a separate estimation ofQ q S andQ q S for each value at each value S i such that R i = 1. Thus, if S is continuous-valued, this would require fitting 2 n i=1 R i regressions. This strategy therefore appears most-appealing in situations where S is binary, in which case its implementation is quite tractable, both in terms of density estimation and in terms of the estimation ofQ q S andQ q S .

We found that the estimates based on a misspecified GLM performed quite poorly in terms of their bias (Table 4) . Not only that, the variance of the estimators was extremely inflated and confidence intervals tended to have conservative coverage probabilities due to standard error estimates that were too large. Table 5 shows estimators of the counterfactual risk ψ(1, 0). We found similar results as with the estimators of indirect effect and proportion mediated. Confidence intervals attained near nominal coverage in larger samples with negligible bias. Table 1 : Bias (scaled by n 1/2 ), standard error ((scaled by n 1/2 , SE), coverage probability of a nominal 95% confidence interval (Cov), and the scaled SE divided by the square root of the efficient variance (Ratio) for the two proposed one-step estimators in different settings. The settings column indicates which nuisance parameters are consistently estimated (g R was consistently estimated in all settings). ψ + n,1 (1, 0) ψ + n,2 (1, 0) Setting Table 5 : Performance of the two one-step estimators of ψ(1, 0) in terms of bias, coverage probability of a 95% confidence interval, and the median estimated standard error divided by the true standard error.

ψ + n,1 (1, 0) ψ + n,2 (1, 0) α ψ(1, 0) Bias Coverage

Est. std. True. std.

Est. std. True. std. 

Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine

The highly adaptive lasso estimator

Using the whole cohort in the analysis of case-control data

Double/debiased/neyman machine learning of treatment effects

A strategic approach to COVID-19 vaccine R&D

Influenza hemagglutination-inhibition antibody titer as a mediator of vaccine-induced protection for influenza B

Non-parametric efficient causal mediation with intermediate confounders

OWS CoVPN COVID-19 vaccine efficacy trial immune correlates statistical analysis plan

Evaluating candidate principal surrogate endpoints

Identification, inference and sensitivity analysis for causal mediation effects

Related causal frameworks for surrogate outcomes

A bayesian approach to surrogacy assessment using principal stratification in clinical trials

Clarifying causal mediation analysis for the applied researcher: Defining effects based on what we want to learn

Direct and indirect effects

Estimation of direct causal effects

Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine

A framework for assessing immunological correlates of protection in vaccine trials

Identifiability and exchangeability for direct and indirect effects

A targeted maximum likelihood estimator for two-stage designs

Semiparametric theory for causal mediation analysis: efficiency bounds, multiple robustness, and sensitivity analysis