The immunogenicity sub-study of a double-blind, placebo-controlled efficacy trial of a high-potency live-attenuated herpes zoster (shingles) vaccine (zoster vaccine live) enrolled 1395 subjects. There were 32 cases of herpes zoster (HZ) between vaccination and end of follow-up. Blood samples obtained at baseline (prior to vaccination), 6 weeks after vaccination, and at 1, 2, 3 years thereafter were tested (among other assays) for VZV antibody titers by glycoprotein ELISA (gpELISA).

Information needed for this analysis (vaccination status, disease status, time to disease or end of follow-up, age, and VZV antibody titers prior to vaccination and at week 6) was available for 1326 subjects. These included 655 zoster vaccine recipients and 671 placebo recipients; 773 individuals were 60–69 years of age, 553 individuals were 70–93 years of age; and 32 individuals experienced HZ after the week 6 visit. Time-to-event was determined for each subject as the time after the week 6 visit to either HZ (for cases) or end of follow-up (for non-cases). Median time-to-HZ was 714 days (time range, 63–1191 days), median time-to-event (including both cases and non-cases) was 949 days (time range, 17–1259 days). Fold rise in VZV antibody titers is derived from the ratio of gpELISA measurements at week 6 divided by that prior to herpes zoster vaccination. Subjects with values of fold rise less than 1 are assigned fold rise value of 1 (log2 fold rise equal to 0; any reduction of assay value is likely the result of measurement error). All placebo recipients are also assigned fold rise value of 1 (log2 fold rise equal to 0; the control group has, by definition, no change of assay; any measured negative or positive change of assay value is likely the result of measurement variability).

Post-vaccination immunogenicity responses are visualized as boxplots by HZ status, vaccination status, and age group. The two-sided Mann-Whitney U test is performed to assess the difference between immunogenicity of herpes zoster vaccine in younger and older subjects.

Anonymized data analyzed in this work were published and described in detail elsewhere29. Briefly, the ancillary study to the phase 3 CYD-TDV trial in Asia (CYD-14)30 was conducted in Cebu, Philippines; a cohort of participants (\(n=\) 611, consisting of 417 CYD-TDV recipients and 194 placebo recipients; mean age, 8 years at baseline; age range, 2–14 years; 90% individuals were seropositive at baseline) were followed for more than 6 years after their third dose29. There were 87 VCD (DENV1, \(n=\,\)12; DENV2, \(n=\,\)43; DENV3, \(n=\,\)16; DENV4, \(n=\) 6; unknown serotype, \(n=\) 10) between the third dose of the vaccine and end of follow-up. Plaque reduction neutralization test (PRNT50) was conducted on blood samples collected shortly post vaccination (around day 28) and then annually and during symptomatic infections29,55,56.

In this work, for six individuals experiencing multiple VCDs (which are here considered competing risks) during the observation period (between the day 28 visit and the end of follow-up), only the first VCD is used. There were 81 VCD cases (DENV1, \(n=\,\)11; DENV2, \(n=\,\)40; DENV3, \(n=\,\)15; DENV4, \(n=\) 5; unknown serotype, \(n=\) 10) in the resulting dataset for models with DENV-Any endpoint. Time-to-event was determined for each subject as the time after the day 28 visit to either VCD (for cases) or end of follow-up (for non-cases). Median time-to-VCD was 385 days (time range, 1–2291 days), median time-to-event (including both cases and non-cases) was 1812 days (time range, 1–2291 days). Serostatus is a binary indicator of pre-vaccination exposure to at least one serotype of DENV. Immunogenicity biomarkers used in this analysis are serotype-specific log2 PRNT50 values measured around day 28 post the third dose of CYD-TDV (median time, 30 days; time range, 19–41 days). These log2 PRNT50 values were available for 564 subjects (392 CYD-TDV recipients and 172 placebo recipients). Average titer is a derived variable, a single summary biomarker that can be assessed as an immune correlate for DENV-Any and DENV1-4 endpoints. When serotype-specific PRNT50 titers are correlated (as here, see Supplementary Fig. 1), the average titer will typically have greater precision than analyses using models specific to each serotype. Average titer is determined for each subject as the arithmetic average of respective four centered and scaled serotype-specific log2 PRNT50 titers. For each serotype, first, centering is performed by subtracting the mean log2 titer (across all subjects) from each log2 titer observation (per subject). Next, scaling is performed by dividing each (centered) log2 titer observation (per subject) by the standard deviation of log2 titers (across all subjects).

Post-vaccination immunogenicity responses are visualized (Fig. 3, Supplementary Fig. 2) as boxplots by serotype-specific VCD status, vaccination status, pre-vaccination serostatus and serotype. The two-sample two-sided t-test is performed to assess the difference between immunogenicity of CYD-TDV in seropositive and seronegative subjects. The assumption of normality was satisfied for DENV1, DENV3, and DENV4, as confirmed visually using quantile-quantile plots and verified through the Kolmogorov-Smirnov test for both seropositive and seronegative groups. To assess the homogeneity of variances between groups, the two-sided Bartlett’s test was conducted. In the case of DENV4, where the assumption of homogeneity of variances was violated, the two-sided Welch’s t-test was employed as an alternative to the standard t-test. Additionally, for DENV2, where the assumption of normality was violated in the seropositive group, the two-sided Mann-Whitney U test was applied.

