Exposure proximal immune correlates analysis

Abstract

1 Introduction

An important goal in vaccinology is to understand the immunological mechanisms associated with the risk of disease and those responsible for protection, known as correlates of risk (CoR), and correlates of protection (CoP), respectively (Qin et al. 2007; Plotkin and Gilbert 2012). In randomized placebo-controlled trials, immune responses are measured shortly after the last dose of the vaccine, and these “peak” immune responses are correlated with disease outcomes using various methods. For example, peak antibody response measured 2 weeks after the last dose of vaccine is used as a covariate in a Cox regression model with time to disease as the outcome. Because disease cases are relatively rare, it is common to measure the immune response in all of the cases and subset of the controls, either by matching or by random sampling of controls. These are known as case–control and case cohort designs, respectively (Breslow et al. 1980; Prentice 1986). Peak immune response correlates analyses describe the early vaccine response on outcome and can be used as a surrogate endpoint for disease in regulatory settings (Gilbert et al. 2022).

However, peak immune response does not tell the complete story about vaccine’s effect on disease risk. Immune responses decay over time, and vaccine-induced protection often wanes. To better understand the time-varying efficacy of vaccines, one can use features of the immune response trajectory to predict disease risk or protection. Challenges include the infrequent and asymmetric sampling of immune responses in cases and controls, and the complexity of disentangling the effects of measured immune responses from unmeasured aspects and time variables. Moreover, vaccination can have a direct effect on disease risk that is not mediated through the measured immune response, and this direct effect itself may decay over time. Understanding the various aspects of a vaccine’s effect in relation to time is critical to guide the development and effective application of vaccine but has not been well studied. In this article, we aim to develop methods to assess the effect of decaying immune response on the risk of disease in the context of Cox regression with sparse sampling of immune response, in a baseline sero-negative population. Additionally, we develop modeling approaches that can help disentangle and investigate various components of a vaccine’s protective effect, utilizing immune response data from vaccine efficacy trials. The proposed model will be useful for prediction of vaccine efficacy.

There is a rich literature on modeling longitudinal covariates with survival data; see Rizopoulos (2012) and Tsiatis and Davidian (2004). Some researchers maximize the joint likelihood of the covariate process and the failure time process. For example, this can be based on a random effect model for the covariate process and survival time (De Gruttola and Tu 1994) or hazard (Wulfsohn and Tsiatis 1997; Rizopoulos et al. 2009) as a function of underlying random effect. Tsiatis et al. (1995) estimated the hazard conditional on covariate measured at the current time using a two-stage method that calibrates the covariate value conditional on observed data and then enters it into the Cox model. Tsiatis and Davidian (2001) developed a conditional score method that does not require distribution assumptions for random effects and obtains least square estimate of random effects based on the covariate trajectory of each individual separately. Recognizing that none of the above joint modeling approaches focused on designs with sub-sampling of longitudinal covariates nested from the full study cohort, Baart et al. (2019) developed a Bayesian estimation method that imputes the missing longitudinal data outside of the case-cohort design. Fu and Gilbert (2017) extended the full conditional score method to covariate sub-sampling, incorporating inverse sampling probability weighting (IPW) or augmented IPW. This work estimates the trajectory of each individual separately using individual-specific data and requires sufficiently repeated longitudinal samples to nonparametrically recover an individual’s trajectory.

Motivated by the need to analyze case-cohort immune correlates study nested within the COVID-19 randomized trials with peak immune response data, together with sparse sampling of longitudinal immune responses, we develop a method to model the effect of the decaying immune response, and/or the effect of the peak immune response and the decaying vaccine efficacy on the risk of disease in the context of Cox regression, adjusting for measurement error in immune response measurements. We consider a two-stage method that separately models the longitudinal trajectory and the time to event. Based on the estimated random effects model, we take the expectation of the hazard over the “missing” true immune responses at the peak time point and at each disease time, taking into account the measurement error in the measurement of immune response. We call the corresponding likelihood an estimated induced partial likelihood (IPL), maximization of which is like a single iteration of the expectation–maximization (EM) algorithm, where the expectation is over the hazard function. We further incorporate inverse sampling probability weighting into the estimated IPL to account for sub-sampling of immune correlates and then derive our estimator as the maximizer of this estimated weighted induced partial likelihood (WIPL).

