Simultaneous Inference of Treatment Effect Modification by Intermediate Response Endpoint Principal Strata with Application to Vaccine Trials

1 Introduction

In vaccine research, identifying biomarkers that can be used as surrogate endpoints for clinical endpoints is an important question to address. A good surrogate can be used to guide the development of the vaccine and predict the vaccine’s protective effect on the clinical endpoint in future settings before conducting efficacy trials. The research in this manuscript is motivated by the need to evaluate immune response biomarkers as effect modifiers of a vaccine’s effect on the clinical endpoint of interest (i. e. clinical vaccine efficacy), which is one way to study biomarkers and develop their utility for predicting vaccine efficacy. Such post-randomization effect modification research has been referred to as principal surrogate evaluation by Frangakis and Rubin [1].

Various research has been conducted for evaluating surrogates under the principal stratification framework [2, 3, 4, 5]. For evaluating continuous immune response biomarkers as principal surrogates in vaccine trials, Gilbert and Hudgens [2] proposed the marginal Causal Effect Predictiveness (mCEP) curve as an estimand for principal surrogacy. The mCEP curve is defined in terms of clinical risks conditional on the potential biomarker if assigned to vaccine and it quantifies how well causal treatment effects on the biomarker predict causal treatment effects on the clinical endpoint. Based on the mCEP curve, a useful biomarker is one that is a strong effect modifier, where the mCEP varies over subgroups defined by the biomarker levels if receiving vaccine [6].

To estimate the mCEP curve in the presence of missing counter-factual biomarker values if receiving vaccine in an efficacy trial, Huang, Gilbert, and Wolfson [7] proposed a pseudo-score type estimator (inspired by Chatterjee, Chen, and Breslow [8]). This estimator utilizes a baseline predictor associated with the biomarker if receiving vaccine as an auxiliary variable to estimate the parameters in the assumed parametric risk model. Huang et al. [7] developed a procedure for inference about the mCEP curve under the condition that the risk model is independent of the baseline predictor after conditioning on the biomarker if receiving vaccine. This condition may be reasonable in some settings. For example, in a HIV vaccine efficacy trial where the biomarker of interest is an immune response to HIV, the trial can be designed to inoculate everyone with an irrelevant vaccine, say rabies, prior to randomization and the immune response to the rabies vaccine at baseline is a good choice for the baseline predictor [9]. In general, however, baseline predictors could still contribute to infection risk in addition to biomarker if receiving vaccine, and their effects need to be properly accounted for valid inference. In this article, we develop methodology to estimate the mCEP curve allowing the risk model to be dependent on the baseline predictor after conditioning on the biomarker if receiving vaccine. Moreover, we develop a new procedure for making simultaneous inference of the mCEP curve over a range of biomarker values, which utilizes the perturbation resampling method to approximate the asymptotic distribution of our estimator. Making simultaneous inference about the mCEP curve is important for understanding biomarker-defined principal strata effect modification over a range of biomarker values but to our knowledge has not been addressed previously. All existing work [2, 7, 10, 11] covered inference at individual points of the mCEP curve only.

The research for this paper was motivated by two phase 3 dengue vaccine efficacy trials: the CYD14 trial conducted in Asia (Capeding et al. [12]) and the CYD15 trial conducted in Latin America (Villar et al. [13]) where the biomarker of interest is neutralizing antibody titer to the vaccine measured at 13 months post randomization. Participants were randomly assigned in 2:1 allocation to receive three injections of a live attenuated tetravalent dengue vaccine (containing one dengue strain each from the four serotypes of dengue) or placebo at months 0, 6, and 12 and were followed for 25 months post first vaccination with active surveillance for occurrence of the primary study endpoint of symptomatic virologically confirmed dengue disease (VCD). Out of the 10,275 (20,869) participants in CYD14 (CYD15), there were 117 (176) VCD cases in the vaccine group and 133 (221) VCD cases in the control group that took place after 13 months post randomization. We demonstrate application of our proposed method using data from CYD14 and CYD15 in Section 5.