An immune response biomarker is a correlate of risk (CoR) if probability of disease is statistically associated with biomarker levels. In the terminology defined by Qin et al.57, CoR uses data exclusively from subjects who have received active vaccination. However, in our analyses, we use data from all subjects (i.e., including those who received placebo control) in time-to-event models to evaluate immunogenicity as a CoR8,9. We thus aim to provide more comprehensive and precise assessment of the association between immunogenicity and the risk of disease.

Where there are multiple immune response biomarkers (serotype-specific log2 PRNT50 values for the dengue CYD-TDV vaccine), multicollinearity is assessed using VIFs (variance inflation factors) and pairwise Pearson’s correlation coefficients. Use of highly correlated immunogenicity predictors in regression is avoided by transforming them into new predictor(s): a single summary biomarker (e.g., average log2 PRNT50 titer for the dengue CYD-TDV vaccine), or principal component(s).

In Cox PH models, given a set of \(L\) predictors \(_,\,_,\,\ldots ,\,_}\) for subject \(i\), the hazard function \(\lambda (t}_)\) has the form:

$$\lambda \left(t|_\right)=_\left(t\right)\mathrm\left(__+__+\cdots +__}\right)=_\left(t\right)\mathrm\left(\right)$$

(1)

with baseline hazard \(_\left(t\right)\) identical for all subjects (has no dependency on \(i\)), and linear predictor \(=__+__+\cdots +__}\). Parameters of the model \(_,\,_,\,\ldots ,\,_\) can be estimated by maximizing the partial log-likelihood58.

Alternatively, a hazard function involving an interaction term (denoted, e.g., \(_\)) between, e.g., predictors \(_\) and \(_\) may be described as

$$\lambda \left(t|X\right)=_\left(t\right)\mathrm\left(+___\right)$$

(2)

The effect of immunogenicity (\(T\)) can be estimated using either a linear term,

$$\lambda \left(t|T\right)=_\left(t\right)\mathrm\left(_T\right)$$

(3)

or a non-linear term (or terms), e.g.,

$$\lambda \left(t|T\right)=_\left(t\right)\mathrm\left(_\log (T)\right)$$

(4)

$$\lambda \left(t|T\right)=_\left(t\right)\mathrm\left(_\sqrt\right)$$

(5)

$$\lambda \left(t|T\right)=_\left(t\right)\mathrm\left(_T+_^\right)$$

(6)

We define the risk curve, \(\rho (T)\), as hazard ratio, a function of immunogenicity (\(T\)), representing the hazard relative to the reference datapoint (\(_}\)).

$$\rho \left(T\right)=\,\frac_}\right)}$$

(7)

In this work, models using the relevant immunogenicity biomarker are first fitted with linear and non-linear term(s) (Eqs. 3 and 6, resp.). The only predictor in these models is the immunogenicity biomarker, \(T\). The best-fitting model is selected based on Akaike information criterion (AIC)59. Next, the resulting model is fitted adjusting for all pre-selected potentially clinically meaningful covariates4,5,60,61,62 and their interactions. The final model for CoR assessment is selected based on AIC (unadjusted or adjusted). The final model’s hazard ratio is visualized as one overall risk curve (for unadjusted model, e.g., for the zoster vaccine) or multiple risk curves (for adjusted models; one per each covariate realization value, e.g., seropositive and seronegative for the dengue CYD-TDV vaccine, or one per each combination of covariates’ realizations) with pointwise 95% CI. The studied immunogenicity biomarker is deemed a CoR if the coefficient involving the biomarker in the final model is different from 0 at pre-specified level of statistical significance (here we adopt \(\alpha =0.05\)).

In the context of competing risks63 (e.g., virologically confirmed dengue disease caused by DENV1-4), the effect of immunogenicity on time-to-event outcome is evaluated analogously using cause-specific Cox PH models31,32 and Fine-Gray subdistribution hazards models33, as further explained below.

The cause-specific Cox PH model has the form of:

$$_\left(t|_\right)=_}}\left(t\right)\mathrm\left(_}}_+__+\cdots +_}_}\right)$$

(8)

This PH model of event type \(c\) at time \(t\) allows effects of the covariates to differ by event types (e.g., DENV1-4), as the subscripted beta coefficients suggest. The cause-specific approach has a limitation in that it assumes noninformative censoring for subjects who experience competing events. This assumption is equivalent to saying competing events are independent, and the sum of individual event probabilities is generally not constrained to 1. However, there is no way to explicitly test whether this assumption is satisfied for any given dataset.

The Fine-Gray subdistribution hazards model aims at modeling the hazard function derived from a cumulative incidence function (CIF, also known as subdistribution hazard):

$$_}\left(t|_\right)=-d\mathrm\}_\left(t|_\right)\}/$$

(9)