While our work is most closely related to the work of Fu and Gilbert (2017) in addressing immune response sub-sampling for assessing association with disease risk and vaccine efficacy, there are several major differences. First, instead of using least squares modeling of each participant’s trajectory as in Fu and Gilbert (2017), we adopt the Gaussian distribution of random effects to accommodate the small number of time points with immune response measurements for each participant, as in the COVID-19 immune correlates studies. Second, our method targets the estimation of vaccine efficacy as a function of the decaying immune response among baseline seronegative (naive) populations, ie the population not exposed to the antigen at study enrollment. Among this population, the immune response among placebo recipients prior to infection is zero, and we utilize the placebo to improve the efficiency of the estimates. In contrast, the Fu and Gilbert method assumes the same random effect model between vaccine and placebo arms, thus, it is applicable to settings where placebo recipients have nonzero immune response, such as in the Dengue vaccine trials for which the method was applied. Moreover, the Fu and Gilbert method was developed for models with decaying immune response only while our method allows more extensive models that include additional components such as peak immune response as well and provides a novel dissection of the role of measured and unmeasured immune responses on disease risk. Lastly, the Fu and Gilbert method models baseline hazard at the scale of the study time (ie time since enrollment) while our method models the baseline hazard at the calendar time scale to accommodate the strong temporal trends in COVID-19 incidence. Our method applies not only to data from the standard vaccine trial period but also to data from the augmented closeout placebo vaccination (CPV) design (Follmann 2006) where placebo participants uninfected at the end of the standard trial period can be re-vaccinated and followed for infection outcome. This augmented design allows placebo recipients timely access to the vaccine when it is no longer proper to keep participants on placebo and allows more participants to contribute to immune correlates analysis (Follmann et al. 2021).

In Section 2, we introduce the problem setting, followed by descriptions of models and the WIPL estimators for immune response trajectories and time-to-event processes. We evaluate our methods through simulation studies in Section 3 and apply them to analyze immune correlates study from a phase III SARS-CoV-2 vaccine trial in Section 4.

2 Materials and Methods

2.1 Setting

We consider the problem setting of assessing immune correlates of risk and vaccine efficacy in a two-arm randomized vaccine trial among baseline-naive population. Let Z be the treatment indicator, with value 0 and 1 representing the placebo and vaccine arms, respectively. Let τ denote the peak time point after vaccination, U denote the time when infection outcome develops, and C denote the censoring time, all measured in days from the start of the clinical trial. Let $S=min(U,C)$ be the minimum of infection time and censoring time, and $δ=I(U≤C)$⁠, where I is the indicator function. We define $T=S−τ$ as the time to infection or censoring since the peak time point. Let W be a vector of baseline covariates. Let x(t) indicate the time-dependent biomarker process at time t after peak time point and let X(t) be corresponding observed biomarker value. For placebo recipients in the baseline-naive population, we have x(t) = 0 prior to infection.

We consider the setting of a two-phase sampling design for an immune correlate study. The first phase includes data from N independent participants in the trial, with information collected on Z_i, W_i, τ_i, S_i, δ_i, and $Ti=Si−τi$⁠. We accommodate staggered entry, allowing τ_i to vary across participants. We denote participants with $δi=1$ as cases and those with δ = 0 as controls, respectively.

The second phase includes a Bernoulli random sample of participants from phase 1, designed to measure an immune response biomarker at prespecified time points $tj:j=1,…,m$ following the peak time point. Let ξ_i be the binary indicator for whether participant i is sampled into the second phase, with the corresponding sampling probability given by $πi=P(ξ=1|Zi,Wi,τi,Si,δi)>0$⁠. For $j=1,…,m$⁠, let $X(tj)$ represent the immune responses measured at times t_j, since the peak time point, for the selected controls and all cases from the vaccine arm. For cases, we additionally obtain an immune response at the time of disease detection. Denote this value by $X(Ti)$ for case i with $Ti=Si−τi$⁠.

Immune responses for vaccine recipients at peak time point (t = 0) are typically measured during the second phase. Our method integrates Cox regression modeling in the immune correlate study with the estimation of immune response trajectories based on sparse longitudinal data (usually not more than three or four per participant). The longitudinal sample may be a part of the immune correlate study or, as in the case study, an external sample.