The remainder of this article is organized as follows. In Section 2, we introduce the problem setting and propose an estimator for the mCEP curve that accounts for effects of the baseline predictor together with the biomarker if receiving vaccine in the risk model. In Section 3, we introduce a perturbation resampling method to approximate the distribution of our estimator for the mCEP curve and make simultaneous inference. In Section 4, we evaluate the finite-sample performance of the proposed estimator through extensive numerical studies; we also evaluate the performance of the proposed resampling procedure for constructing confidence bands and compare it with alternative bootstrap procedures. In Section 5, we present an analysis of the two phase 3 dengue vaccine efficacy trials, CYD14 and CYD15, using our proposed methods. Finally, in Section 6 we end the article by discussing some potential limitations of our methods and making suggestions for future research.

2 Methods

We consider data from a two-arm vaccine efficacy trial that randomizes n participants to either the vaccine arm or the placebo arm, with Z being the binary indicator of assignment to the vaccine arm. Let W be the baseline covariates measured in everyone such as demographics, S be the intermediate response endpoint measured at a fixed time τ post randomization, Y be the indicator of clinical endpoint occurrence (e. g. infection with a pathogen) after τ over a fixed follow-up period, and Yτ be the indicator that a subject has the clinical endpoint before τ. The variable S is only measured among those who remain free of Y through time τ and is thus undefined if Yτ=1. To define the estimand of interest, we use potential outcomes, where all post-randomization measurements are considered under either z = 0 or z = 1 for each individual. Let S(z), Yτ(z), Y(z) be the potential outcomes if the subject receives treatment z, for z = 0 or 1. If Yτ(z)=1, S(z) is undefined and we set S(z) = *.

Gilbert and Hudgens [2] defined the mCEP curve as

where

for z = 0, 1 and the function h(x, y) is a known contrast function satisfying h(x, y) = 0 if and only if x = y. A common choice of the contrast function h is the vaccine efficacy function:

which we refer to as the VE curve and it constitutes the estimand of interest in this paper. The VE curve measures the percentage of reduction in the clinical endpoint rate for the subgroup of vaccine recipients at-risk for the clinical endpoint at τ under both treatment assignments and with immune response S(1)=s1 compared to what it would have been had they been assigned to the placebo arm. Large variation in the VE curve indicates strong effect modification. The VE curve can be a useful tool for the future development of vaccines by providing a ranking of immune response biomarkers by their strength of effect modification.

However, the mCEP curves are not identifiable from the standard assumptions made in randomized vaccine efficacy trials due to missing potential outcomes. To address this problem, two augmented trial designs have been proposed by Follmann [9]. The first utilizes baseline immunogenicity predictors (BIPs) to develop an imputation model for the unobserved immune response biomarker values based on the estimable relationship between baseline covariates and biomarker values. For example, in dengue vaccine trials, study participants are vaccinated either with dengue vaccine (z = 1) or placebo (z = 0) and immune response to the dengue vaccine (S) is measured at month 13 post randomization in the vaccine group. Therefore the potential outcome S(1) is observable for the vaccine group but missing for the placebo group. If we can identify a baseline covariate(s) W that is predictive of S(1), and W is measured in subjects in both the vaccine and the placebo groups, then a model predicting S(1) from W fit from the vaccine group can be used to predict S(1) for the placebo group. The second augmented design vaccinates all or a fraction of placebo recipients who remain free of the clinical endpoint at the closeout of the trial, and the immune response biomarker Sc is measured at time τ after vaccination. These values Sc are then used to substitute for the missing biomarker values S(1) as if they originally had been assigned to the vaccine arm. For example, in dengue vaccine trials, if a CPV component were to be added, a subset of placebo patients who are free of the dengue disease primary endpoint at the end of the trial follow-up period will receive the dengue vaccine at study closeout and their immune response will be measured at month 13 post vaccination. Follmann named the second augmented design “closeout placebo vaccination” (CPV). A major advantage of including a CPV component is that it makes it possible to evaluate the risk model assumptions. We call the BIP design with no CPV component the BIP-only design and the BIP design with a nonzero CPV component the BIP+CPV design.

Frequently a two-phase sampling design is used. In the first phase, baseline covariates W and the clinical outcome data Y and Yτ are measured for everyone, and in the second phase, the potential outcome S(1) (and Sc if a CPV component is added) is measured in a subcohort of study participants sampled according to a random mechanism. Sampling of S(1) can depend on other phase-I variables such as Y. For example, in our motivating dengue studies, a two-phase case-control sampling using a BIP-only design could carry out in a way such that in the first phase, baseline covariates W such as age, gender, and country along with the clinical endpoint of dengue disease Y and Yτ are measured for everyone while in the second phase, all cases in the vaccine arm (defined by Y = 1, Yτ=0, and Z = 1) have S(1) measured and a random subset of controls in the vaccine arm (defined by Y = 0 and Z = 1) are sampled for S(1) measurement. On the other hand, a two-phase case-control sampling using a BIP+CPV design could carry out in a way such that in the first phase, W, Y and Yτ are measured for everyone while in the second phase, all cases in the vaccine arm and a random subset of controls in the vaccine arm have S(1) measured, and a random subset of placebo recipients who are uninfected of dengue at the end of the trial (defined by Y = 0 and Z = 0) are sampled for Sc measurement.

