{'fraction_non_alphanumeric': 0.03810291068045175, 'fraction_numerical': 0.017007563890256455, 'mean_word_length': 4.699345293767537, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient nonparametric inference on the effects of stochastic interventions under two-phase sampling, with applications to vaccine efficacy trials

The advent and subsequent widespread availability of preventive vaccines has altered the course of public health over the past century. Despite this success, effective vaccines to prevent many high-burden diseases, including HIV, have been slow to develop. Vaccine development can be aided by the identification of immune response markers that serve as effective surrogates for clinically significant infection or disease endpoints. However, measuring immune response is often costly, which has motivated the usage of two-phase sampling for immune response sampling in clinical trials of preventive vaccines. In such trials, measurement of immunological markers is performed on a subset of trial participants, where enrollment in this second phase is potentially contingent on the observed study outcome and other participant-level information. We propose nonparametric methodology for efficiently estimating a counterfactual parameter that quantifies the impact of a given immune response marker on the subsequent probability of infection. Along the way, we fill in a theoretical gap pertaining to the asymptotic behavior of nonparametric efficient estimators in the context of two-phase sampling, including a multiple robustness property enjoyed by our estimators. Techniques for constructing confidence intervals and hypothesis tests are presented, and an open source software implementation of the methodology, the txshift R package, is introduced. We illustrate the proposed techniques using data from a recent preventive HIV vaccine efficacy trial.


Introduction
Ascertaining the population-level causal effects of exposures is often a motivating goal in scientific research. Such effects are commonly formulated via summaries of the distribution of counterfactual random variables, which describe arXiv:2003.13771v1 [stat.ME] 30 Mar 2020 the values a measurement would have taken if, possibly counter-to-fact, a particular level of exposure were assigned to the unit. Often, the exposure of interest is continuous-valued -for example, the dose of a drug, amount of weekly exercise, or level of an immune response marker induced by a vaccine. The latter serves as our motivating example as we consider data generated by a phase IIb trial of a vaccine to prevent infection by human immunodeficiency virus (HIV), the HIV Vaccine Trials Network's (HVTN) 505 efficacy trial (Hammer et al., 2013). In addition to evaluating the overall efficacy of the vaccine, a key secondary question of the trial was to evaluate the role of vaccine-induced immune responses in generating protective efficacy against HIV (Janes et al., 2017). Identification of immune response markers causally related to protection is critical both for further developing a biological understanding of the action mechanism of a vaccine's effects and for insights to guide the development of future vaccines.
To study such relationships, it is natural to consider a dose-response curve that summarizes participants' risk of HIV infection as a function of the level of a particular immune response marker. A causal formulation of such a doseresponse analysis would consider a (possibly infinite) collection of counterfactual outcomes, each representing the HIV infection risk that would have been observed if all individuals' immune responses had been set to a particular level. Studying how the proportion of infected individuals varies as a function of the level of an immune response marker could provide insights into causal mechanisms underlying the vaccine's effects. Unfortunately, several difficulties, both theoretical and practical, arise when considering such a dose-response approach. From a statistical perspective, nonparametric estimation and inference on the causal dose-response curve is challenging and requires non-standard techniques (Díaz and van der Laan, 2013;Kennedy et al., 2017;van der Laan et al., 2018). More importantly, such an approach may require consideration of counterfactual variables that are scientifically unrealistic. Namely, it may be impossible to imagine a world where every participant exhibits high immune responses, simply due to phenotypic variability of participants' immune systems. This fact calls into question both the validity and practicality of the counterfactual dose-response analysis strategy in the context of evaluating the effects of immune responses.
An alternative framework for assessing the causal effects of continuous-valued exposures is rooted in considering counterfactual outcomes resulting from stochastic interventions (Díaz and van der Laan, 2012;Haneuse and Rotnitzky, 2013;VanderWeele and Hernan, 2013;Young et al., 2014). Whereas static interventions assign the same fixed level of an exposure to all observed units, stochastic interventions consider instead setting the exposure level equal to a random draw from a particular distribution. This approach provides a more flexible means of defining counterfactual random variables. Indeed, static interventions can be viewed as a special case of stochastic interventions in which the intervention mechanism is drawn from a degenerate distribution with all mass placed on only a single exposure level. In order to define scientifically meaningful counterfactuals, care must be taken in defining the particular distribution from which the exposure is drawn. A popular strategy is to draw exposure levels from a modified version of the true exposure distribution -the natural distribution of the exposure that occurs under no intervention. For example, one may consider an intervention that draws the post-vaccination level of immune response markers from a distribution that is similar to the naturally observed (post-vaccination) distribution of immune response markers but that has been shifted upwards (or downwards) for some or all participants. Counterfactuals defined by such an intervention may be better aligned with a plausible future intervention, such as a refined version of the current vaccine that provides improved immune response in some or all participants. Evaluating the population-level risk of HIV infection under such interventions is scientifically useful for several reasons. Firstly, this measure of risk provides a scientifically relevant mechanism by which to rank-order immune response markers in terms of their importance with respect to HIV risk. Such information could be used in defining go/no-go criteria in future early-phase HIV vaccine development pipelines. Secondly, the risk measure may also provide a way to predict next-generation vaccine efficacy based on the induced immune response profile, elucidating whether the immune response markers induced by a candidate vaccine are sufficiently promising to advance the vaccine to the clinical trial stage.
The above discussion highlights the need for rigorous methodology to identify and estimate population-level causal effects of stochastic interventions; recent work has provided several approaches (Díaz and van der Laan, 2012;Haneuse and Rotnitzky, 2013;Young et al., 2014;Kennedy, 2019). Both Díaz and van der Laan (2012) and Haneuse and Rotnitzky (2013) propose conditions for the identification of the causal effect under a stochastic intervention and detail several estimators of these quantities, including strategies relying on inverse probability weighting, outcome regression, and doubly robust estimation. These approaches alone are insufficient for application to studies like the HVTN 505 trial, where a two-phased, case-control sampling design was used to measure participants' immune response profiles. Under this design, all participants were HIV-negative at the week 24 visit, with eligible cases being vaccine recipients diagnosed with HIV-1 infection by the month 24 visit; 100% of cases were sampled from among those with samples available for measurement of immune response markers at the week 26 visit, while a random sample of eligible HIV-uninfected vaccine recipients was taken (Janes et al., 2017). This sampling design severely complicates the estimation of causal effects. Rose and van der Laan (2011), among others, discuss strategies for efficient estimation under two-phase sampling designs, placing emphasis on an inverse probability of censoring weighted modification that may be coupled with targeted minimum loss estimation (or similar frameworks) to account for the study design.
Their approach is shown to yield an asymptotically linear estimator so long as the probability of inclusion in the second-phase sample is known by design or estimable via nonparametric maximum likelihood. This latter requirement precludes usage of their proposed estimators in situations where, for example, sampling probabilities depend on continuous-valued covariates. While the literature has tended to use "two-phase sampling" to denote outcomedependent Bernoulli or without-replacement sampling based on discrete covariates, recent efforts have extended the concept to allow the use of continuous-valued covariates in construction of the second-phase sample (e.g., Chatterjee et al., 2003;Chatterjee and Chen, 2007;Gilbert et al., 2014). Consequently, Rose and van der Laan (2011) sketched a more complicated procedure for generating efficient estimators in such settings, one requiring augmentation of their simple reweighting scheme. To our knowledge, this approach has neither been evaluated via simulation nor in applied data analysis.
In the present work, we develop estimators of the mean counterfactual outcome under a stochastic intervention when the exposure is measured via two-phase sampling. We provide several contributions to the disparate literatures on two-phase sampling and stochastic interventions. To the former, we (i) formalize the assumptions needed for efficient nonparametric inference under two-phase sampling; (ii) characterize a multiple robustness of the estimators that arises from the second-order remainder of a linear expansion of the target parameter; and (iii) provide the first comparison of the practical performance of these estimators in rigorous numerical experiments. Our contributions to the literature on stochastic interventions are (i) a novel estimator of a conditional density that is valid under two-phase sampling, while achieving a fast convergence rate, a crucial development for generating efficient estimators of the mean counterfactual; and (ii) an extension of nonparametric inference on mean counterfactuals under stochastic interventions using projections onto nonparametric working marginal structural models. Finally, we provide open source software packages, txshift  and haldensify , for the R statistical programming environment (R Core Team, 2020), that facilitate implementation of our proposed estimators and their nuisance regression functions, respectively.
The remainder of this work is organized as follows. In Section 2, we introduce the target parameter (2.1) and discuss prior work on stochastic interventions (2.2) and two-phase sampling (2.3). Section 3 includes discussion of strategies for estimating nuisance regression functions (3.1) and the formulation of two estimators (3.2) that achieve the nonparametric efficiency bound for estimation of the counterfactual mean under a stochastic intervention. Numerical experiments comparing the proposed estimators to simpler variants are presented in section 4. Section 5 describes an analysis of data from the HVTN 505 trial using our proposed methodology. We conclude by summarizing the results of our statistical and scientific investigations and discussing avenues for future investigation.
2 Preliminaries and Background 2.1 Notation, data, and target parameter Consider data generated by typical cohort sampling: the data on a single observational unit is denoted X = (W, A, Y ), where W ∈ W is a vector of baseline covariates, A ∈ A a real-valued exposure, and Y ∈ Y an outcome of interest. Initially, we assume access to n independent copies of X, using P X 0 to denote the distribution of X. We assume a nonparametric statistical model M X for P X 0 . We denote by q 0,Y the conditional density of Y given (A, W ) with respect to some dominating measure, q 0,A the conditional density of A given W with respect to dominating measure µ, and q 0,W the density of W with respect to dominating measure ν. We use p X 0 to denote the density of X with respect to the product measure. This density evaluated on a typical observation x may be expressed To define a counterfactual quantity of interest, we introduce a nonparametric structural equation model (NPSEM) to describe the data-generating process of X (Pearl, 2000). Specifically, we assume the data are generated by the following system of structural equations: where f W , f A , and f Y are deterministic functions, and U W , U A , and U Y are exogenous random variables.
The NPSEM provides a parameterization of p X 0 in terms of the distribution of the random variables (X, U ) modeled by the system of structural equations. Importantly, it also implies a model for the distribution of counterfactual random variables generated by specific interventions on the data-generating process. For example, a static intervention replaces f A with a value a in the support of A. By contrast, a stochastic intervention modifies the value A would naturally assume, f A (W, U A ), replacing it with a draw from a post-intervention distributionq 0,A (· | W ), where the zero subscript is included to emphasize that this distribution may depend on P X 0 . A static intervention may be viewed as a particular type of stochastic intervention in whichq 0,A (· | W ) places all mass on a single point. Díaz and van der Laan (2012) described a stochastic intervention that draws A from a distribution such thatq 0,A (a | W ) = q 0,A (a − δ(W ) | W ) for a user-supplied shifting function δ(W ) and a given a ∈ A. Haneuse and Rotnitzky (2013) showed that estimating the effect of this intervention is equivalent with that of an intervention that modifies the value A would naturally assume according to a regime d(A, W ). Importantly, the regime d(A, W ) may depend on both the covariates W and the exposure level A that would be assigned in the absence of the regime; consequently, this has been termed a modified treatment policy (MTP). Both Haneuse and Rotnitzky (2013) and Díaz and van der Laan (2018) considered an MTP of the form d(a, w) = a + δ(w) for δ(w) = γ ∈ R if a + γ ≤ u(w) and d(a, w) = a if a + γ > u(w), where u(w) is the maximum value in the support of q 0,A (· | W = w). This intervention generates a counterfactual random variable Y A+δ(W ) := f Y (A + δ(W ), W, U Y ) whose distribution we denote P δ 0 ; we seek to estimate ψ 0,δ := E P δ 0 {Y A+δ(W ) }, the mean of this counterfactual outcome. In the context of the HVTN 505 trial, this parameter corresponds to the counterfactual risk of HIV-1 infection had the levels of immune response markers of vaccinated participants been increased by γ units relative to the level induced by the current vaccine. This quantity may reflect the immune response marker distribution of a next-generation HIV vaccine with improved immunogenicity relative to the vaccine evaluated in HVTN 505. While the magnitude of shifting could generally be allowed to vary with participant characteristics, in our analysis, we consider an intervention that uniformly shifts all participants' immune responses by γ units, that is, d(a, w) = a + γ for all a ∈ A. Note that for HVTN 505, the parameter of interest is defined only for the vaccine group, making A a post-vaccination marker measuring an HIV-specific immune response. Importantly, it is not conceivable to define the target parameter for placebo recipients since only HIV-negative participants are enrolled in the trial and A is only defined if measured prior to HIV infection; consequently, all relevant placebo recipients have value zero for the marker A, and it is not meaningful to apply d(a, w) to shift the distribution of A.
Analysis of the HVTN 505 trial is complicated by its two-phase sampling design, a technique commonly used for assessing vaccine-mediated immune response in efficacy trials (Haynes et al., 2012). In general two-phase sampling, we do not observe X on all participants. Instead, we observe O = (W, C, CA, Y ) ∼ P 0 , where P 0 is the distribution of O and C is an indicator that an observation is included in the second-phase sample. In particular, C i = 1 if A is measured on the i th observation and C i = 0 otherwise. By convention, CA denotes that unobserved values of A are set to zero; this arbitrary labeling has no effect on our subsequent developments. In the context of vaccine trials, the probability of inclusion in the second-phase sample often depends on W and Y ; for each w and y we define g 0,C (y, w) := P(C = 1 | Y = y, W = w). Consequently, P 0 ∈ M = {P P X 0 ,g 0,C : P X 0 ∈ M X , g 0,C } -that is, P 0 is the distribution of O implied by the pair {P X 0 , g 0,C }. For example, in HVTN 505 all infected participants with samples available for marker measurement at week 28 had immune responses measured, i.e., g 0,C (1, w) = 1 for all w; however, only a subset of non-infected participants had immune responses measured. We will assume access to an i.i.d. sample O 1 , . . . , O n , denoting its empirical distribution by P n . The statistical estimation problem is thus to develop efficient nonparametric estimators of ψ 0,δ based on n independent realizations of the observed data unit O.

Identifying the counterfactual mean under a stochastic intervention
Díaz and van der Laan (2012) established that ψ 0,δ is identified by , the conditional mean of Y given A = a and W = w, as implied by P X 0 . For the statistical functional given in equation (1) to correspond to the causal estimand of interest, several untestable assumptions are required does not depend on a j + δ(w j ) for i = 1, . . . , n and i = j (Rubin, 1978(Rubin, , 1980. • Strong ignorability (no unmeasured confounding): A i ⊥ ⊥ Y ai+δ(wi) i | W = w i , for i = 1, . . . , n.
• Positivity: a i ∈ A =⇒ a i + δ(w i ) ∈ A | W = w i for all w ∈ W and i = 1, . . . n.
The positivity assumption required to establish equation (1) is unlike that required for static or dynamic interventions. In particular, it does not require that the post-intervention exposure density place mass across all strata defined by W .
Instead, for Q 0,Y to be well-defined, we require that the post-intervention exposure mechanism be bounded away from zero, i.e., P P X 0 {q 0,A (A − δ(W ) | W )/q 0,A (A | W ) > 0} = 1, which is satisfied by our definition of δ(W ). Díaz and van der Laan (2012) further provided the efficient influence function (EIF) of ψ 0,δ with respect to a nonparametric model, using this object as a key ingredient in the construction of their proposed estimators. The EIF, evaluated on a typical full data observation x, is where H(a, w) = 1(a < u(w)) (3)

Correcting for two-phase sampling
Since the earliest discussion of two-phase sampling (Neyman, 1938), a rich literature on such designs has emerged (e.g., Manski and Lerman, 1977;White, 1982;Flanders and Greenland, 1991;Pickles et al., 1995;Chatterjee et al., 2003;Breslow et al., 2009a,b;Dai et al., 2009;Gilbert et al., 2014;Kundu and Chatterjee, 2019). Early proposals of estimation strategies focused on parametric models of the sampling mechanism (Breslow and Cain, 1988), while Robins et al. (1994Robins et al. ( , 1995 introduced the first semiparametric estimator to incorporate inverse probability of sampling weights. Lawless et al. (1999) proposed a discretization-based procedure that uses stratum membership (across a small number of discrete strata) to construct the second-phase sample. Breslow and Holubkov (1997) proposed the use of maximum likelihood estimation (MLE) when a logistic regression may be used to link the outcome to all covariates, while Breslow et al. (2003) established the nonparametric MLE and its asymptotic properties for this class of problem. Subsequent proposals include pseudoscore estimators (Chatterjee et al., 2003), which may accommodate the use of continuous covariates in the sampling mechanism via kernel smoothing approaches but then require discreteness of covariates in the second-phase sample (Chatterjee and Chen, 2007); re-calibration under semiparametric regression (Chen and Chen, 2000;Fong and Gilbert, 2015); and targeted minimum loss estimation (Rose and van der Laan, 2011).
We build on the results of Rose and van der Laan (2011), who provide a study of nonparametric efficiency theory in two-phase sampling designs. In particular, they provide a representation of the EIF of a target parameter of the full data distribution P X 0 when the observed data O 1 , . . . , O n are generated by a two-phase sampling design. Their general developments suggest that, in the present problem, the EIF evaluated on a typical observation from the observed data o = (w, c, ca, y), is where given Y = y and W = w of the full data EIF amongst observations in the second-phase sample. Rose and van der Laan (2011) proposed two estimation strategies. The first -which we call the reweighted estimator -incorporates inverse probability weights (IPW) based on known or estimated values of the second-phase sampling probability g 0,C to an EIF-based estimation procedure such as one-step estimation or targeted minimum loss (TML) estimation. The estimator is shown to be asymptotically linear and efficient in designs where the sampling mechanism is known or can be estimated using nonparametric maximum likelihood. The second estimator described in this proposal is for situations where the sampling design is unknown and must be estimated using nonparametric regression, such as when the censoring process is outside the control of investigators; however, the authors did not provide a formal study of the theoretical properties of this estimator nor any numerical evaluations. Owing to its complexity, examples of this approach are extremely limited in the literature (Brown, 2014). We aim to fill in these gaps by providing a formal theorem establishing conditions under which this estimator will achieve asymptotic efficiency, as well as a simulation study demonstrating its practical performance in the context of estimating the counterfactual mean under a stochastic intervention.

Methodology
We now describe the implementation of our proposed estimators. We utilize two frameworks for estimation: the one-step framework (Pfanzagl and Wefelmeyer, 1985;Bickel et al., 1993) and TML estimation (van der Laan and Rubin, 2006;Rose, 2011, 2018). Both estimators develop in two stages, with crucial differences appearing in the second stage. In the first stage, we construct initial estimators of key nuisance quantities, while in the second stage we perform a bias-correction based on the estimated EIF. The one-step framework updates an initial substitution estimator by adding the empirical mean of the estimated EIF. By contrast, the TML estimation framework uses a univariate logistic tilting model to build a targeted estimator of Q 0,Y , which is subsequently used to construct an updated substitution estimator.

Estimating nuisance parameters
Our general strategy for estimating nuisance parameters relies on first using the entire observed data set to estimate the second-phase sampling probabilities, g 0,C . Subsequently, inverse probability of sampling weights based on these estimates are used to generate estimates of relevant full data quantities using data available only on observations in the second-phase sample. These quantities include the outcome regression Q 0,Y , the exposure density q 0,A , and the joint distribution of covariates and exposure, which we denote by Q 0,AW . Finally, estimates of full data quantities are used to estimate G 0 , the conditional mean of the full data EIF given Y and W amongst observations included in the second-phase sample.
Excepting Q 0,AW , which we estimate using an inverse probability of sampling weighted empirical distribution, we describe both parametric and flexible, data adaptive estimators. The data adaptive estimators are more parsimonious with our theoretical developments, which pertain to nonparametric-efficient estimation; nevertheless, our developments hold equally well for parametric working models. In Theorem 1, we detail assumptions on the stochastic behavior of estimators of these nuisance functions and relate these to the behavior of the resultant estimator of the target parameter.
An estimator of the sampling mechanism g 0,C could be derived from any classification method (e.g., logistic regression), in which P P0 (C = 1 | Y, W ) is estimated using the full sample; however, nonparametric or semiparametric estimation may be preferable depending on the availability of information about the two-phase sampling design.
To generate an estimate Q n,AW of the full data joint distribution of (A, W ), we use a stabilized inverse probability weighted empirical distribution. For a given (a, w), .
To estimate Q 0,Y , one may again use any classification or regression model, where Y is the outcome and functions of A and W are included as predictors. In fitting this model, inverse probability of sampling weights C i /g n,C (Y i , W i ) for i = 1, . . . , n, are included to account for the two-phase sampling design. Any valid regression estimator may be leveraged for this purpose, so long as the implementation of the estimator respects the inclusion of sample-level weights; in practice, we recommend the use of a semiparametric or nonparametric estimator. We denote by Q n,Y (a, w) the estimate evaluated on a data unit with A = a, W = w.
The simplest strategy for estimating the generalized propensity score q 0,A is to assume a parametric working model and use standard regression techniques to generate suitable estimates of the density. For example, one could operate under the working assumption that A given W follows a Gaussian distribution with homoscedastic variance and mean p j=1 β j φ j (W ), where φ = (φ j : j) are user-selected basis functions and β = (β j : j) are unknown regression parameters. In this case, a density estimate would be generated by fitting a linear regression of A on φ(W ) (e.g., minimizing inverse probability of sampling weighted least squares) to estimate the conditional mean of A given W , paired with an estimate of the variance of A. In this case, the estimated conditional density is given by the density of a Gaussian distribution evaluated at these estimates.
The relative dearth of available estimators of a conditional density motivated our development of a novel estimator that accounts for two-phase sampling designs. We detail this approach in the Supplementary Materials and provide an implementation of our proposal in the haldensify R package . Going forward, we denote by q n,A (a | w) the estimated conditional density of A given W = w, evaluated at a ∈ A.
The final nuisance parameter that must be estimated is G 0 , the conditional mean of the random variable D F (P X 0 )(O) given (Y, W ) amongst those included in the second-phase sample. To estimate this quantity, we generate a pseudooutcome as follows. First, define the substitution estimator, as well as the auxiliary term Using these quantities, for all i such that C i = 1, we compute A simple estimation strategy for G 0 is to adopt a parametric working model and fit, for example, a linear regression of the pseudo-outcome D F n,i on basis functions of Y and W . Importantly, since G 0 is defined as a conditional expectation with respect to the observed data distribution, we need not include inverse probability of sampling weights in this regression estimate. While a parametric working model for G 0 is permissible, given the complexity of the object, correct specification of this model is likely challenging and we recommend more flexible approaches. We let G n (Y i , W i ) denote the value of the chosen regression estimator evaluated on the i th observation i = 1, . . . , n.

One-step estimator
Based on the estimated nuisance functions detailed above, efficient estimators may be constructed using either of the one-step or targeted minimum loss estimation frameworks. The one-step estimator adds the empirical mean of the estimated EIF to the initial plug-in estimator, The resultant augmented one-step estimator ψ + n,δ relies on the nuisance functions estimators (Q n,Y , g n,A , G n , g n,C ). Theorem 1 details sufficient assumptions on these estimators for ensuring that the one-step is asymptotically efficient.

Targeted minimum loss estimator
An asymptotically linear TML estimator of ψ 0,δ may be constructed by using inverse probability of sampling weights to update the initial estimator of the outcome mechanism Q n,Y to an estimator Q n,Y . A plug-in estimator is then constructed, ψ n,δ := Q n,Y (a + δ(w), w)dQ n,AW (a, w) .
This updated estimator Q n,Y is constructed as follows.
1. Define a working logistic regression model for the conditional mean of C given {Y, W }, logit(g n,C,ξ ) = logit(g n,C ) + ξ(G n /g n,C ) : ξ ∈ R. An estimate ξ n of the parameter ξ of this model is computed using maximum likelihood, and we define the updated estimator g n,C := g n,C,ξn .
2. Next, define a working logistic regression model for the conditional mean of Y given {A, W }: logit(Q n,Y, ) = logit(Q n,Y ) + H n : ∈ R. An estimate n of is obtained using weighted logistic regression with weights The outlined procedure includes two targeting steps. The first of these steps constructs an update of the initial estimator of the second-phase sampling probability, g n,C , based on the initial estimate of G n . This step ensures that the revised estimate satisfies in a single step when a universal least favorable submodel (van der Laan and Gruber, 2016) is used, though an iterative procedure may be used to achieve the same result. In the second step, the updated outcome regression Q n,Y is generated based on the conditional density estimate q n,A ; the inclusion of weights in the regression ensures that When the first step of this procedure is omitted, the resultant TML estimator is equivalent to the reweighted estimator of Rose and van der Laan (2011). The additional step allows our estimator to attain asymptotic linearity in a broader set of circumstances. In particular, the reweighted estimator requires that the sampling weights are known or can be estimated at a parametric rate. Our approach allows more flexible estimators of sampling weights to be used.

Asymptotic analysis of efficient estimators
We establish the asymptotic efficiency of our estimators in the theorem below. The theorem depends on a several regularity conditions, which are discussed in the Supplementary Materials. The theorem is provided in the context of the TML estimator, but, with a similar set of assumptions, the same result holds for the one-step estimator; for brevity, we omit this analogous theorem. In the sequel, D F (Q 0,Y , q 0,A ) and D F (P X 0 ) are used interchangeably since D F depends on P X 0 primarily through estimates of the nuisance functions Q 0,Y and q 0,A . Theorem 1 (Asymptotic linearity and efficiency of the TML estimator ψ n,δ ). Assuming conditions C1-C6, An immediate corollary of Theorem 1 is that ψ n,δ is asymptotically efficient, since it is an asymptotically linear estimator with influence function equal to the efficient influence function. Moreover, the central limit theorem implies that the scaled, centered estimator converges in distribution to a mean-zero Gaussian random variable with variance The proof of Theorem 1 is given in the Supplementary Materials. The conditions of the theorem are standard in semiparametric inference problems, essentially requiring a sub-parametric rate of convergence of each of the nuisance estimators to their true counterparts, a Donsker class condition on the EIF evaluated at the estimated nuisance parameters, and L 2 (P 0 )-consistency of this same object.
With respect to this first condition, we note that the highly adaptive lasso (HAL) regression estimator has been shown to achieve a sufficiently fast rate of convergence so as to satisfy the requirements of the theorem (van der Laan, 2017; Benkeser and van der Laan, 2016; Bibaut and van der Laan, 2019; Coyle et al., 2019) under the assumption that the true regression function is right-hand continuous with left-hand limits and bounded sectional variation norm. This provided further motivation for our development of a HAL-based conditional density estimator. We note that our simulation studies and analysis of the HVTN 505 trial data utilize HAL to increase the applicability of our theorem.
With respect to the Donsker condition, we note that this assumption may be avoided by making use of cross-validation (also known as cross-fitting) in the process of estimating nuisance parameters (Klaassen, 1987;Zheng and van der Laan, 2011;Chernozhukov et al., 2016). Such an estimator enjoys the same asymptotic properties as our proposed non-sample-splitting counterpart while eschewing reliance on the Donsker condition.

Multiple robustness of efficient estimators
The EIF of our estimators enjoys a multiple robustness property, which allows our estimators to achieve consistency even in situations where certain combinations of nuisance parameters are inconsistently estimated.
Lemma 1 (Multiple robustness of the EIF). Let (G, g C , Q Y , q A ) denote the limits of the nuisance estimators (G n , g n,C , Q n,Y , q n,A ) in probability. Suppose either of the following two conditions hold In the case of the one-step estimator, the initial nuisance function estimates g n,C and Q n,Y are used instead. The lemma implies that our efficient estimators will be asymptotically consistent if at least one of (G 0 , g 0,C ) and one of (Q 0,Y , q 0,A ) are consistently estimated.

Confidence intervals and hypothesis tests
Theorem 1 established the limiting distribution of our efficient estimators. From the limit distribution, inference for either estimator may be attained in the form of Wald-type confidence intervals and corresponding hypothesis tests.
These procedures are asymptotically justified under the conditions of theorem 1. Importantly, while multiple robustness implies that consistent estimation of ψ 0,d is possible under inconsistent estimation of some nuisance parameters, the validity of confidence intervals and hypothesis tests requires consistent estimation of all nuisance parameters.

Summarization via working marginal structural models
Estimation of ψ 0,δ for a single pre-specified shift δ may be unsatisfactory in some contexts, as it does not provide information concerning a dose-response relationship between exposure and outcome. Thus, to develop an understanding of a dose-response pattern in the context of stochastic interventions, it may be informative to estimate the counterfactual mean outcome across several values of δ. In the context of HVTN 505, we consider estimation of a grid of counterfactual means ψ 0 = (ψ n,δ1 , . . . , ψ n,δ K ) and examine how the risk of HIV infection varies with choice of δ over a fixed grid, i.e., δ k ∈ {δ 1 , . . . , δ K }. After estimating the counterfactual mean for each δ k in the grid, a summary measure of the effect of the stochastic interventions on the mean counterfactual outcomes may be constructed by projection onto a working marginal structural model (MSM). For example, we might consider a (possibly weighted) least-squares projection on the linear working model m β (δ) = β 0 + β 1 δ. In this simple case, the parameter β 1 corresponds to the linear trend in mean counterfactual outcomes as a function of the shift values in the grid.
More generally, we can define β(δ) = argmin β∈R d δ∈{δ1,...,δ K } h(δ){ψ 0,δ − m β (δ)} 2 , for a user-selected weight function h(δ). We note that adjustment of the weight function, as well as the functional form of m β (δ), allow for a wide variety of working models to be considered. Alternatively, β(δ) can be viewed as the solution of The goal is to make statistical inference on the parameter β. We note that this approach does not assume a linear dose-response curve, but rather uses a working model to summarize the relationship between exposure and outcome (Neugebauer and van der Laan, 2007). This approach is distinct from that of Haneuse and Rotnitzky (2013), whose proposal involving MSMs pertains specifically to parametric models.
To estimate β, we assume access to the TML or one-step estimates ψ n = (ψ n,δ1 , . . . , ψ n,δ K ) for each δ k in the prespecified grid. The estimate β n of β is the solution in β of the equation 0 = U (β, ψ n ). To derive the limit distribution of β n , let D 0,ψ denote a vector whose k th entry is the EIF associated with parameter ψ 0,δ k . The delta method implies that the influence function of β n is D β = [− d dβ U (β, ψ 0 ) −1 ] d dψ0 U (β, ψ 0 )D ψ0 , and that n 1/2 (β n − β) converges in distribution to a mean-zero Gaussian random variable with variance Σ = E P0 {D β (O) 2 }. The empirical covariance matrix of D 0,ψ evaluated at observed nuisance parameter estimates may serve as an estimate of Σ.

Simulation Studies
The proposed estimators were evaluated using two simulation experiments. In the first, we compare our proposed estimators to alternative estimators proposed in the literature. To highlight the benefits offered by our approach over the simple reweighted estimator of Rose and van der Laan (2011), we focus on how estimation of g 0,C influences the estimator performance. In particular, we consider both the usage of standard logistic regression and the highly adaptive lasso in the construction of g n,C .
In a second simulation study, we analyze the performance of our estimators using a data-generating mechanism that is inspired directly by data from the HVTN 505 trial. Here, we focus on comparing the relative performance of our efficient, augmented one-step and TML estimators to assess their capabilities for use in real-world data analysis.
In both sets of simulation studies, we consider estimation of ψ 0,δ for several values of δ. The simulation studies were performed using the txshift software package  for the R language and environment for statistical computing (R Core Team, 2020), available at https://github.com/nhejazi/txshift.

Simulation #1: Comparison of estimator variants under differences in sampling mechanism estimation
In this setting, we evaluate the estimators on data simulated from the following data-generating mechanism: where expit(x) = {1+exp(−x)} −1 . In this setting, the outcome has a fairly high event rate, with P(Y = 1 | A, W ) ≈ 0.415. We used this data-generating process to sample n i.i.d. observations for n ∈ {100, 400, 900, 1600, 2500} and used the resultant data to estimate the target parameter with each of the estimators considered. This was repeated 1000 times. We considered estimation of ψ 0,δ for δ ∈ {−0.5, 0, 0.5}, where the corresponding true values of the target parameter were approximately {0.501, 0.415, 0.333}.
We compared the reweighted estimators of Rose and van der Laan (2011) to our proposed estimators. For reference, we also present the results of a naive estimate that ignores the two-phase sampling design. In each of these three cases, we consider one-step and TML estimators, giving six estimators in total. Each of these six estimators was constructed by estimating the exposure mechanism q n,A and outcome mechanism Q n,Y via maximum likelihood based on correctly specified parametric models, while g n,C was constructed using either logistic regression or the highly adaptive lasso. Although a parametric regression model adequately captures the true functional form of the sampling mechanism in this data-generating process, one may not have access to such information about the sampling mechanism in practice. In this case, nonparametric regression is preferable. This motivated our desire to examine the relative performance of the reweighted versus our proposed estimators under such a setting. Based on theory, when g n,C is estimated using logistic regression, we would expect both the reweighted estimator and our proposed estimator to be asymptotically linear, which would be supported by observing that the bias of the estimators disappears faster than n −1/2 -rate. On the other hand, when g n,C is nonparametrically estimated (in our simulation, via the highly adaptive lasso), the reweighted estimators should not achieve asymptotic linearity, while our proposed estimators should. The naive estimators, which make no adjustment for the two-phase sampling design, were expected to perform poorly -even though this is a "best case" scenario for these estimators in the sense that the two nuisance parameters are correctly estimated using parametric models.
We compared all estimators in terms of their bias (scaled by n 1/2 ), mean squared-error (scaled by n), and coverage of 95% Wald-style confidence intervals.  When the sampling mechanism is estimated via a correctly specified parametric model (upper panel), the reweighted and our proposed estimators behave as expected, with low bias and stable MSE. However, the reweighted estimators display coverage exceeding 95%, while our proposed estimators achieve nominal coverage. This occurs because the influence function provided in Rose and van der Laan (2011), which is the basis for the standard error estimates used to build these confidence intervals, does not include the first-order contribution due to estimation of the sampling mechanism. The result is a conservative standard error estimate. Unsurprisingly, we found that the naive estimator performed poorly in all sample sizes, highlighting the importance of appropriately accounting for sampling design.
When the sampling mechanism was estimated using HAL (lower panel), the reweighted estimators do not attain asymptotic linearity, as evidenced by the scaled bias and MSE increasing with sample size. On the other hand our proposed estimators have small bias and MSE approaching the efficiency bound, thus demonstrating the benefits of the additional effort required to produce our estimators over the simpler reweighted estimators.
For consistency with the observed rate of HIV infection in the HVTN 505 trial, the outcome was made to have a relatively rare event rate, with P(Y = 1 | A, W ) ≈ 0.059. We considered a setting where we observed n = 1400 i.i.d. observations, approximately matching the sample size of the vaccinated arm in HVTN 505. Estimator performance was assessed and reported by aggregating across 1000 repetitions. Under this data-generating mechanism, the true parameter values were approximated as ψ 0,δ = {0.0627, 0.0617, 0.0609, 0.0598, 0.0589, 0.0580, 0.0571, 0.0561, 0.0554} for δ ∈ {−2.0, −1.5, −1.0, −0.5, 0.0, 0.5, 1.0, 1.5, 2.0}, respectively. We also considered an analogous setting, but where the effect of the treatment on the outcome was removed. The results were largely similar and are presented in the Supplementary Materials.
The proposed estimators were constructed by using the highly adaptive lasso to estimate the sampling mechanism g n,C , the exposure mechanism q n,A , the outcome mechanism Q n,Y , and the pseudo-outcome regression G n . We compared the proposed one-step and TML estimators in terms of their bias and MSE. The results are summarized in Figure 2. Inspection of the bias and MSE indicates that both of our proposed estimators display adequate performance, with small finite-sample bias and mean-squared error. In terms of bias and MSE, the performance of the two estimators is nearly indistinguishable. From these numerical investigations, we concluded that our augmented estimators are both well-suited to estimating ψ 0,δ in the context of a re-analysis of data from the HVTN 505 trial.

Application to the HVTN 50Trial
The HVTN 505 trial was a randomized control trial that enrolled 2504 HIV-negative participants and randomized participants 1-to-1 to receive an active vaccine or placebo. The one-year incidence of HIV-1 infection was about 1.8% per person-year in the vaccine arm and 1.4% per person-year in the placebo arm, during primary follow-up for HIV-1 acquisition (between week 28 and month 24; the same period as was used for assessment of immune correlates). Blood was drawn post-vaccination and immune responses measured via intracellular cytokine staining of preserved HIV-1-stimulated peripheral blood mononuclear cells for all HIV-1 cases diagnosed between week 28 and month 24 and a random sample of uninfected controls (Janes et al., 2017). The two-phase sampling of vaccine recipient controls without-replacement sampled five controls per case within each of eight baseline covariate strata defined by categories of body mass index and race/ethnicity (White, Hispanic, Black). Janes et al. (2017) provide the first immune correlates analysis of the HVTN 505 data based on the constructed two-phase sample, focusing on two cellular immune responses: the Envelope CD4+ and CD8+ polyfunctionality scores. Subsequently, Fong et al. (2018) analyzed an expanded set of immune markers, including six primary humoral immune response variables beyond the two first studied by Janes et al. (2017). Both sets of analyses found the CD4+ and CD8+ polyfunctionality scores to be informative of HIV-1 infection status by month 24 of the study (end of follow-up).
We examined how a range of posited shifts in standardized polyfunctionality scores of the CD4+ and CD8+ immune markers would impact the mean counterfactual risk of HIV-1 infection in vaccine recipients. We considered standardized polyfunctionality scores, so that our pre-specified grid of shifts δ ∈ {−2.0, −1.5, −1.0, −0.5, 0.0, 0.5, 1.0, 1.5, 2.0} can be interpreted as shifts on the scale of standard deviation (sd). We present results based on our TML estimator; results for the one-step estimator were similar. In order to summarize the manner in which the mean counterfactual risk of HIV-1 infection varies with shifts in the polyfunctionality scores, we consider projection onto a linear working MSM, as discussed in section 3.4. Our augmented TML estimator ψ n,δ for the mean counterfactual risk of HIV-1 infection requires the construction of initial estimators of the nuisance functions {g 0,C , G 0 , Q 0,Y , q 0,A }.
We used flexible estimation strategies for each of the nuisance parameters. The conditional probability of inclusion in the second-phase sample was estimated using HAL, adjusting for age, sex, race/ethnicity, body mass index, and a behavioral risk score for HIV-1 infection. The density q n,A of the CD4+ or CD8+ polyfunctionality scores, conditional on the same set of covariates, was estimated using our proposed HAL-based density estimator . The outcome regression Q 0,Y was estimated using super learner (van der Laan et al., 2007), which was used to build an ensemble model as a weighted combination of prediction functions from a diverse collection of candidate algorithms. Further details, including choices of the learning algorithms making up the library and the weights assigned to each by the super learner ensemble, are given in the Supplementary Materials. Finally, the pseudo-outcome regression, G n , was performed via HAL.
We note that our choice of estimation strategy -specifically, the use of HAL regression for g 0,C and HAL-based conditional density estimation for q 0,A -ensures a sufficiently fast rate of convergence of these nuisance parameter estimates to their true values under only mild assumptions. Importantly, this allows for the second case described in lemma 1 to be fulfilled, ensuring consistency of the resultant TML estimates. This consistency of the estimates would hold even in the case that the HAL regression and super learner ensembles used for G 0 and Q 0,Y , respectively, proved inconsistent; however, as both of these nuisance parameters are also estimated via flexible modeling approaches, it would be expected that the conditions of theorem 1 are made to hold. Thus, we expect that our TML estimates of the counterfactual risk of HIV-1 infection across the grid in δ will be consistent and (quite possibly) nonparametricefficient. Taking advantage of the robustness built into our estimation strategy lends an additional degree of reliability as to the results of our analysis of the HVTN 505 data. Finally, we recall that the estimate of β 0 , the slope of the marginal structural model through the TML estimates of the counterfactual risk in δ, is endowed with the same consistency and efficiency properties as the individual TML estimates.
Results of applying our estimation procedure separately to both the CD4+ and CD8+ polyfunctionality scores are presented in Figure 3.  Examination of the point estimates and confidence intervals of ψ 0,δ in Figure 3 reveals that downshifts in the CD4+ polyfunctionality score led to a small increase in estimated HIV-1 infection risk (Figure 3, top panel). For example, a shift of two standard units lower in the CD4+ polyfunctionality score was found to double the risk of HIV-1 infection. The estimated slope parameter of the working MSMβ TMLE pointed to an estimated decrease in risk of about -0.3% per standard unit of CD4+ polyfunctionality change.
The estimated result of shifts in the polyfunctionality score of the CD8+ immunogenic marker displayed a markedly much stronger relationship with the risk of HIV-1 infection (Figure 3, lower panel). While positive shifts of the standardized CD8+ polyfunctionality score beyond those observed in the trial do not appear to have a strong effect on the risk of HIV-1 infection, shifts that lower the CD8+ polyfunctionality score display a strongly negative linear trend, indicating that decreases in the activity of the CD8+ marker would adversely affect the risk of HIV-1 infection. At the largest negative shift considered, the counterfactual risk of HIV-1 infection is over four times that of the risk observed in the HVTN 505 trial.
Overall, the results of our analyses support the conclusions of Janes et al. (2017) and Fong et al. (2018), further indicating that modulation of the CD4+ and CD8+ polyfunctionality scores may reduce the risk of HIV-1 infection, with CD8+ polyfunctionality playing a particularly important role. Notably, our analysis differs from the previous efforts in two ways: our estimates (i) are based on a formal causal model, which provides an alternative estimand to summarize relationships between immunogenic response and risk of HIV-1 infection, and (ii) leverage machine learning to allow the use of flexible modeling strategies while simultaneously delivering robust inference.

Discussion
Our analysis of the HVTN 505 trial could be improved in several respects. First, there was participant dropout observed in the trial, which our analysis ignored. A more robust analysis could leverage available covariate information to account for potentially informative missingness. Beyond this issue, there are several other directions for potentially interesting extensions. It would be of interest to extend our estimation strategy to other effects based on stochastic interventions, including the population intervention (in)direct effects (Díaz and Hejazi, 2020) defined in mediation analysis and the incremental propensity score interventions of Kennedy (2019). Extending our estimation strategy to such settings, and its application in the analysis of other vaccine efficacy trials, will be the subject of future research.

Supplementary materials
Conditional density estimation based on the highly adaptive lasso In section 3.1, several of the challenges associated with conditional density estimation were briefly expressed. The auxiliary covariate given in equation (3) of the EIF is a ratio of such conditional densities. In the construction of our proposed estimators, the exposure mechanism q 0,A is a density of the intervention A, conditional on the observed covariates W , which must be evaluated at both the observed values and the post-intervention counterfactual values. As consistent estimation of the the generalized propensity score is an integral part of our proposed methodology, we outline here a conditional density estimator, built around the HAL regression function, that achieves the convergence rate required by our formal theorem. We note that proposals for the data adaptive estimation of such quantities are sparse in the literature (e.g., Zhu et al., 2015). Notably, Díaz and van der Laan (2011) gave a proposal for constructing semiparametric estimators of such a target quantity based on exploiting the relationship between the hazard and density functions. Our proposal builds upon theirs in several key ways: (i) we adjust their algorithm so as to incorporate sample-level weights, necessary for making use of inverse probability of censoring weights; and (ii) we replace their use of an arbitrary classification model with one based on the HAL regression function. While our first modification is general and may be applied to the estimation strategy Díaz and van der Laan (2011) propose, our latter contribution requires adjusting the penalization aspect of HAL regression to respect the use of a loss function appropriate for prediction on the hazard scale. To make these contributions widely accessible, we introduce the haldensify R package , available at https://github.com/nhejazi/haldensify To build an estimator of a conditional density, Díaz and van der Laan (2011) considered discretizing the observed a ∈ A based on a number of bins T and a binning procedure (e.g., including the same number of points in each bin or forcing bins to be of the same length). We note that the choice of the tuning parameter T corresponds roughly to the choice of bandwidth in classical kernel density estimation; this will be made clear upon further examination of the proposed algorithm. The data {A, W } are reformatted such that the hazard of an observed value a ∈ A falling in a given bin may be evaluated via standard classification techniques. In fact, this proposal may be viewed as a re-formulation of the classification problem into a corresponding set of hazard regressions: where the probability that a value of a ∈ A falls in a bin [α t−1 , α t ) may be directly estimated from a standard classification model. The likelihood of this model may be re-expressed in terms of the likelihood of a binary variable in a data set expressed through a repeated measures structure. Specifically, this re-formatting procedure is carried out by creating a data set in which any given observation A i appears (repeatedly) for as many intervals [α t−1 , α t ) that there are prior to the interval to which the observed a belongs. A new binary outcome variable, indicating A i ∈ [α t−1 , α t ), is recorded as part of this new data structure. With the re-formatted data, a pooled hazard regression, spanning the support of A is then executed. Finally, the conditional density estimator for α t−1 ≤ a < α t , may be constructed. As part of this procedure, the hazard estimates are mapped to density estimates through rescaling of the estimates by the bin size (α t − α t−1 ).
In its original proposal, a key element of this procedure was the use of any arbitrary classification procedure for estimating P(a ∈ [α t−1 , α t ) | W ), facilitating the incorporation of flexible, data adaptive estimators. We alter this proposal in two ways, (i) replacing the arbitrary estimator of P(a ∈ [α t−1 , α t ) | W ) with HAL regression and (ii) accommodating the use of sample-level weights, making it possible for the resultant conditional density estimator to achieve a convergence rate with respect to a loss-based dissimilarity of n −1/4 under only mild assumptions. This is an important advance that is needed for the asymptotic analysis of our proposed estimators, as per section 3.2.3. As a secondary advance, our procedure alters the HAL regression function to use a loss function tailored for estimation of the hazard, invoking 1 -penalization in a manner consistent with this loss.

Proof of Theorem 1
We now examine a proof of the theorem establishing conditions for the weak convergence of our efficient estimators. Building upon section 2, note that the full data parameter may be expressed as a mapping Ψ F : M X → R and that Ψ F (P X 0 ) ≡ Ψ F (Q 0,Y ), since the parameter mapping depends on P X 0 only through the functional Q 0,Y . We recall that the EIF D F (Q 0,Y , g 0,A ) coincides with the canonical gradient of the parameter mapping Ψ F : M X → R, since a regular asymptotically linear estimator with influence function equal to the canonical gradient is asymptotically efficient (Bickel et al., 1993;van der Laan and Robins, 2003). Note further that the observed data parameter is defined such that Ψ(P 0 ) ≡ Ψ F (P X 0 ), where the observed data parameter mapping Ψ : M → R is pathwise differentiable at a distribution P in the statistical model with a gradient given by the EIF: The class of all gradients of Ψ at P is given by {D(G 0 , g 0,C , D F (P X )) : D F (P X )} where D F (P X ) varies over all gradients of the full data parameter Ψ F at distributions P X ∈ M X . As D F varies, G 0 also necessarily varies since it is defined as the conditional mean of D F given {C = 1, Y, W }. In particular, if the full data model M X is nonparametric, then there is only one full data gradient, which is the canonical gradient (or EIF) of Ψ at P .
In the sequel, let f 2 = E{f (O) 2 } 1/2 denote the L 2 (P 0 ) norm of a P 0 -measurable function f , define G 0 := , and let G n be the estimate of G 0 . An exact characterization of the second-order remainder for the full data parameter Ψ F (P X ) − Ψ F (P X 0 ) is given by ] while the exact second-order remainder for the observed data parameter Ψ(P ) − Ψ(P 0 ) is analogously defined Combining these definitions, the exact second-order remainder for the observed data parameter may be expressed in terms of the second-order remainder of its full data counterpart: Assume the following conditions: C1. E Pn D(G n , g n,C , D F (Q n,Y , q n,A )) = 0.
C6. Let F v (M ) be the class of cadlag functions f on a cube [0, τ ] ⊂ R d (for some integer d), for which the sectional variation norm f v is bounded by a universal constant M < ∞. Assume that D(G n , g n,C , D F (Q n,Y , q n,A )) ∈ F v (M ) with probability tending to 1 (n.b., the definition F v (M ) can be replaced by any Donsker class).
Proof [Theorem. 1: asymptotic linearity and efficiency of the TML estimator ψ n,δ ] Under conditions C1-C6, we have E Pn D(G n , g n,C , D F (Q n,Y , q n,A )) = o P (n −1/2 ); moreover, by definition of R 2 (P, P 0 ): Combining with condition C1, we have Conditions C5 and C6 imply that the sum of the first two terms equals E Pn [D(G 0 , g 0,C , D F (Q 0,Y , q 0,A )(O)] + o P (n −1/2 ). Condition C4 ensures the third term equals o P (n −1/2 ). Conditions C2 and C3 imply that the fourth term equals o P (n −1/2 ), which completes the proof. An analogous proof holds for the one-step estimator ψ + n,δ when the initial estimates Q n,Y and g n,C are used instead of their revised counterparts. Here, we do not require condition C1, as our proposed estimation procedure guarantees that it will be attained.
Simulation #1: Comparison of estimator variants under the shifts δ = −0.5 and δ = 0 In the results reported in section 4.1, for our first simulation study under the shift δ = 0.5, we noted excellent performance of our proposed estimator variants under several standard metrics, including √ n-bias, n-MSE, and the coverage of confidence intervals. To further assess the quality of performance of our augmented estimators, we examine the same six estimator variants under the shifts δ ∈ {−0.5, 0}. Details of the simulation study have been previously described in section 4.1. As before, results are reported based on aggregation across 1000 repetitions. The results of these numerical investigations are reported in Figures 4 and 5. Both the reweighted estimators of Rose and van der Laan (2011) and our augmented estimators display the expected level of performance when the sampling mechanism is estimated via a correctly specified parametric model, standing out against the performance of their unadjusted counterparts. As in the case of δ = 0.5, the reweighted estimators display coverage exceeding the desired 95% level, while their augmented analogs cover at exactly the nominal level. This suggests again that the reweighted estimators exhibit an inflated variance relative to that of our augmented estimators. The lower panel of each figure visualizes differences in the performance of the estimators when the sampling mechanism is estimated flexibly via HAL. These results reveal a significant discrepancy in the performance of the reweighted and augmented estimators, with our augmented one-step and TML estimators outperforming their reweighted counterparts in terms of √ n-bias, n-MSE, and coverage. In this second case, we note that the TML estimators appear to exhibit a small degree of estimation instability at n = 100, displaying larger n-MSE in both the reweighted and augmented variants. We conjecture that this performance is due to an instability induced by the targeting procedure -in practice, this could be ameliorated by alterations to the convergence criterion. Overall, the results of these numerical investigations do not differ substantially from those presented in section 4.1, establishing that our augmented estimators improve upon alternative estimation strategies when nonparametric estimation of g 0,C is performed.
Simulation #2: Comparing proposed estimators in a scenario like the HVTN 505 trial (no effect of treatment) To further assess the performance of our augmented estimators in a scenario like the HVTN 505 trial, we replace the structural equation for the outcome Y with a draw from the Bernoulli distribution parameterized as Bernoulli(p = expit(−2.8 − 0.0013W 1 − 0.0016W 2 + 0.0678W 3 + 0.039W 4 )) to achieve a setting in which there is no effect of the treatment on the outcome. In this case, the outcome process is relatively rare, with P(Y = 1 | A, W ) ≈ 0.053. The structural equations for other components of the observed data O are unaltered. In this setting, the true parameter value was approximated to be ψ 0,δ = 0.0526 regardless of the value of δ. We present the results of assessing our augmented estimators in Figure 6. As the ground truth of the effect of shifting the post-vaccination activity of the CD4+ and CD8+ immunogenic markers on the risk of HIV-1 infection is unknown in the HVTN 505 trial, we assess our estimators in this scenario so as to be sure of there performance when no effect of treatment is present. Estimation of all nuisance parameters is performed using exactly the same techniques outlined previously in section 4.2.
In terms of bias and MSE, our proposed estimators display adequate performance: the augmented one-step and TML estimators are effectively indistinguishable in terms of MSE, while, in terms of bias, the augmented TML estimator appears to outperform the one-step estimator uniformly across all values of δ, likely on account of the enhanced finite-sample performance of TML estimators (van der Laan and Rose, 2011). As with the numerical investigations presented in section 4.2, these results suggest that our augmented estimators will perform well enough to allow accurate estimation of ψ 0,δ when applied to data from the HVTN 505 trial.
Ensemble models for Q n,Y in the HVTN 505 re-analysis As noted in section 5, for both the CD4+ and CD8+ analyses, the estimated outcome mechanism Q n,Y was constructed from an ensemble model based on the super learner algorithm (van der Laan et al., 2007). The cross-validation selector that forms the basis of the super learner has been shown to exhibit unique theoretical guarantees, including asymptotic equivalence to the oracle selector (van der Laan et al., 2004(van der Laan et al., , 2006van der Vaart et al., 2006), that make its use preferable over other ensemble learning approaches.
A rich library of candidate classification algorithms was considered in the super learner ensemble model for Q n,Y . These included 1 -penalized lasso regression (Tibshirani, 1996;Friedman et al., 2009); 2 -penalized ridge regression (Tikhonov and Arsenin, 1977;Hoerl and Kennard, 1970;Friedman et al., 2009); three elastic net regressions weighting the 1 penalty at α ∈ {0.25, 0.50, 0.75} and the 2 penaly at (1−α) (Zou and Hastie, 2003;Friedman et al., 2009); random forests (Breiman, 2001) with 50, 100, and 500 trees based on its implementation in the ranger R package (Wright et al., 2017); extreme gradient boosted trees with 20, 50, 100, and 300 fitting iterations (Chen and Guestrin, 2016); multivariate adaptive polynomial spline regression (Kooperberg et al., 1997;Stone et al., 1994); multivariate adaptive regression splines (Friedman et al., 1991); generalized linear models with Bayesian priors on parameters; a multilayer perceptron (Rosenblatt, 1961); and the highly adaptive lasso (van der Laan, 2017; Benkeser and van der Laan, 2016; Coyle et al., 2019). The implementation of the super learner algorithm in the sl3 R package (Coyle et al., 2020) was used, and weights assigned to each learning algorithm by the super learner are given in the Supplementary Materials in Tables 1 and 2 for the CD4+ and CD8+ analyses, respectively. In the super learner ensemble model for the outcome regression in the CD4+ analysis, the three best learning algorithms were a GLM with Bayesian priors on parameter estimates, a polynomial spline regression model, and an elastic net regression model that favored the 1 (lasso) penalty over the 2 (ridge) penalty. Another variant of elastic net regression, which favored the 2 penalty over the 1 penalty, and ridge regression were also given nontrivial weights in the ensemble model. A variant of extreme gradient boosting and the highly adaptive lasso received low weights, while all other candidate algorithms in the library were assigned weights of zero. In the CD8+ analysis, the three best learning algorithms, chosen by the super learner ensemble model for the outcome regression, were a GLM with Bayesian priors on parameter estimates, ridge regression, and lasso regression. Other algorithms that were assigned nontrivial weights included an elastic net regression model that favored the 1 (lasso) penalty over the 2 (ridge) penalty and a random forest with 500 trees. A variant of elastic net regression, random forests with 50 and 100 trees, and a multilayer perceptron all received low weights, while all other candidate algorithms in the library were assigned weights of zero.