Moreover, we consider designs incorporating a possible augmented CPV component (Follmann 2006) in the study: at the end of the standard trial follow-up period, a fraction of placebo recipients who remain uninfected may receive the vaccine and then be followed up for infection outcome. Among those participants receiving the closeout placebo vaccination, immune response biomarkers are measured at a set of prespecified time point from a random sample and from all cases post-CPV. For cases, immune responses can again be measured at the time of infection. Let ξ_k be the sub-sampling indicator for the $kth$ participant from the CPV component, where the corresponding sampling probability is $πk=P(ξ=1|Zk,Wk,τk,Sk,δk)>0$⁠. Here, Z_k = 1 indicates vaccination, τ_k represents the peak time after re-vaccination, S_k is the minimum between the calendar time to event and the last follow-up time, δ_k is the censoring indicator post-CPV, and W_k represents the baseline covariate measured at study enrollment. This CPV component has been incorporated into multiple COVID-19 trials (Follmann et al. 2021), making it crucial to evaluate a method that includes data from CPV.

2.2 Model

Evaluating the role of vaccine-induced immune responses on disease risk is complex, as vaccination induces a suite of immune responses, yet only a few features are measured. Antibodies are produced in addition to long-lived memory B and T cells, which remain dormant until they are quickly mobilized for an anamnestic response following fresh exposure to the antigen, such as SARS-CoV-2. Thus, the circulating antibody level, which we measure by $X()$ is only one aspect of the protective effect of vaccination. Furthermore, the decay in circulating antibodies over time generally occurs faster than that of cell-based features (Goel et al. 2021). Therefore, the true vaccine effect depends on both seen and unseen factors that diminish at different rates. Given the limited data and the fact that we only measure $X()$⁠, incorporating this complexity is challenging.

Here, u represents the calendar time to infection, τ is the peak time point after vaccination, and $x(u−τ)$ is the true antibody response at time $u−τ$ post-peak (namely, the “exposure-proximal” correlate). This model captures the effect of the immune response at the time of antigen exposure on the risk of infection. We use calendar time as the time index for this model, rather than time since randomization, to better address secular variation in the attack rate due to waves or seasonality.

While this model is attractive in terms of dissecting mechanisms, model (3) may be statistically difficult to fit since $x(t),x(0)$⁠, and t are all correlated with each other for $t=u−τ$⁠. Furthermore, with an unspecified baseline hazard, it may be difficult to identify β₃. Indeed, for a simultaneous entry trial $u−τ$ is equivalent to u, the term $λ0(u) exp {β3(u−τ)}$ is eliminated from the partial likelihood. The same lack of identifiability occurs if we model the hazard as a function of t (the time post-peak), rather than calendar time.

We note that, as always, the model is relevant only within the observed range of covariates, ie for x(t) within the observed range, and thus interpretation of the $1−exp(β0)$ may not be meaningful if x(t) does not approach 0.

The coefficients of various interaction terms in the Cox model are related to the contributions of different aspects of the vaccine’s effect. Specifically, nonzero β₁, β₂, and β₃ indicate the effects of the exposure-proximal immune response, the peak immune response, and the vaccine’s decay over time through mechanisms other than those specific immune responses. Figure 1 illustrates the model under four scenarios, using VE over time t since the peak time point. The three curves correspond to individuals with low, median, and high antibody responses that decay over time. The model is rich in terms of allowing various aspects of the immune system to influence risk. These models can support decisions on whether to administer booster shots. Figure 2 displays two different models for waning vaccine efficacy. The thick solid curve represents a model with $β1=β2=0$⁠, indicating that waning efficacy depends on unmeasured factors independent of antibodies. By day 292, the VE has dropped to 0.50, suggesting a recommendation for everyone to get boosted around 292 days would be sensible. In contrast, the dashed curves represent a model with $β2=β3=0$⁠, where the drop in vaccine efficacy depends on decaying antibody levels, which vary from person to person. The thick dashed curve denotes an individual at the antibody level at the 50th percentile while the dashed lines represent antibody levels at the 10th and 90th percentiles. These individuals have their VE drop to 0.50 on days 219, 318, and 416, post the peak time point, respectively. Using an antibody-based trigger for boosters would be more precise but also more complex to implement and would require frequent testing. Nonetheless, it might be a compelling strategy for some, especially those who with a poor response to vaccination.