When Gilbert et al. [14] examined the power for detecting vaccine efficacy modification using a parametric estimated likelihood method, the results were seemingly counterintuitive. In some scenarios where W was strongly correlated with S(1) (and Sc if a CPV component is added), the BIP-only design was more powerful than the BIP+CPV design. Huang et al. [7] investigated this result and concluded that the decreased efficiency caused by the CPV sampling was due to the fact that the CPV component is included in the step of maximizing the likelihood but not in the step of estimating the conditional distribution of the biomarker. To address the inconsistent use of the CPV component, Huang et al. [7] proposed a pseudo-score type estimator suitable for both the BIP-only design and the BIP+CPV design to identify the parameters β in the assumed parametric risk model riskZS(1),W≡PY(Z)=1|S(1),W,Yτ(1)=Yτ(0)=0=gβ;S(1),Z,W. Furthermore, in order to estimate the mCEP curve, specifically the VE curve as a function of S(1), Huang et al. [7] considered the scenario that the risk model is independent of W after conditioning on S(1): riskZS(1),W=gβ;S(1),Z. Then it follows that the VE curve estimator can be represented as a function of the risk model directly:

In general, however, clinical risks under Z might not be independent of W after conditioning on S(1). Under these scenarios, the estimation of the VE curve requires information not only on the risk model, but also on the distribution of the baseline covariate W conditional on the immune response biomarker if receiving vaccine S(1). Next, we propose an estimator for the VE curve applicable to both the BIP-only design and the BIP+CPV design that is built upon the pseudo-score estimator by Huang et al. [7] but allows the risk model to depend on W after conditioning on S(1).

3 Simultaneous inference

Drawing simultaneous inference for the VE curve over a range of biomarker values is of interest in biomarker-defined principal strata effect modification analysis. For example, to evaluate whether conditional vaccine efficacy departs from a specific value, ve, for all S’s in the range of sl,su, the null hypothesis to be tested is H0: VE(s1)=ve for all s1∈sl,su. To evaluate whether the VE curve of one biomarker S equals the VE curve of another biomarker S^′ for values in the range of sl,su, the null hypothesis to be tested is H0: VE(S=s1)=VE(S′=s1) for all s1∈sl,su. Similarly, to evaluate whether the VE curve for trial 1 equals the VE curve for trial 2, the null hypothesis is H0: VE(s1|W1=1)=VE(s1|W1=0) for s1∈sl,su, where W1 is the indicator of trial 1 (i. e. W1=1 for trial 1 and W1=0 for trial 2). Such simultaneous inference typically involves estimation of the distribution of a process, which is often not tractable explicitly. Furthermore, explicit variance estimation might not be feasible and/or reliable. In this article, we propose a perturbation resampling method Parzen, Wei, and Ying inspired by Parzen, Wei, and Ying [18] to approximate the distribution of our estimator for VE(s1) and to draw simultaneous inference.

The construction of simultaneous confidence bands requires approximating the distribution of a Gaussian process Wˆs1≡nVEˆ(new)(s1)−VE0(s1). We propose a perturbation resampling procedure that provides a valid estimate for the distribution of Wˆs1, based on which we construct the pointwise confidence intervals and simultaneous confidence bands for the VE curve. Because VE ranges from negative infinity to 1, we perform our estimation on the log scale of relative risk (RR), where RR(s1)=1−VE(s1). To be specific, the perturbation estimation can be carried out using the following resampling procedure:

1. Generate n random realizations of ε from a known distribution with mean of 1 and variance of 1 to create E≡ϵi,i=1,2,...,n.