The CIF, accounting for competing risks, estimates the marginal probability for each competing event as a function of its cause-specific probability and overall survival probability (i.e., ensures the probabilities add up to 1). The CIF for event type \(c\) at time \(_\) is the cumulative sum up to that time of incidence probabilities over all times for the events of interest (i.e., type \(c\)), which is expressed as:

$$}_\left(_}_\right)=\mathop\limits_^\hat_}\left(_}_\right)=\mathop\limits_^\hat\left(_}_\right)\times \hat_}(_}_)$$

(10)

where \(\hat\left(_|_\right)\) denotes the estimate of overall (any event) probability of surviving \(_\) (the last event-of-interest time before \(_\)) given predictor vector \(_\) and \(\hat_}(_|_)\) represents the estimate of hazard at ordered event time \(_\) for the event-type of interest given \(_\).

This approach overcomes the limitation of the cause-specific Cox PH approach by taking into account the informative nature of censoring due to competing risks (i.e., events other than the one of interest which alter the probability of experiencing the event of interest). The CIF models developed by Fine and Gray are analogous to the Cox PH model:

$$_}\left(t|_\right)=_}}\left(t\right)\mathrm\left(__+__+\cdots +__}\right)$$

(11)

The results from fitting these models have a similar interpretation regarding the effects of predictors in the model as can be derived from the Cox PH model approach for competing risks data.

An immune response biomarker that is predictive of VE is termed a correlate of protection (CoP)7. Here, the Prentice criterion24 of conditional independence is applied to test if immunogenicity fully mediates the vaccine effect, i.e., if the effect of vaccination status on time-to-disease endpoint is insignificant after controlling for the impact of immunogenicity.

To assess a CoP, an additional predictor, the vaccination status, is added to the final CoR model (selected as described in the previous section). The studied immunogenicity biomarker is deemed a CoP if the coefficient involving the biomarker remains significantly different from 0 and the coefficient involving the vaccination status is not significantly different from 0 at pre-specified level of statistical significance (we adopt \(\alpha =0.05\)).

For an efficacious vaccine, a CoP should be also a CoR. However, as VE can differ across subgroups, a vaccine that is efficacious in some subgroups can also fail to be efficacious in others. In such a case, an immune biomarker could be a CoR in some subgroups only (and fail to be a CoR in others), but if it fully mediates the vaccine effect in the covariate-adjusted CoR model, we call it a CoP.

For illustration, assume that only one binary baseline covariate (e.g., age group or serostatus), \(X\), is considered and VE is to be estimated in two covariate-defined subgroups (as in the presented examples). Let \(_^}\), \(_^}\) be the immunogenicity biomarker measurement for \(i\)-th vaccinated subject (\(i=1,\ldots ,\,I\)) in the first covariate-defined subgroup (where \(X=0\)) and \(j\)-th vaccinated subject (\(j=1,\ldots ,\,J\)) in the second covariate-defined subgroup (where \(X=1)\), respectively. Let \(_^}\), \(_^}\) be the immunogenicity biomarker measurement for \(m\)-th control subject (\(m=1,\ldots ,\,M\)) in the first covariate-defined subgroup (where \(X=0\)) and \(n\)-th control subject (\(n=1,\ldots ,\,N\)) in the second covariate-defined subgroup (where \(X=1\)), respectively. The immunogenicity-based VE in subgroups of interest can be estimated using the subgroup-specific risk curves, \(\rho \left(T\right|X=0)\), \(\rho \left(T\right|X=1)\), obtained by covariate-adjusted models (Eqs. 1 and 2) and subgroup-specific immunogenicity data, \(_^}\), \(_^}\), \(_^}\), \(_^}\), as:

$$}\left(T|X=0\right)=(1-})\cdot 100=\left(1-\frac_^\rho (_^}}X=0)}_^\rho (_^}}X=0)}\right)\,\cdot\, 100$$

(12)

$$}(T|X=1)=(1-})\cdot 100=\left(1-\frac_^\rho (_^}}X=1)}_^\rho (_^}}X=1)}\right)\,\cdot\, 100$$

(13)

The VE point estimates in subgroups are obtained using the coefficient values, e.g., \(_,\,_,\,_},\,\ldots ,\,_\), that were fit to the original observed time-to-event, immunogenicity, and demographic data. For a given set of data, the 95% CI associated with estimated VE needs to account for the uncertainty regarding the \(_,\,_,\,_},\,\ldots ,\,_\) parameters and variability in the observed data. This can be done via parametric resampling of the posterior distribution for parameters and bootstrapping the observed data in the vaccinated and control groups. The bootstrap resampling of observed data is performed on subjects (1000 times): each time a subject is selected, his/her immunogenicity value and all his/her demographic characteristics (covariate values) are used (to account for the titer distribution and its uncertainty) in the estimation of VE. To account for uncertainty in the model parameters, for each of the 1000 bootstrapped data sets a set of model parameter values is randomly drawn from the parameters’ probability distribution (a multivariate normal distribution characterized by the coefficients’ point estimates and the covariance matrix derived from the fit). The 2.5th and 97.5th percentiles of these 1000 VE values are used to establish the 95% CI.

Further information on research design is available in the Nature Research Reporting Summary linked to this article.