Fig. 1

Four scenarios for the effect of antibody on the risk of disease. The top two panels denote $β2,β3=0,0$ (left) and $β1,β3=0,0$ (right). The bottom panels denote $β3=0$ and $β1,β2=0,0$⁠. The solid line denotes the median immune response which decays over time, while the dashed lines denote individuals with higher and lower responses.

Open in new tabDownload slide

Fig. 2

Estimated vaccine efficacy over time for two models. For the thick solid curve, VE depends on time since vaccination, for the black curves VE depends on each individual’s decaying antibody. The three dashed curves denote someone at the 10th, 50th, and 90th percentiles of the antibody distribution. The vertical lines denote when a person would be recommended for a booster dose as VE has gone below 50%.

Open in new tabDownload slide

2.3 Estimation

Equations (1) and (4) specify a joint model for a covariate process and a time to event process. Considering missingness due to both unobserved underlying true immune responses and participants not-sampled in the immune correlates study, there are different ways to approach estimation. These include the maximization of a joint model or estimated partial likelihood type approaches (Zhou and Pepe 1995), where the parameters of the covariate distribution are first estimated and used in the likelihood for the time-to-event process. In this work, we use the latter framework by first estimating the parameters of the immune response process and then estimating the expectation of the hazard conditional on observed data to replace the hazard in the partial likelihood. This two-step estimation procedure allows for the use of either data from participants in the immune correlate study or an external dataset to estimate the immune response process.

First, we obtain $Θ^$⁠, the estimator of Θ with linear mixed effect model using longitudinal immune response data. These data are exclusively from vaccine recipients and considers time points from the “peak” time until just before the time of infection. Note that the longitudinal sample could be from our second-phase immune correlate study nested within the vaccine trial, or it could be an external data as in the case study.

More details about the functional form of (9) can be found in Supplementary Appendix A. By substituting equation (9) into $WIPLv(β)$ and $WIPLVP(β0,β)$ as specified in (7) and (8), we can obtain the estimators for $β$ or $(β0,β)$ by maximizing corresponding likelihood.

Wang et al. (2001) studied an estimated induced partial likelihood estimator in a simpler problem setting, calibrating the true underlying univariate covariate with its surrogate, measured with error, and demonstrated asymptotic normality of the estimators under standard regularity condition. Here, we use nonparametric bootstrap sampling to estimate the standard error of the WIPL estimators and to construct corresponding confidence intervals.

Note that if immune response measure has little measurement error, ie $σ2$ from equation (1) is close to 0, the second term in equation (9) above is nearly zero, and the proposed approach approximately follows a regression calibration method for addressing measurement error (Prentice 1982). In this method, the covariate with measurement error (ie X_i) is simply replaced by the expectation of the underlying truth, ie $E(x|Θ^,Vi)$⁠. However, as demonstrated in Wang et al. (2001) and corroborated by our simulation studies, the regression calibration estimator can be substantially biased in extreme cases; in contrast, the bias in WIPL estimators is minimized through the inclusion of a second-order approximation.

Maximization of the estimated weighted induced partial likelihood resembles a single step of the EM algorithm, where the expectation is calculated over the hazard. With a rare disease, estimation based on a fully specified joint model of the time-to-event and biomarker processes should be quite similar to the approach described above, as the estimates of the biomarker process are largely unaffected by the time to event information.

In Supplementary Appendix B, we sketch the Fu and Gilbert (FG) (Fu and Gilbert 2017) estimator, a comparative method proposed to estimate the infection hazard in model conditional solely on the exposure-proximal correlate x(t). There are several key differences between the FG estimator and our proposed estimator: (i) The FG estimator was developed on the study time scale, rather than the calendar time scale as in our estimator. (ii) The FG method employs a least squares procedure to estimate the immune response trajectory for each individual separately. Therefore, it is only applicable to studies where the longitudinal immune responses are collected from participants within the immune correlates study sample for hazard modeling, rather than from an external sample. (iii) The FG estimator assumes that immune responses over time across study arms follow the same time effect model with same measurement error structure. This assumption is reasonable in settings where participants exhibit a positive immune response at baseline (e.g. in the Dengue Vaccine study) but is unsuitable for the baseline-naive study settings we are examining, where participants in the placebo arm have zero immune responses.