2. Use E to obtain the perturbed estimator βˆ(ϵ) by solving U(ϵ)(β,FN(ϵ),πˆ(ϵ))=0, where

   U(ϵ)(β,FN(ϵ),πˆ(ϵ))=∑δi=1Uβ(Yi|Si,Zi,Wi)⋅ϵi+∑δi=0∫ϵi⋅Uβ(Yi|s,Zi,Wi)P(Yi|s,Zi,Wi)Pˆ(ϵ)(δ=1|s,Wi)dFN(ϵ)(s|Wi,δ=1)∫P(Yi|s,Zi,Wi)Pˆ(ϵ)(δ=1)|s,Wi)dFN(ϵ)(s|Wi,δ=1),FN(ϵ)(s|w,δ=1)=∑δi=1I(Si≤s,Wi=w)⋅ϵi∑δi=1I(Wi=w)⋅ϵi,πˆ(ϵ)(Si,Wj)=Pˆ(ϵ)(δ=1|Si,Wj)=∑z=01∑y=01Pˆ(ϵ)(δ=1|y,z,Wj)P(y|Si,z,Wj)P(z),Pˆ(ϵ)(δ=1|y,z,Wj)=∑kI(δ=1,Yk=y,Zk=z,Wk=Wj)⋅ϵk∑kI(Yk=y,Zk=z,Wk=Wj)⋅ϵk.

3. Use E to obtain the perturbed estimator γˆ(ϵ) by solving

   Ψ(ϵ)(γ)=1n∑i=1nδiπ0iht(s1i)(gji−pji)⋅ϵi=0

   when π0 is known, or

   Ψ(ϵ)(γ)=1n∑i=1nδiπˆ(ϵ)(yi,zi,wi)ht(s1i)(gji−pji)⋅ϵi=0

   when π0 is unknown, where πˆ(ϵ)(yi,zi,wi)=Pˆ(ϵ)(δ=1|yi,zi,wi), which is defined in Step 2.

4. With βˆ(ϵ) and γˆ(ϵ), we obtain the perturbed version of logRRˆ(s1) as:

   logRRˆ(ϵ)(s1)=logrisk1(ϵ)(s1)risk0(ϵ)(s1)=log∑j=1Dgβˆ(ϵ);S=s1,Z=1,W=wj⋅μjγˆ(ϵ);s1∑j=1Dgβˆ(ϵ);S=s1,Z=0,W=wj⋅μjγˆ(ϵ);s1.

Repeat Steps 1–4 B0 times to obtain B0 realizations of logRRˆ(ϵ)(s1), denoted by logRRˆ(b)(s1),b=1,2,...,B0. The empirical distribution of WˆlogRR(s1)(b)≡nlogRRˆ(b)(s1)−logRRˆ(s1) conditional on the observed data can be used to approximate the distribution of WˆlogRR(s1)≡nlogRRˆ(s1)−logRR0(s1). In the Supplemental Material Section C, we provide theoretical justification for why the distribution of WˆlogRR(s1)(b) given the observed data Oi=(Zi,Wi,Yiτ,Yi,δi,Siδi)′, i = 1,..., n can be used to approximate the unconditional distribution of WˆlogRR(s1). With the above resampling method, one may calculate the sample standard deviation, σˆlogRR(s1) of the B0 realizations logRRˆ(b)(s1),b=1,2,...,B0. A 100(1−α)% pointwise confidence interval and simultaneous confidence band for logRR(s1),s1∈ζ may be constructed as logRRˆ(s1)±Z1−α/2σˆlogRR(s1) and logRRˆ(s1)±Q1−ασˆlogRR(s1), respectively, where Zη is the 100ηth percentile of the standard normal distribution N(0, 1) and Qη is the 100ηth percentile of sups1∈ζσˆlogRR(s1)−1WˆlogRR(s1)(b),b=1,2,...,B0. Finally, the Wald 100(1−α)% pointwise and simultaneous confidence bands for VE(s1) are obtained by transformation of the symmetric bounds from the logRR scale back to the VE scale.

4 Simulation studies

We conducted simulation studies to examine the finite-sample performance of our proposed estimator VEˆ(new)(s1) and the perturbation resampling procedure for inference. For comparison, we also studied VEˆ(A8)(s1), the estimator that makes the extra assumption (A8) that after accounting for Z and S the conditional risk does not depend on W, as well as a traditional bootstrap procedure for making pointwise and simultaneous inference. Specifically, we assessed (1) the bias of VEˆ(new)(s1) compared to VEˆ(A8)(s1), (2) the distribution of the estimated standard errors of logRRˆ(s1) obtained by the perturbation resampling procedure compared to the traditional bootstrap procedure, and (3) the empirical coverage probabilities of the resulting pointwise and simultaneous confidence bands based on perturbation compared to the bootstrap.

Simulated data followed a 2:1 vaccine:placebo randomized two-arm trial with 3000 subjects (2000 vaccine recipients and 1000 placebo recipients). The variable S was generated from a normal distribution with mean 3 and standard deviation 1. Simulated values of S less than 0 were truncated to 0. W was generated to have four categories, 1, 2, 3, 4, from the following multinomial model conditional on S: lnP(W=1|S)P(W=4|S)=−1.99+0.89S; lnP(W=2|S)P(W=4|S)=−4.80+1.84S; lnP(W=3|S)P(W=4|S)=−9.79+3.29S. The parameters in the multinomial model were chosen such that there were about equal numbers of subjects in each of the four categories and the correlation between W and S was 0.5. The risk model P(Y=1|S,Z,W) was assumed to take the form P(Y=1|S,Z,W)=Φ(β0+β1Z+β2S+β3SZ+β4W) with Φ denoting the cdf of the standard normal distribution. We chose the risk model parameters under three different scenarios: (1) the probability of Y = 1 in the placebo arm (r0) equals 0.090 and the probability of Y = 1 in the vaccine arm (r1) equals 0.042, (2) r0=0.055 and r1=0.020, (3) r0=0.0090 and r1=0.0068. We call these three scenarios the Non-Rare case, the Med-Rare case, and the Rare case, respectively, and they reflect the intention-to-treat cohort arm-specific probability of dengue disease in CYD14, CYD15, and of HIV infection in RV144, a phase 3 HIV vaccine efficacy trial conducted in Thailand (Rerks-Ngarm et al. [19]), respectively. We also studied the performance of our estimators under different values of two ratios for two-phase sampling: rv, the average ratio of sampled controls to cases in the vaccine arm; and rp, the average ratio of sampled controls in the placebo arm to cases in the vaccine arm. Under the BIP-only design, rp=0 and S was treated as missing for all placebo recipients. We considered three values for rv in our simulation: 5, 10, and All, where rv=All denotes the scenario where all vaccine recipients at-risk at τ had S measured. Under the BIP+CPV design, we considered three settings: rv=5 and rp=5; rv=10 and rp=10; and rv=All and rp=All, the last of which denotes the scenario where all vaccine recipients had S measured and all event-free placebo recipients were included in the CPV component. We show the results from the BIP+CPV design in this section. Analyses from the BIP-only design yielded similar conclusions and corresponding results have been included in the Supplemental Material Section D.

4.1 Bias of VEˆ(A8)(s1) and VEˆ(new)(s1)

For each of 1000 simulated data sets, the estimates VEˆ(A8)(s1) and VEˆ(new)(s1) were computed. Figure 1 displays the true VE curve, the average VEˆ(A8)(s1), and average VEˆ(new)(s1) over the 1000 simulations for different sampling ratios in the Non-Rare, Med-Rare, Rare case, respectively. It also reports the average sampling proportions of S for each treatment arm π‾v,π‾p, with the sampling proportions in each of the 1000 simulated data sets calculated as: πv=∑i=1nI(δi=1,Yi=0,Zi=1)∑i=1nI(Yi=0,Zi=1) and πp=∑i=1nI(δi=1,Yi=0,Zi=0)∑i=1nI(Yi=0,Zi=0). It can been seen in Figure 1 that the VEˆ(A8)(s1) estimator could significantly underestimate the true VE for small values of S. For example, the bias in VEˆ(A8)(s1) at s1=0.5 for rv=5 and rp=5 was –73 % in the Rare case, –67 % in the Med-Rare case, and –37 % in the Non-Rare case. On the other hand, the bias for VEˆ(new)(s1) is generally negligible across all scenarios. These results confirmed the superior performance of our proposed estimator over the estimator making the assumption of (A8), riskZS,W=riskZS=gβ;S,Z, when the assumption is violated.

Figure 1:

Estimated vaccine efficacy curve using our proposed method without assumption A8 and using the [7] method with assumption A8, compared to the true VE curve for checking the bias of these two estimators based on 500 simulated datasets for the Rare case where the probability of Y = 1 for Z = 0 (r0) equals 0.090 and for Z = 1 (r1) equals 0.042, the Med-Rare case where r0=0.055 and r1=0.020, and the Non-Rare case where r0=0.0090 and r1=0.0068 with a BIP+CPV design.