2.3.1 Accommodating CPV design

Our proposed WIPL methods for standard vaccine trials, as discussed in previous section, can be readily generalized to accommodate the augmented study design with a CPV component. The CPV design contributes the following observations for construction of the partial likelihood:

• Suppose there are N_Vparticipants assigned to the vaccine arm at the beginning of the trial. For vaccine recipient $i∈1,…,NV$⁠, we observe $Zi=1,τi,Si,δi,ξi,ξiXi, Ti=Si−τi$⁠, and $π^i$ as in a standard trial. We denote this set of observations as $V$⁠.

• Suppose there are N_P participants assigned to the placebo arm at the beginning of the trial. For each placebo recipient $j∈NV+1,…,NV+NP$⁠, let C_j represent the planned re-vaccination time if uninfected (which serves as the censoring time). We observe $Zj=0,τj,Sj,δj,ξj,ξjXj$ with $ξj=1,Xj=0$ and $π^j=1$⁠. Here, S_j is the minimum between the event time and C_j, $Tj=Sj−τj$⁠, and δ_j indicates whether an infection occurs before C_j. We denote this set of observations as $P$⁠.

• Suppose there are N_C participants receiving CPV, contributing an additional N_C observations for $k∈{NV+NP+1,…,NV+NP+NC}$⁠. The $kth$ observation includes $Zk=1,τk,Sk,δk,ξk,ξkXk$ and $π^k$⁠, where τ_k corresponds to peak time point after re-vaccination, S_k is the minimum between the event time and follow-up time during the CPV period, $Tk=Sk−τk$⁠, δ_k is the censoring indicator during the CPV period, ξ_k denotes whether the participant was subsampled for immunogenicity measurement $Xk$ during the CPV period, and $π^k$ represents estimated sampling probability for immunogenicity measurement in CPV. We denote this set of observations as $C$⁠.

2.4 Adjustment for ties

Previously we considered estimation under the assumption that no ties are present among the failure times, allowing the data to be uniquely sorted by failure time. However, in practice, ties are often present, primarily due to failure time being reported only to the nearest day, as seen in our motivating Moderna trial data. In standard Cox regression, various approaches to adjust for ties in failure time have been considered (Kalbfleisch and Prentice 2011). Maximizing the average of the partial likelihood over all possible ways of breaking the ties is intuitive but computationally challenging. Two common approximations of the average likelihood in the standard Cox model setting are: the Breslow approximation (Breslow 1974), which has the simplest form but results in a consistently larger denominator in the partial likelihood than appropriate; and Efron’s approximation (Efron 1977), which offers improved accuracy compared to Breslow’s approximation and is moderately easy to work with.

3 Simulations

The immune correlates study has a random 1,500 vaccinees and all vaccine cases, for which antibody levels are measured at three time points (day 0, 6 months, 1 year post-peak time point) or four time points (the three previous times plus month 3 post-peak time point) prior to the follow-up date. For cases, antibody levels were also measured at the time of infection.

Supplementary Table 1 compares the β₁ estimators based on vaccine recipients data among HF, FG, NRC for several scenarios. For HF and NRC, the LME model fitting utilizes all measures prior to infection, with each vaccine recipient weighted by the inverse probability of being sampled for immune response measurement, as implemented in the R WeMix package. Immune responses for placebo recipients (without re-vaccination) are set to zero. The reference scenario includes a random intercept and slope effect (std(a₁)=0.0015) with $σ=0.145$⁠. In this scenario, immune responses are measured at three time points: day 0, 6 months, and 1 year post-peak time point. We also explore variation in the scenarios that feature increased measurement error (⁠$σ=0.3$⁠), random intercept only (ie std(a₁)=0), and sparsity of time point sampling. This sparsity includes measurements at only day 0 and month 6 for 50% of vaccinees in the immune correlates study, or having longitudinal measurements at four time points. Based on vaccinee data only, both the HF estimator and the FG estimator show minimal bias in the β₁ estimate, whereas the NRC estimator exhibits larger bias (e.g. 6.2% in the reference scenario, which further increases to 11.5% when the measurement error increases to $σ=0.3$⁠). The HF estimator is more efficient compared to the FG estimator, particularly with increased sparsity in longitudinal time points. For instance, HG is 14% more efficient than FG in the reference scenario. The efficiency gain increases to 25% when 50% vaccinees in the immune correlates study lacks a year 1 post-peak measurement, but drops to 5% when the immune correlates set includes measurement at four time points. A significant efficiency gain for HG relative to FG is observed when the observed immune response trajectory follows a random intercept model: 65% and 115% when the immune correlates set includes all or half measuring immune response at year 1 post-peak, respectively.

In Supplementary Appendix C (Supplementary Tables 1–3), we present further comparison results between the HG and NRC estimators in settings that allow for staggered study entry with uniform accrual over 4 months; the analysis occurs 15 months after the start of the study. These comparisons are based on data from either the vaccine-only or the vaccine+placebo estimators, under various scenarios with measurement error σ set at 0.145, 0.3, and a hazard ratio model including x(t) or x(0) and t in addition. While the proposed HF estimator remains approximately unbiased in all models, the NRC estimator—whether based on the vaccine alone or the vaccine+placebo—can exhibit significantly larger bias, particularly as the measurement error increases. For instance, with longitudinal measures at three time points in a model including x(t) and t, the vaccine-only NRC estimator’s biases for β₁ and β₂ are 7.3% and –7.5%, respectively when $σ=0.145$⁠, which escalates to 15.4% and –13.1% when $σ=0.3$⁠; the vaccine+placebo NRC estimator’s biases for β₀, β₁, and β₂ are 1%, 3.12%, and –2.61, respectively when $σ=0.145$⁠, changing to –15.01%, 5.19%, and 0.75% when $σ=0.3$⁠.

We further investigate the proposed HF estimators in more details by comparing analysis models with different levels of complexity in a study with staggered entry. Initially, we set $β1=β2=β3=0$ and explore three distinct designs:

• Design A: a standard trial with $λ0(u)=0.116$ and $β0=−3$ such that on average there are 1,778 and 94 cases from placebo and vaccine arm, respectively at months 15 after the start of the study, with longitudinal measurements at day 0, month 6, and year 1 postpeak time point.

• Design B: a standard trial with $λ0(u)=0.015$ and $β0=−2.2$ such that on average there are 242 and 27 cases from placebo and vaccine arm, respectively at months 15 after the start of the study, with longitudinal measurements at day 0, month 6, and year 1 postpeak time point.

• Design C: a design of the same follow up as Design A but with a CPV component at 9 months after enrollment (expecting VE to show around month 9), all un-infected placebo recipients are re-vaccinated and followed until reaching the 15 month after study enrollment. Among CPV participants, 1,500 random sampled participants and all cases have measurement of antibody at day 0 and at 3 months post the peak time point after the re-vaccination.

The top part of Table 1 presents the results from the reference scenario, under Design A. This scenario models the coefficients β₁, β₂, and β₃ together, and only uses X(t) value at prespecified time points for model fitting, ie explicitly excluding X(T) information from cases. For both the vaccine-only or vaccine+placebo approaches, the estimates for all coefficients appear unbiased. Comparing the variance rows of scenario 1, we see that there is no efficiency gain using the placebo group for estimating β₁ and β₂, the coefficients of $x(t),x(0)$⁠, respectively. However there is a substantial efficiency gain for estimating β₃, the coefficient for t. This gain is understandable as as without a placebo group, it is challenging to distinguish the baseline hazard $λ0(u)$ from $exp(β3t)$⁠. Specifically, if all participants enroll at the same time and τ is uniform across the cohorts, β₃ becomes nonestimable. The inclusion of the placebo group significantly aids in separating these components, illustrated by the relative efficiency is 1.707/0.453 = 3.77.

In finalizing the model, it is conceivable that β₃ will be near zero, suggesting that a reduced model incorporating only x(0) and x(t) should be evaluated. Thus, the second scenario focuses on estimates based on this reduced model for Design A. Compared to the reference scenario (which includes β₃), and utilizing all available data, significant efficiency gains are observed with 0.158/0.041 = 3.85 for β₁ and 0.326/0.190 = 1.72 for β₂, respectively.