4.2 Standard error estimator

To examine the finite-sample performance of our proposed resampling procedure, we calculated the estimated standard errors of logRRˆ(s1) obtained by the perturbation resampling procedure. For each of the 1000 simulated datasets, B0=500 perturbations were used and E=ϵi,i=1,2,...,n were generated from the exponential distribution with rate 1. The results were not sensitive to the choice of the distribution of E and B0=500 was generally sufficient to approximate standard errors well. We also estimated standard errors of logRRˆ(s1) based on 500 nonparametric bootstrap samples. For comparison, we calculated the Monte Carlo standard errors as a benchmark for the correct standard errors. The results are summarized in Figure 2. Generally the perturbation resampling procedure yielded standard error estimates closer to the Monte Carlo standard errors than the nonparametric bootstrap procedure. When the disease endpoint is rare, perturbation-based standard errors were remarkably smaller than the bootstrap-based standard errors.

Figure 2:

Estimated standard errors of logRRˆ(s1), solid for the perturbation resampling approach and dashed for the bootstrap approach, for the Rare case, the Med-Rare case, and the Non-Rare case with a BIP+CPV design.

4.3 Coverage probabilities of pointwise and simultaneous confidence bands

We present coverage probabilities of the 95 % pointwise confidence intervals (CIs) for logRR(s1) for different s1 values in Figure 3. In all scenarios, perturbation methods yielded coverage levels closer to the nominal 95 % across most s1 values compared to the bootstrap methods. The bootstrap CIs over-covered the truth, especially in the Rare case, due to the fact that the bootstrap-based standard error estimates were remarkably large, yielding wide confidence intervals and high bootstrap CI coverage, while the perturbation-based standard error estimates continued to perform well as shown in Section 4.2. The empirical coverage levels of the 95 % simultaneous confidence bands from the perturbation/bootstrap methods are also reported in Figure 3 as “sim.cover.per”/“sim.cover.boot”. Perturbation and bootstrap yielded similar simultaneous confidence band coverage, with coverage percentage close to the nominal 95 % in the Med-Rare and the Non-Rare case. Based on these results, we conclude that the perturbation methods demonstrate a clear advantage over the the bootstrap methods, especially when the disease endpoint is rare, by producing smaller estimated standard errors and therefore narrower pointwise confidence intervals and simultaneous confidence bands while maintaining nominal coverage.

Figure 3:

Empirical coverage probabilities of 95 % pointwise confidence intervals and simultaneous confidence bands about VE(s1), for the Rare case, the Med-Rare case, and the Non-Rare case with a BIP+CPV design.

Lastly, we compare efficiency of pseudo-score type estimators for β allowing β4≠0 (denoted by βˆ(new)) to the pseudo-score type estimators that set β4=0 (denoted by βˆ(A8)) when the true risk model riskZS,W=Φ(β0+β1Z+β2S+β3SZ+0W), therefore the conditional risk does not depend on W after accounting for Z and S, and the assumption (A8) made in VEˆ(A8) is correct. The relative efficiency, defined as the ratio of sampling variance of βˆ(new) and sampling variance of βˆ(A8), was 1.25, 1.08, 1.28, and 1.07 for the Med-Rare Case with rv=10 and rp=10 for β0, β1, β2, β3, respectively. Relative efficiency for a BIP-only design is included in the Supplemental Material Section D. There is some efficiency loss by using the proposed approach allowing β4≠0 when assumption (A8) holds, but the efficiency loss is minor especially for estimating the interaction between S and Z.

In summary, our proposed procedure for estimating the VE curve had negligible bias and our proposed perturbation resampling method yielded confidence intervals with proper coverage probabilities. The perturbation-based standard error estimators in general were smaller than the bootstrap-based standard error estimators especially when the disease endpoint is rare.

5 Dengue example

We demonstrate application of our proposed method using data from CYD14 and CYD15. Based on a Cox model, estimated vaccine efficacy against VCD due to any serotype diagnosed after Month 13 and by Month 25 was 56.5 % (95 % confidence interval, 43.8 to 66.4) for CYD14 and 60.8 % (95 % confidence interval, 52.0 to 68.0) for CYD15. Neutralizing antibody titers were measured to each of the four parental dengue strains at Month 13 (time τ) in case-control samples. For illustration of our methods, we let S be the individual’s sample average response to the four serotype-specific log10-transformed antibody titers (i. e. “average titer”) and include the subjects’ age and country in the baseline covariate vector W, and Y is the indicator of the same dengue disease endpoint noted above. See Moodie et al. [20] for the reporting of the full analysis to a clinical audience.