In Scenario 3, Design A is simplified even further by estimating only β₁. This model reduction focuses exclusively on the impact of x(t). The efficiency gains from this further reduction is 0.041/0.030 = 1.37.

Scenario 4 explores the utilization of the proximal antibody measurement, X(T), from cases under the assumption that X(T) is measured immediately prior to exposure. Compared to the reference scenario and using data from both vaccine and placebo arms, there are modest gains in efficiency of $0.158/0.125,0.326/0.290,0.453/0.396=1.26,1.12,1.14$⁠, for $β1,β2,β3$⁠, respectively.

In Scenario 5, we consider Design B, which has fewer events compared to Design A. While still maintaining minimal bias, this scenario shows lower efficiency compared to the reference scenario.

Scenario 6 considers Design C, where a CPV component is added to Design A. This modification induces a different distribution of antibody trajectories over time and increases the number of vaccine cases while decreasing the number of placebo cases. We observe similar variances between the vaccine and vaccine+placebo settings for the estimates of $β1,β2$⁠, and β₃, with a modest efficiency gain of 0.461/0.430 = 1.07 for the estimation of β₃ when placebo recipients are included. Compared to the reference scenario, the most notable difference is in the variance of $β^3$⁠, which goes from 1.707 to 0.461 when only data from the vaccine arm is used. Adding the CPV component also leads to appreciable efficiency gains of vaccine-only estimators for reduced models. For example, the efficiency gain for the estimation of β₁ using only vaccinees is 0.128/0.042 = 3.05 for models excluding β₃ (ie Scenario 7 vs Scenario 2) and 0.063/0.029 = 2.17 for models with β₁ only (ie Scenario 8 vs Scenario 3). Overall, having a CPV component leads to a slight efficiency gain for methods utilizing both vaccine and placebo recipients.

Some conclusions can be drawn from the simulation.

• Incorporating the placebo group significantly improves efficiency. For the model including β₁ and β₂, the variance ratio for $β^1$ is 0.128/0.041 or 3.1 in Design A.

• It is challenging to estimate all three parameters $(β1,β2,β3)$ simultaneously, even with the inclusion of the placebo group. The variance ratio for $β^1$ in a three parameter model compared to the model with β₁ alone is 0.158/0.030 or 5.3 in Design A.

• Adding the antibody response upon presentation of disease, X(T), results in a modest increase in efficiency. However, given the potential for significant changes in X(T) that could introduce larger bias, and considering the minor efficiency gain, we do not recommend including X(T) for modeling purpose.

• Including antibody measurements after a placebo crossover leads to a modest gain in efficiency, primarily due to the increased number of vaccine cases with measured immune responses.

Next, for settings with non-null x(t) and/or x(0) and t effects, we investigate performance of the HF estimators under working models with varying level of complexity. Specifically, (i) for settings with a non-null x(t) effect (nonzero β₁), we assess the estimation using a model with x(t) alone, and models with x(t) and x(0) or x(t) and t. (ii) For settings with non-null x(t) and x(0) effect (nonzero β₁ and β₂), we assess the estimation using model with $x(t),x(0)$⁠; models with x(t) or x(0) alone; and a model with $x(t),x(0),$ and t. (iii) For settings with non-null x(t) and t effects (nonzero β₁ and β₃), we assess the estimation using a model with x(t) and t; models with x(t) or t alone; and a model with $x(t),x(0)$⁠, and t. Table 2 and Supplementary Table 4 present the bias and variance for model parameter estimates using vaccine-only data and vaccine+placebo data in a standard trial, with immune responses measured at three and four pre-specified time points, respectively. Table 3 and Supplementary Table 5 further present performance for estimating $VE(x(t),x(0),t)$⁠, the vaccine efficacy corresponding to various combination of x(t), x(0), and t values. Overall, the proposed HF estimator exhibits minimal bias in estimating model parameters and VE when important covariates are included, although the inclusion of additional covariates with null effect decrease the precision of the estimates. On the other hand, the omission of important covariates can lead to severe bias in VE estimates. Similar observations can be made for HF estimators based on data with a CPV component (Supplementary Tables 6 and 7).