A probit risk model was assumed conditional on Z, S, and W: riskZS,W=Φ(β0+β1Z+β2S+β3SZ+β4W). Our proposed methods were conducted on the data set of 9–16 year olds pooling across both CYD14 and CYD15 to assess how VE varied with Month 13 neutralizing antibody titers if assigned to receive the vaccine. Justification for pooling the data across the trials is provided in Moodie et al. [20], with main justification that the protocols for these two trials were essentially the same. In Figure 4 we present the point estimates and pointwise and simultaneous confidence intervals for the VE curve, showing that VE against VCD increases with average titer. Estimated VE reached 42.1 %, 85.6 %, and 96.4 %, respectively, at titer 100, 1000, and 10,000.

Furthermore, let LB(s1) denote the lower bound and UB(s1) denote the upper bound of the simultaneous band for VE(s1). Then a test of the hypothesis H0:VE(s1)=ve0 for all s1∈ζ with size α is obtained by rejecting H0 when and only when, at least for one s1, the statement

is false. It is evident that the size of the test is α, since by the definition of the simultaneous bands constructed in Section 3, 1−α is the chance for the statement (12) to be simultaneously true for all s1∈ζ. The fact that a horizonal line cannot be drawn between the 100(1−α)% simultaneous confidence bands without intersecting both the lower and the upper limits in Figure 4 indicates at least for one s1, the statement (12) is false. Therefore, we reject H0 and conclude that VE(s1) significantly varies with s1.

Figure 4:

Estimated vaccine efficacy against virologically confirmed dengue disease of any serotype through Month 25 by Month 13 average titer if receiving vaccine with 95 % pointwise confidence intervals and simultaneous confidence bands in 9–16 year olds in the two Phase 3 dengue vaccine efficacy trials combined (CYD14 and CYD15).

6 Discussion

The identification and evaluation of response markers as surrogate endpoints and as effect modifiers are important goals in many biomedical research areas. A response type biomarker shown to be a strong modifier of clinical treatment efficacy can be used to help with trial design, trial implementation, exploration of biological mechanisms of clinical treatment efficacy, and for selecting biomarker study endpoints for evaluating treatments/vaccines in new trials.

In this paper we have proposed a procedure for estimating the marginal CEP curve. Our proposed methods have the advantage of allowing clinical risks under Z being dependent on W after conditioning on S. We also developed procedures for obtaining pointwise and simultaneous confidence intervals about the marginal CEP curve via perturbation resampling. In addition, we have shown that pointwise and simultaneous interval estimates via perturbation resampling are more accurate and tighter than those obtained by the bootstrap, especially when the disease endpoint is rare. Theoretically, if sample sizes are large enough, and there are a large number of events, bootstrap and perturbation resampling procedures should both yield correct inference. However, through simulation we discovered that 500 perturbation resamples approximate standard errors better than 500 bootstrap resamples. When the number of events is small, and some events are associated with extreme S values, the bootstrap resampling procedure will sometimes delete these high influential points, and sometime repeat them. Therefore, the influence from these points will be wiped out in some bootstrap samples while compounded in others, yielding unstable estimates for different samples. On the other hand, the perturbation resampling procedure keeps the same data points but ‘perturbs’ their influence through continuously generated weights, thus yielding more stable estimates. Our simulation results also showed that the precision gain from perturbation compared to bootstrap is more dramatic when the disease endpoint is rare. Therefore, the perturbation procedure is favorable given finite sample size and/or rare events but the relative advantage of perturbation should diminish with a large number of events as the relative influence of the most influential points should be diminished.

While we have defined the baseline covariates W to be categorical such that we estimated F(S|W,δ=1) nonparametrically and employed a multinomial model for P(W|S), our proposed definition and estimation procedure can be extended to scenarios where W is continuous by adopting a nonparametric smoothing technique to estimate F(S|W,δ=1) and using a parametric, semiparametric, or nonparametric model to estimate P(W|S). Furthermore, this article assumed the baseline covariates W are phase-I variables that are available in all subjects. Extending our estimator to settings where some Ws are only available from a subset of participants is of interest for future research. The R code implementing the proposed methods is available upon request.