PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Stat Med. Author manuscript; available in PMC 2017 September 20.
Published in final edited form as:
Stat Med. 2016 September 20; 35(21): 3745–3759.
Published online 2016 March 31. doi:  10.1002/sim.6952
PMCID: PMC4965346
NIHMSID: NIHMS773751

Power/Sample Size Calculations for Assessing Correlates of Risk in Clinical Efficacy Trials

Abstract

In a randomized controlled clinical trial that assesses treatment efficacy, a common objective is to assess the association of a measured biomarker response endpoint with the primary study endpoint in the active treatment group, using a case-cohort, case-control, or two-phase sampling design. Methods for power and sample size calculations for such biomarker association analyses typically do not account for the level of treatment efficacy, precluding interpretation of the biomarker association results in terms of biomarker effect modification of treatment efficacy, with detriment that the power calculations may tacitly and inadvertently assume that the treatment harms some study participants. We develop power and sample size methods accounting for this issue, and the methods also account for inter-individual variability of the biomarker that is not biologically relevant (e.g., due to technical measurement error). We focus on a binary study endpoint and on a biomarker subject to measurement error that is normally distributed or categorical with two or three levels. We illustrate the methods with preventive HIV vaccine efficacy trials, and include an R package implementing the methods.

Keywords: Case-cohort design, Case-Control design, Immune response biomarkers, Measurement error, Principal stratification, Two-phase sampling design, Vaccine efficacy trial

1. Introduction

Commonly, clinical efficacy trials randomize study participants to receive a treatment or control preparation (e.g., placebo) at one or more visits, and follow these participants for occurrence of the primary clinical study endpoint. The primary objective assesses treatment efficacy against the clinical endpoint, and a common secondary objective assesses the association of intermediate response endpoints (e.g., biomarkers) measured after the administration of treatment with primary endpoint occurrence in the active treatment group. Applications of this secondary objective include developing prognostic biomarkers and providing information for other analysis objectives such as surrogate endpoint and mediation assessment. Typical statistical approaches for assessing such correlates of risk (CoRs) have included logistic or Cox proportional hazards regression models that account for the sampling design that was used for measuring the biomarkers (e.g., [14]).

For power calculations to detect CoRs in a cohort such as an active treatment group, many methods have been developed for case-cohort studies (e.g., [5]), case-control studies (e.g., [67]), and the generalization of case-control studies to two-phase sampling studies (e.g., [8]). However, the available approaches typically do not account for the level of clinical treatment efficacy overall and in biomarker response subgroups, precluding interpretation of the results in terms of potential correlates of efficacy/protection. We develop an approach to CoR power/sample size calculations that accounts for this issue, which is important because if the power calculations are based solely on the biomarker-outcome association in the active treatment group, then one could design a case-control study to, say, have 90% power to detect a biomarker-outcome odds ratio of 0.5, but not realize that this power is achieved under a tacit assumption that the endpoint rate is higher in the active arm than the control arm for the subgroup with lowest biomarker responses. By specifying overall treatment efficacy and biomarker-specific treatment efficacies as input parameters, our approach makes transparent in the power calculations the link between the CoR effect size in the active treatment arm and the corresponding difference in biomarker-specific treatment efficacies.

In addition, our approach accounts for the component of inter-individual variability of the biomarker that is not biologically relevant (e.g., due to technical measurement error of the device employed to measure a biological response), which is important because the degree of measurement error of the biomarker heavily influences power of the CoR analysis, such that accounting for this issue is needed to obtain accurate power calculations. In our approach the user inputs a parameter ρ defined as the estimated fraction of the biomarker’s variance that is potentially biologically relevant for protection, and displays how power and sample size requirements vary with ρ.

Our approach can be used for a general binary clinical endpoint model with case-cohort, case-control, or two-phase sampling of the biomarker, using without replacement or Bernoulli sampling. We illustrate the approach with a logistic regression model and case-control without replacement sampling. For rare event studies (e.g., with cumulative endpoint rate less than 10%), we found in simulations that the power for the logistic regression model tends to be very similar to that for a Cox regression model [9]; thus in this setting the approach may provide sufficiently accurate power results for time-to-event CoR analysis. The simplification afforded by using a binary outcome is helpful for focusing attention on the two issues listed above.

Related research has developed power calculators of testing procedures for assessing the association of a true biomarker subject to measurement error and a sub-sampling design with an outcome (e.g., [1012]). Here, we depart from this research objective by developing a power calculator of testing procedures for assessing the association of a measured/observed biomarker that has components of variability thought to be not possibly associated with the outcome. Whereas the former testing procedures incorporate bias-correction techniques, leveraging, for example, validation sets or replicate biomarker measurements, our power calculator may be used with a large number of available hypothesis testing procedures from the case-cohort/case-control/two-phase sampling statistical methods literature (going back to Horvitz and Thompson [13]), where the methods do not need bias-correction techniques. Thus, the contribution of this work is to provide more interpretable and accurate power calculations for the common scientific endeavor to understand power for detecting the association of a measured/observed biomarker with the outcome. Moreover, previous work has developed power calculation formulas for associating a measured biomarker subject to measurement error with a dichotomous outcome; for example [1416] considered a normally distributed biomarker following a classical measurement error model, with application to logistic regression correlates analysis.

While the newly proposed power calculator applies for general randomized controlled two-group clinical trials, for definiteness we focus on preventive vaccine efficacy trials, which randomize study participants to receive a candidate vaccine or placebo at one or more visits, and follow these participants for occurrence of clinically significant infection with the pathogen under study [17]. The primary objective assesses vaccine efficacy (V E) defined as the multiplicative reduction (vaccine versus placebo) in the rate of the primary endpoint, and a secondary objective assesses the association of immune response biomarkers measured shortly after vaccination with the primary endpoint. For settings where some trial participants were previously infected with the pathogen (e.g., influenza) this analysis is done for each of the vaccine and placebo groups or pooling over the groups, and for settings where trial participants have not been previously infected with the pathogen (e.g., HIV), such that the immune response biomarker does not vary in the placebo group [18], this analysis is done either pooling over the vaccine and placebo groups or in the vaccine group only. In the vaccine field such analyses have been named CoR analyses (e.g., [1920]) and for definiteness we focus on assessing a CoR in the vaccine group. The approach is illustrated with power calculations for the RV144 HIV vaccine efficacy trial after the primary analysis was conducted, and with sample size calculations for the prospective design of a sequel HIV vaccine efficacy trial being planned by the HIV Vaccine Trials Network.

Section 2 describes the study set-up, parameters of interest, and identifiability assumptions. Section 3 describes the power and sample size calculation approach. Section 4 illustrates the power/sample size calculator with the two examples, and Section 5 concludes with discussion. Supporting Materials Appendix A discusses how to unbiasedly characterize the biomarker distribution accounting for the sampling design, Supporting Materials Appendix B provides selected mathematical details of the power calculation methods, and Supporting Materials Appendix C addresses the important topic of how to estimate the noise level of the biomarker. Supporting Materials Appendix D presents supplementary figures for the two illustrations and Supporting Materials Appendix E summarizes how to use the R package.

2. Study Set-Up, Parameters of Interest, Identifiability Assumptions

2.1. Randomized Clinical Trial for Assessing Vaccine Efficacy

We consider a double-blind clinical trial that randomizes participants to vaccine or placebo, with Z the indicator of assignment to vaccine and W baseline covariates. Let S be the immune response biomarker measured at a fixed time τ post-randomization, which we assume to be continuous or trichotomous, with the case of dichotomous S covered as a special case. Participants are followed for occurrence of the primary clinical study endpoint, clinically significant infection with the pathogen, with followup through time τmax, with T the time from randomization until the study endpoint and Y [equivalent] I[T ≤ τmax] the binary outcome of interest. Let Yτ [equivalent] I[T ≤ τ] and Vτ be the indicator that a subject attends the visit at τ. Fitting to the motivating application, we focus on settings where it is only interesting to study the association of S with Y for subjects who did not experience the event before the biomarker is measured. Therefore, subjects with (1 − Yτ)Vτ = 1 are the subgroup observed to be at-risk at τ who could potentially have S measured for the association study.

Because S is expensive to measure, a case-cohort, case-control, or two-phase sampling design is often used; let R be the indicator that S is measured. Let Δ be the indicator that Y is observed, i.e., Δ = 0 if the subject drops out before time τmax and before experiencing the event, and Δ = 1 otherwise. Let L [equivalent] (R(z), R(z)S(z), Yτ(z), Vτ(z), Δ(z), Δ(z)Y (z)) be the potential outcomes if assigned treatment z, for z = 0, 1, where S(z) is defined if and only if Yτ (z) = 0, such that S(z) = * if Yτ (z) = 1. (Note that Yτ (z) = 1 and Vτ (z) = 0 each imply R(z) = 0.) The observed data for a subject are O [equivalent] (Z,W,R,RS, Yτ, Vτ, Δ, ΔY). The CoR power calculations are based on the N vaccine recipients observed to be at-risk at τ (those with Z(1 − Yτ)Vτ = 1), and test for whether P(Y = 1|S = s1, Z = 1, Yτ = 0) varies in s1. To understand our approach, it is critical to note that the CoR power calculations do not need the potential outcomes formulation, as they are based solely on the observable random variables O. The potential outcomes are used to define biomarker-specific vaccine efficacy and hence provide a way to relate CoR effect sizes to vaccine efficacy effect sizes.

To facilitate building this relationship, we assume the vaccine has no effect on the study endpoint before the biomarker sampling time τ: P(Yτ (1) = Yτ (0)) = 1; this assumption will be more credible and less influential for τ near baseline. This assumption is useful by ensuring that the biomarker-specific vaccine efficacy parameters measure causal effects of vaccination, and for equating the CoR parameter P(Y = 1|S = s1, Z = 1, Yτ = 0) to P(Y = 1|S = s1, Z = 1, Yτ (1) = Yτ (0) = 0), which links the CoR and V E parameter types (as described below). Henceforth all unconditional and conditional probabilities of Y (z) = 1 tacitly condition on Yτ (1) = Yτ (0) = 0.

2.2. Vaccine Efficacy Parameters: Trichotomous Biomarker

We suppose that each of the N vaccine recipients is in one of three latent/unknown baseline subgroups X, the “lower protected” (X = 0), the “medium protected” (X = 1), or the “higher protected” (X = 2). Define the x-specific outcome risks as

riskzlat(x)P(Y(z)=1|X=x) for x=0,1,2 and z=0,1,
(1)

such that the vaccine efficacy for latent subgroup x is VExlat1RRxlat with RRxlatrisk1lat(x)/risk0lat(x), for x = 0, 1, 2.

Define PxlatP(X=x) for x = 0, 1, 2, and define the marginal risks riskz [equivalent] P(Y (z) = 1) for z = 0, 1. Then the overall vaccine efficacy V E equals

VE=1RR=1risk1risk0=1P0latrisk1lat(0)+P1latrisk1lat(1)+P2latrisk1lat(2)P0latrisk0lat(0)+P1latrisk0lat(1)+P2latrisk0lat(2).
(2)

We also define risks and vaccine efficacies for subgroups defined by S(1) or by (X, S(1)):

riskz(s1)P(Y(z)=1|S(1)=s1),riskzrisk(x,s1)P(Y(z)=1|X=x,S(1)=s1)
(3)

for x = 0, 1, 2, s1 = 0, 1, 2 and z = 0, 1, and

VE(s1)1RR(s1)=1risk1(s1)/risk0(s1)

VElat(x,s1)1RRlat(x,s1)=1risk1lat(x,s1)/risk0lat(x,s1).

The observed biomarker response s1 = 0 represents a “low” response in some fashion and s1 = 2 a higher response, with s1 = 1 an intermediate response. For example, s1 = 0 could be a negative/non-response and s1 = 2 a response above a pre-specified putative correlate of protection threshold. If S were measured without error, then X = S such that V E(s1) = V Elat(x, s1) and the latent variable formulation would not be needed; we use it to allow measurement error to create differences in V E(s1) versus V Elat(x, s1), with greater differences for noisier biomarkers (developed next).

2.3. Accounting for Measurement Error in the Biomarker

To incorporate assay noise into the power/sample size calculations, we define protection-related sensitivity/specificity and false positive/negative parameters as:

SensP(S(1)=2|X=2),SpecP(S(1)=0|X=0),
(4)

FP0P(S(1)=2|X=0),FN2P(S(1)=0|X=2),
(5)

FP1P(S(1)=2|X=1),FN1P(S(1)=0|X=1).
(6)

The probability an observed at-risk vaccine recipient has a low or high response, P0 [equivalent] P(S = 0|Z(1 − Yτ)Vτ = 1) or P2 [equivalent] P(S = 2|Z(1 − Yτ)Vτ = 1), equals

P0=Spec*P0lat+FN1*P1lat+FN2*P2lat,
(7)

P2=Sens*P2lat+FP1*P1lat+FP0*P0lat.
(8)

We consider two approaches to the trichotomous biomarker power calculations. Approach 1 takes as inputs (Sens, Spec, FP0, FN2, FP1, FN1), whereas Approach 2 uses an additive measurement error model for a normally distributed continuous-readout biomarker S*, and defines the values of S by S = 0 if S* ≤ θ0, S = 2 if S* > θ2, and S = 1 otherwise, with θ0 and θ2 two user-specified constants with θ0 < θ2. In particular, for Approach 2 we consider a normally distributed latent ‘true’ biomarker X* and link S* to X* by an additive classical measurement error model

S*=X*+e,X*~N(0,σtr2),e~N(0,σe2),
(9)

with X* independent of e, implying S*~N(0,σobs2) with σobs2=σtr2+σe2. Here ρ1σe2/σobs2 is the fraction of the variability of S* that is potentially biologically relevant for protection, and is specified to reflect the quality of the biomarker. The ‘true’ trichotomous biomarker X is defined by two percentiles of X* that are determined mathematically by model (9) and the two percentiles θ0 and θ2 (see Supporting Materials Appendix B). Figure 1 illustrates the set-up for Approach 2.

Figure 1
Division of participants into three latent subgroups with low, medium, and high levels of vaccine efficacy VE0lat,VE2lat and VE2lat, with prevalences P0lat,P1lat and P2lat, respectively. Under Approach 2 of Step 7 in Section 3.1 the latent normal measurement ...

The above set-up handles a dichotomous biomarker as a special case, by setting P1lat=P1=0, in which case only the Sens and Spec parameters are needed for the calculations [because FN2 = 1 − Sens and FP0 = 1 − Spec; see equations (4)–(8)]. The R code handles the dichotomous biomarker as a special case.

2.4. Vaccine Efficacy Parameters and Model: Continuous Biomarker

The formulation for a continuous biomarker is similar, where now the latent subgroups are defined by the true unobservable biomarker X* in model (9) above. Now

VElat(x*)1risk1lat(x*)/risk0lat(x*),VE(s1)1risk1(s1)/risk0(s1),

with riskzlat(x*)P(Y(z)=1|X*(1)=x*) and riskz(s1) [equivalent] P(Y (z) = 1|S*(1) = s1) for x* and s1 varying over the continuous support of X*(1) and S*(1), respectively.

For the power calculations we specify a fraction PlowestVElat of subjects with the lowest X*(1) values ≤ ν to all have the same specified lowest level of vaccine efficacy V Elowest:

VElowestVElat(X*(1)ν)=1risk1lat(ν)/risk0lat(ν).
(10)

For example, V Elowest may be set to 0 and PlowestVElat defined as the fraction of subjects without a positive vaccine-induced immune response. The constant ν is determined by PlowestVElat, V Elowest, and the measurement error model (9): ν=ρσobsΦ1(PlowestVElat), where Φ−1(·) is the inverse of the standard normal cdf.

For x* ≤ ν, risk1lat(x*) is modeled as a constant following (10),

risk1lat(x*)=(1VElowest)risk0lat(ν)for x*ν,
(11)

and, for x* > ν, risk1lat(x*) is modeled via a logistic regression model

logit(risk1lat(x*))=αlat+βlatx*for x*>ν.
(12)

Using model (11)(12) that specifies a lowest value of vaccine efficacy is useful because the alternative simpler model that would specify (12) for all x would force V E(x) to be negative for the lowest values of x. In many applications this is undesirable as enhanced risk of disease caused by vaccination may be considered unlikely and the most relevant power calculations would dissallow this possibility. (Albeit the power calculator works for V Elowest specified negative.)

Model (11)(12) combined with (9) and the assumption risk0lat(x)=risk0 as stated in Section 2.7 below implies that

VE=1[PlowestVElatrisk1lat(ν)+νlogit1(αlat+βlatx*)ϕ(x*/ρσobs)dx*]/risk0,
(13)

where ϕ(·) is the standard normal pdf. This formula will be used later for implementing the power calculations.

2.5. Correlate of Risk (CoR) Hypotheses and Estimands of Interest

We address the scientific objective to assesses a CoR among vaccine recipients. For trichotomous S this entails testing the following null versus alternative hypotheses

H0:risk1(s1=2)=risk1(s1=1)=risk1(s1=0) vs.
(14)

H1:risk1(s1=2)risk1(s1=1)risk1(s1=0)
(15)

with ‘<’ for at least one of the two inequalities in H1. For continuous S* this tests

H0:risk1(s1) is consant in s1vs.H1:risk1(s1)risk1(s1) for all s1<s1
(16)

with ‘<’ for some s1<s1. While for data analysis 2-sided tests would typically be used, the power calculations are clearer to interpret by testing for the 1-sided alternative H1 of lower clinical risk in vaccine recipients with increasing s1.

2.6. Methods of Analysis with Bernoulli and Without Replacement Sampling

Two main approaches to selecting the subset of subjects for whom to measure the biomarkers are:

  • Prospective case-cohort [1]: Select a simple or stratified random sample from all randomized vaccine recipients, and augment the sample with all study endpoint cases that were not randomly sampled; and
  • Retrospective case-control or 2-phase sampling [4]: Conditional on final case status and possibly a discrete stratification covariate measured in all subjects, select a fixed number of vaccine recipients from each case status × covariate stratum.

Our sample size calculations consider both approaches. The first approach has advantages including that the randomly sampled subjects can be used for unbiased assessment of the distribution of the biomarker in the study population, absolute risks can be assessed in biomarker subgroups, and the association of biomarkers with multiple study endpoints can be straightforwardly assessed. The second approach does not facilitate the latter two goals, and some re-weighting is required to use the sampled subjects for unbiased assessment of the distribution of the biomarker in the study population (see Supporting Materials Appendix A). An advantage of the second approach is that waiting until the primary analysis is completed before selecting the controls allows accounting for the vaccine efficacy results for optimizing the biomarker sampling design. This affords opportunities to improve efficiency of the analysis [4].

2.7. Identifiability Assumptions

In addition to assuming iid random variables (Li, Xi*, Xi) and (Oi, Xi*, Xi) for i = 1, …, N, we assume the standard set of assumptions that have been used in correlates of risk and protection studies: SUTVA, ignorable treatment assignment (Z [perpendicular] L|W), equal early clinical risk (P(Yτ (1) = Yτ (0)) = 1), and random censoring (Y (z) [perpendicular] Δ(z) for z = 0, 1). We also assume S is missing at random (MAR): R depends only on the observed data O. To the extent the investigator controls the biomarker sampling design MAR is guaranteed to hold, although it could be in question due to happenstance missingness caused by not attending the visit at τ. Moreover, we focus on the scenario that after accounting for the latent category (and any baseline covariates W included in the CoR analysis) the measured biomarker in vaccine recipients does not affect risk, i.e., risk1lat(x*,s1)P(Y(1)=1|X*(1)=x*,S*(1)=s1)=risk1lat(x*) for all s1 and x*, and similarly for risk as a function of trichotomous X and S.

We develop the power calculations for the relatively simple scenario of homogeneous risk in the placebo group, where risk0lat(x*,s1)=risk0(s1)=risk0 for all s1 and x* and similarly for risk as a function of trichotomous X and S. In general, risk0(x*, s1) and risk0(s1) are not identifiable (because S(1) is a counterfactual random variable for subjects assigned Z = 0), and power calculations could be conducted under many scenarios for these functions. However, the special case is very helpful for power calculations because risk0 can be specified based on the observed or projected incidence in the trial. Because the CoR data analysis itself would control for known baseline prognostic factors W, the scenario in which the power calculations are accurate is risk0lat(x*,s1)=risk0(s1)=risk0 after conditioning on W.

2.8. CoR Effect Sizes RRt and RRc as a Function of Vaccine Efficacies

From (14), (15), and (16), analysis of the vaccine group data provides inference on the relative risks RRt [equivalent] risk1(2)/risk1(0) for a trichotomous biomarker and RRc [equivalent] risk1(s1)/risk1(s1 − 1) for a continuous biomarker. We refer to RRt and RRc as the user-specified “CoR effect sizes” of the power calculations. From the assumptions of Section 2.7, RRt and RRc are identified from the observed data measured from the subset of vaccine recipients with R = 1, because they imply risk1(s1) = P(Y = 1|Z = 1, R = 1, s1) [21]. Therefore under the assumptions the power calculations for testing H0 can be based on the set of vaccine recipients with S (or S*) measured at τ.

For a trichotomous biomarker, straightforward calculation shows that RRt is linked to the latent V E parameters via the equation

RRt=RR0lat*FP0+RR1lat*FP1+RR2lat*SensRR0lat*Spec+RR1lat*FN1+RR2lat*FN2.
(17)

This formula makes the estimable RRt interpretable in terms of a gradient in vaccine efficacies, where RRt=RR2lat/RR0lat for a noise-free biomarker with 1 − Sens = 1 − Spec = FP0 = FP1 = FN2 = FN1 = 0 (illustrated in Figure 5 below). Otherwise, under H1, RRt is closer to 1.0 than RR2lat/RR0lat.

Figure 5
For the RV144 scenario with overall V E = 0.26, CoR effect size RRt = risk1(2) /risk1(0) versus the ratio RR2lat/RR0lat measuring relative vaccine efficacy for the higher protected and lower protected latent subgroups, for four scenarios of ρ. ...

For a continuous biomarker S* following model (9), RRc is linked to the latent vaccine efficacy parameters via an equation that depends on s1. Because RRc depends on s1, it is not particularly useful to index power calculations by RRc. Instead, we interpret RRc as the effect size for a noise-free biomarker (ρ = 1). Under the logistic model (12), RRc is the relative risk per standard deviation increase in X* in the region above ν, where we use the approximation of a relative risk by an odds ratio.

3. Power and Sample Size Calculations

3.1. Without Replacement Sampling

Of the N vaccine recipients observed to be at-risk at τ, let ncases (ncontrols) be the number of observed cases (controls) from whom S (or S*) is measured, where cases have ΔY = 1 and controls have Δ(1 − Y) = 1. If the power calculations are done at the design stage then N, ncases, and ncontrols are projected numbers.

For a trichotomous S, the algorithm for the power calculations is as follows:

  1. Specify the overall vaccine efficacy between τ and τmax, V E. For example, if the power calculations are done before the trial then V E may be set to the protocol-specified design alternative and if afterwards to the estimated V E.
  2. Specify risk0, either based on the protocol-specified placebo group endpoint rate projection or as an estimate after the trial (e.g., estimated as n1/(n1 +n2) where n1 is the number of observed placebo cases between τ and τmax and n2 is the number who reached time τmax free of the outcome). The estimator of risk0 should be defensibly unbiased accounting for participant dropout.
  3. Select a grid of vaccine efficacies for the lower protected subgroup, VE0lat. One useful grid ranges from overall V E specified in Step 1 (the null hypothesis) to 0 (the maximal alternative hypothesis not allowing harm by vaccination).
  4. Select a grid of vaccine efficacies for the medium protected subgroup, VE1latVE0lat; one useful default choice that we use in the illustration is VE1lat=VE.
  5. Specify P0lat (which in many applications may be specified as the rate of “negative” response) and specify P2lat. These values determine P1lat=1P0latP2lat. Based on equation (2) assuming risk0(2) = risk0(0) = risk0,
    VE=P0latVE0lat+P1latVE1lat+P2latVE2lat.
    This formula determines VEhiVElat to be
    VE2lat=[VE*(P0lat+P2lat)P0lat*VE0lat]/P2lat.
    If one varies VE0lat from V E to 0 as suggested in Step 3, then by the above equation VE2lat varies from V E to VE*(P0lat+P2lat)/P2lat.
  6. Specify values P0 and P2, which determines P1 = 1 − P0P2. Thus in total the user specifies (VE, risk0, VE0lat,VE1lat,P0lat,P2lat, P0, P2) Typically it is of interest to set P0=P0lat and P2=P2lat as a key scenario for computing power, as in the illustration.
  7. Approach 1: Following equation (7), specify two of the three parameters (Spec, FN2, FN1), which determines the remaining parameter. Similarly, following equation (8), specify two of the three parameters (Sens, FP0, FP1), which determines the remaining parameter. One choice specifies Sens and Spec and sets FN2 = FP0 = 0, because a biomarker of reasonable quality should have FN2 and FP0 close to zero.
    Note that Approach 1 provides a completely general approach to studying a trichotomous or dichotomous biomarker as a CoR, without making any use of the normal measurement error model (9).
    Approach 2: Specify σobs2 and ρ, which together with the values (P0, P2, P0lat,P2lat) fixed in Step 6 determine θ0 and θ2 and determine (Sens, Spec, FP0, FN2, FP1, FN1) (see Table 1 in the illustration). Typically σobs2 is set to 1.0 as S* can always be scaled to have unit variance.
    With Approach 2, it is helpful to plot the CoR effect size RRt versus the latent protection effect size RR2lat/RR0lat via formula (17) for different values of ρ (illustrated in Figure 5).
  8. Simulate a large number of data sets under the above true parameter values. The procedure fixes the total number of cases ncases and the total number of controls ncontrols = Kncases for an input counting number K, such that S (or S*) are the randomly generated variables.
  9. For each simulated data set compute the Wald test statistic for H0 from a logistic regression model (with S the covariate of interest) that uses inverse probability weighting to account for the marker sub-sampling design. Any method for valid testing of the biomarker-outcome association using MAR case-cohort/case-control/ two-phase sampling may be used in this step.
  10. Compute the power as the fraction of the simulated data sets where the Wald test statistic yields 1-sided p ≤ α/2 for specified α with direction favoring H1.
  11. Given the specified V E and risk0, repeat the power calculations varying ρ from 1.0 (noise-free biomarker) to a value less than one reflecting a worst case for noise level of the biomarker. Figure 4 in the next section illustrates these calculations.
    Figure 4
    Vaccine group CoR power calculations for a trichotomous biomarker S for the completed RV144 HIV vaccine efficacy trial with ncases = 41 and ncontrols = 205 (1-sided α = 0.05/2-level Wald test), for four scenarios of ρ. The x-axis is risk ...
  12. Repeat the power calculations for different controls:case ratios K, and, if done before the trial, possibly for different values of V E and risk0.

For Step 7 Approach 2, under the assumptions of Section 2 and specified values for σobs2, ρ, and (P0lat,P2lat), the remaining inputs Sens, Spec, FP0, FN2, FP1, FN1, (θ0, θ2) are determined by solving the equations (4)(6) and (7)(8); the R code does this using stochastic integration (Supporting Materials Appendix B).

Step 8 proceeds as follows. First, for each of the N vaccine recipients, determine the numbers that are in the three latent subgroups as Nx=Pxlat*N rounded to the nearest integers, for x = 0, 1, 2. Second, determine the latent class membership of each of the ncases cases by a realization of a trinomial random variable with success probabilities (P(X = 0|Y = 1, Yτ = 0, Z = 1), P(X = 1|Y = 1, Yτ = 0, Z = 1), P(X = 2|Y = 1, Yτ = 0, Z = 1)), where P(X = x|Y = 1, Yτ = 0, Z = 1) is expressed in terms of the Platx and risk1(x) via Bayes rule. This determines the number of cases ncases(x) in each category x = 0, 1, 2 satisfying ncases=x=02ncases(x). Third, within each subgroup x, simulate Si for the entire set of N subjects as a trinomial random variable. For x = 0 the response probabilities are (Spec, 1 − FP0Spec, FP0); for x = 1 the response probabilities are (FN1, 1 − FP1FN1, FP1); and for x = 2 the response probabilities are (FN2, 1−SensFN2, Sens). This determines the number of controls ncontrols(x) in each category x = 0, 1, 2 by subtracting off ncases(x), satisfying the constraint ncases=x=02ncases(x) where ncontrols is fixed at K * ncases. Fourth, finalize the analysis data set by specifying Ri = 1 or Ri = 0 for each of the N subjects.

For Step 9, we use the tps(·) function in the R package osDesign that implements the 2-phase logistic regression method of Breslow and Holubkov [22], entering S as an ordered score variable with levels S = 0, 1, 2 and conducting a one degree of freedom Wald test. Alternatively a generalized two degree of freedom Wald test could be used. In addition, alternative analysis methods could be used that leverage correlations between the biomarker and auxiliary covariates measured in everyone, potentially increasing power [4]. However, it will often be advantageous to base the power calculations on the simpler method both for the utility of having conservative power calculations and because the strength of correlation of the auxiliaries must be fairly high to yield a material power gain (often not available in practice).

For a continuous normally distributed biomarker S* scaled to have mean 0 and σobs2=1, the same simulated data sets (using Approach 2 in Step 7) can be used for the power calculations, with process as follows.

  1. Steps 1 and 2 as for the trichotomous biomarker case.
  2. Fix PlowestVElat, V Elowest, and ρ [set membership] (0, 1]. Under the logistic regression model (12) for risk1lat(x), V Elat(x) = 1 − logit−1lat + βlatx)/risk0 for x > ν and V Elat(x) = V Elowest for x ≤ ν. The fixed values (V E, risk0, PlowestVElat, V Elowest) mathematically determine αlat and βlat by the above equation and equation (12) (solutions in Supporting Materials Appendix B).
  3. Given the specified (V E, risk0, PlowestVElat, V Elowest), plot the VE curve V Elat(x) verus x given a true CoR effect size RRc = explat). V Elat(x) is calculated using model (11)(12), where given fixed βlat, βlat,αlat=logit(risk1lat(ν))βlatν. Repeat the analysis for multiple values of RRc (see the illustration in Figure 3).
    Figure 3
    For the RV144 scenario with overall V E = 0.26, the vaccine efficacy curve V E(x) for the true biomarker X* = x for six different scenarios of the true CoR relative risk effect size RRc in the vaccine group and values of V Elowest (all curves assume ρ ...
  4. Data sets are simulated by first specifying risk1lat(x*) following models (11)(12) on a fine grid of x* values. Second the true latent biomarkers Xi* for the ncases cases are sampled from P(X* = x*|Y = 1, Yτ = 0, Z = 1) calculated using Bayes rule. Similarly the Xi* values for the ncontrols = K * ncases controls are sampled from P(X* = x*|Y = 0, Yτ = 0, Z = 1).
  5. Steps 9–12 as above, except that Step 9 now uses S* as the covariate of interest in the logistic regression model, again implemented with tps(·).

3.2. Bernoulli Sampling

Under Bernoulli sampling, of the N vaccine recipients observed (or projected) to be at-risk at τ, ncases (ncontrols) is the expected number of observed cases (controls) from whom S and S* are measured, i.e., ncases and ncontrols are random. For a trichotomous biomarker, the power analysis proceeds as described in Section 3.1, except Step 8 uses Bernoulli sampling (classic case-cohort sampling [1]). In particular, for each of the N vaccine recipients, determine the case status Y conditional on X* = x* as a realization of a Bernoulli random variable with success probability risk1lat(x*). For a continuous biomarker, Step 8 is altered by determining the case status Y conditional on X* = x* as a realization of a Bernoulli random variable with success probability risk1lat(x*).

4. Illustration of the Power Calculations

We first illustrate the correlates power calculations for the RV144 preventive HIV vaccine efficacy trial of a candidate vaccine versus placebo that was conducted in the general population in Thailand [23]. For this example the power calculations are conducted after the trial was completed (e.g., [2425]). The RV144 trial randomized 8198 (8197) HIV uninfected individuals to receive vaccine (placebo) and followed them for the primary endpoint of HIV infection over 42 months. Subjects received immunizations at Week 0, 4, 8, 24, and immune response biomarkers measured at τ = Month 6 (Week 26 visit) were assessed in vaccine recipients as CoRs of HIV infection by τmax = 42 months. Relevant for the CoR power calculations, estimated overall V E to prevent infections after τ through τmax was 0.26. Of the N = 7703 vaccine recipients observed to be at risk at Month 6, biomarkers were measured in the 41 subjects who were observed to subsequently experience the HIV infection endpoint, and in a frequency matched 5:1 controls:cases allocation random sample of 205 observed controls (i.e., without replacement two-phase sampling). Based on these data several papers have reported significant continuous, trichotomous, and dichotomous CoRs in the vaccine group, with initial paper Haynes et al. [25]. These analyses were done using two-phase logistic regression [22] and two-phase Cox regression [26], which gave almost identical answers.

For the power analysis with a continuous biomarker following model (12) we assume PlowestVElat=0.40, such that the 40% of vaccine recipients with lowest X* responses had vaccine efficacy V Elowest. We varied V Elowest from 0 to the overall V E estimate of 0.26. We estimated risk1 as n1/(n1 + n2) where n1 is the number of vaccine recipients observed to be at-risk at τ = 6 months who were diagnosed with HIV infection by the end of follow-up τmax = 42 months (n1 = 41) and n2 is the number of vaccine recipients observed to be at-risk at τ who completed follow-up HIV negative (n2 = 7662). Then we estimated risk0 as risk^1/(10.26)=0.0072.

Figure 2 shows the power curves for ρ = 1, 0.9, 0.7, 0.5. As expected power decreases with the degree of noise. The interpretation of the plot may be aided by annotating it with results from previous efficacy trials that identified CoRs. In particular, suppose a previous trial reported an estimated RR^ per sd increment in observed S*. Under the measurement error model (9), for ρ = 1 this equates to RR^ per sd increment in X*, and for fixed ρ < 1 this equates to RR^1/ρ per sd increment in X*. Typically there will be uncertainty as to the level of ρ in the previous trial, such that affixing the estimated relative risk per sd increment in X* to each curve provides a scenario analysis of the power available to detect the previously identified CoR under a spectrum of noise levels. In Figure 2 we use the gp70-V1V2 binding antibody correlate observed in RV144 [25], which had RR^=0.57 per sd increment in S*. If this biomarker is assumed to have no measurement error (ρ = 1), power is 0.19, whereas under substantial measurement error (ρ = 0.7), power drops to 0.13. In additional simulations, the power curves are higher if overall V E is higher (not shown).

Figure 2
Power to detect a normally distributed biomarker S* as an inverse correlate of risk (CoR) of HIV infection in vaccine recipients for the completed RV144 HIV vaccine efficacy trial with ncases = 41 and ncontrols = 205 (1-sided α = 0.05/2-level ...

To help interpret the power results in Figure 2, Figure 3 shows the V E(x) curves for six different scenarios of the true CoR relative risk effect size RRc (ρ = 1) and values of V Elowest for the RV144 scenario with estimated overall V E of 0.26. The null hypothesis RRc = 1 corresponds to a flat curve V E(x) = V E, and increasing departures from the null hypothesis H0 correspond to increasingly variable and steep VE curves. This figure shows that for the scenario risk0(s1, x) = risk0 and no measurement error, an association of the biomarker with infection risk in the vaccine group (a CoR) is equivalent to an association of the biomarker with V E. For interpreting Figure 2, if we focus on the ρ = 0.9 curve with effect size RRc = 0.53 and V Elowest = 0.04 (green solid curve), V E varies substantially in X* but power is low, only about 0.14.

Figure 4 shows the power curves (top panels) based on the same simulated data sets following the recipe given in Section 3.1 (using Approach 2 in Step 7), for a trichotomous biomarker with P0=P2=P0lat=P2lat set to 0.1, 0.2, 0.3, or 0.4 with RR1lat=RRoverall and RR2lat tied to RR0lat through the relationship expressed in Step 5 of Section 3.1. The results show that power majorly increases with P0 = P2, which is intuitively expected given that greater sample sizes at the poles of lowest and highest VE should yield the greatest power.

To help interpret the power results of Figure 4, Figure 5 shows the relationship between the CoR effect size RRt and the relative risk ratio RR2lat/RR0lat for the four values of ρ, with Table 1 showing how ρ maps to Sens, Spec, FP0, FN2, FP1, FN1 for each set of input parameters used in Figure 4. Figure 5 shows that for a noise-free biomarker with ρ = 1, RRt=RR2lat/RR0lat such that a CoR in the vaccine group is equivalent to the relative vaccine efficacy parameter, whereas for imperfectly measured biomarkers with ρ < 1, RRt>RR2lat/RR0lat such that the CoR effect size is closer to the null than the relative vaccine efficacy parameter. We illustrate a co-interpretation of Figures 4 and and55 for the ρ = 0.9 marker in Figure 5 and P(S = 0) = 0.4 (bottom-right panels). There is about 25% power to detect a CoR with effect size RRt = 0.60 (Figure 4), which corresponds to 25% power to detect RR2lat/RR0lat0.50 (Figure 5). Supporting Materials Figure 1 shows an ROC curve (sensitivity versus one minus specificity) as P2lat=P2=1P0=1P0lat ranges from 0.10 to 0.90. Our overall conclusion for this example is as follows: Because estimated overall V E was low (at 0.26), the assumption of V E ≥ 0 for all biomarker subgroups constrains the possible CoR effect sizes to a limited range hence yielding low power of the CoR analysis; in contrast if V E were allowed to be negative for some subgroups then power would be greater.

Table 1
Mapping of the measurement error model (9) with σobs2=1 indexed by ρ to Sens, Spec, FP0, FN2, FP1, FN1

Our second example considers calculations being used to plan the sample size of a Phase 3 HIV vaccine efficacy trial under design by the HIV Vaccine Trials Network. This trial randomizes HIV negative individuals to vaccine or placebo in a 1:1 allocation and follows subjects for HIV infection during a τmax = 36 month follow-up period. We assume 4% annual HIV incidence in the placebo group and 5% annual dropout incidence, as well as overall V E = 0.50. The immune response biomarkers to assess as CoRs are measured at month τ = 6.5. All vaccine group subjects diagnosed with HIV between month 6.5 and 36 have biomarkers measured, as do a random sample of HIV uninfected controls with controls:cases ratio 1:1, 3:1, 5:1, or 10:1. Figure 6 shows the trichotomous biomarker power curves versus the number of infections in the vaccine group (and the total sample size observed to be at risk at τ) to detect a CoR effect size of risk1lat(0)=0.065,risk1lat(1)=0.043,risk1lat(2)=0.022, for ρ fixed at value 0.9 that may be a realistic scenario for a biomarker assessed as a CoR. (Under the constant placebo risk assumption these calculations assume VE0lat=0.25,VE1lat=0.50,VE2lat=0.75) The calculations are for the scenarios P0=P0lat=P2lat=P2 ranging from 0.10 to 0.50. The results show that power sharply increases with the prevalences P0lat=P2lat and increases with the controls:cases ratio, with only incremental gain moving from 5:1 to 10:1. Based on this analysis, to achieve 90% power to detect a CoR with P0 = P2 = 0.30, one choice would be the 5:1 allocation design, requiring 2800 total vaccine recipients observed to be at-risk at 6.5 months.

Figure 6
Vaccine group CoR power curves versus total sample size (Number of participants observed to be at risk at time τ) for a trichotomous biomarker S measured at time τ to plan a 2-arm HIV vaccine efficacy trial with equal allocation randomization ...

5. Discussion

We developed an approach to power and sample size calculations for a typical “correlates of risk” (CoR) data analysis in a randomized controlled clinical efficacy trial for testing an association of an observed biomarker measured in a sub-sample (via a case-cohort, case-control, or two-phase sampling design) of the active treatment group with a clinical endpoint. The contribution of this work is to integrate into the calculations two issues – the level of treatment efficacy across biomarker subgroups and the fraction ρ of the variability of the biomarker that is potentially biologically relevant for protection (ρ1σe2/σobs2). The first issue is important because, if ignored, a statistician may design the sample size of a CoR study not realizing the tacit assumptions being made about treatment efficacy. A particularly egregious mistake would be powering a study to detect a CoR with no recognition that achieving the desired power requires that treatment efficacy be negative for some biomarker subgroups, rendering the CoR study underpowered if treatment efficacy is never negative. Our approach provides a way to explicitly explore the relationship of the CoR effect size with treatment efficacy, including a way to specify the lowest treatment efficacy at a fixed value such as zero. The second issue is important because the degree of measurement error ρ heavily influences power [1416], such that accounting for ρ is needed for accurate power calculations, and may be useful for screening out biomarkers for which the CoR study would be underpowered given an unacceptably low value of ρ.

For the continuous biomarker calculations and for the Approach 2 trichotomous biomarker calculations, we have assumed a classical additive normal measurement error model for the observed continuous biomarker S*, the veracity of which should be tested. In general, in the planning of biomarker CoR studies it is important to conduct biomarker assay laboratory validation studies to estimate ρ; we discuss approaches to this in Web Appendix C.

Our power calculator applies for a univariate biomarker, yet studying the association of multiple biomarkers with outcome is an important application. The calculator for a trichotomous biomarker may be useful for trials that collect possibly high-dimensional multivariate biomarkers, and for which unsupervised clustering based on the biomarkers yields a cluster of “putatively not protected” subjects and a cluster of “putatively protected” subjects. In this scenario, the power calculator may be applied with all other subjects constituting the third cluster. In addition, the calculator for a normally distributed biomarker may be used for studying power to detect a linear combination of multiple biomarkers as a CoR.

Our CoR power and sample size calculations are for the scenario that the biomarker is not associated with the clinical endpoint in the control group after accounting for baseline covariates W that would be controlled for in the CoR data analysis. This assumption is not needed for the CoR calculations because they use data from the active treatment group only. However, this assumption is used as a way to interpret the CoR power calculations in terms of biomarker-specific treatment efficacy, providing a mapping from the CoR calculations (in terms of risk gradients in the active treatment group) to gradients in treatment efficacy. Additional calculations may be conducted under alternative scenarios, where the approach here could be extended to allow functions risk0(x, s1) other than risk0(x, s1) = risk0. While the main application of the methods is more interpretable and accurate CoR power and sample size calculations, a second application is power and sample size calculations for assessing modification of treatment efficacy by the biomarker, i.e., assessing the vaccine efficacy curve V E(s1) directly, which is conducted in the principal stratification framework [18, 27]. Supporting Materials Figures 2 and 3 show such power curves for our two illustrative examples.

Supplementary Material

Supp Info

Acknowledgments

Research reported in this publication was supported by the National Institute Of Allergy And Infectious Diseases (NIAID) of the National Institutes of Health (NIH) under Award Numbers R37AI054165 and UM1AI068635. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. The authors thank the participants, investigators, and sponsors of the RV144 trial, including the U.S. Military HIV Research Program (MHRP); U.S. Army Medical Research and Materiel Command; NIAID; U.S. and Thai Components, Armed Forces Research Institute of Medical Science; Ministry of Public Health, Thailand; Mahidol University; SanofiPasteur; and Global Solutions for Infectious Diseases.

Footnotes

Web-based Supporting Materials

Title: Appendices Appendix A describes a technique for unbiased characterization of the biomarker distribution accounting for the biomarker sampling design. Appendix B provides selected mathematical details of the power calculation methods. Appendix C discusses how to estimate the noise level of the biomarker under study. Appendix D presents supplementary figures for the two illustrations. Appendix E summarizes how to use the R package implementing the methods.

References

1. Prentice R. A case-cohort design for epidemiologic cohort studies and disease prevention trials. Biometrika. 1986;73:1–11.
2. Barlow W. Robust variance estimation for the case-cohort design. Biometrics. 1994;50:1064–1072. [PubMed]
3. Kulich M, Lin D. Improving efficiency of relative-risk estimation in case cohort studies. Journal of the American Statistical Association. 2004;99:832–844.
4. Breslow N, Lumley T, Ballantyne C, Chambless L, Kulich M. Using the whole cohort in the analysis of case-cohort data. American Journal of Epidemiology. 2009;169:1398–1405. [PMC free article] [PubMed]
5. Cai J, Zeng D. Sample size/power calculation for case-cohort studies. Biometrics. 2004;60:1015–1024. [PubMed]
6. Dupont WD, Plummer WD., Jr Power and sample size calculations: a review and computer program. Controlled Clinical Trials. 1990;11:116–128. [PubMed]
7. Garc00ED;a-Closas M, Lubin JH. Power and sample size calculations in case control studies of gene-environment interactions: comments on different approaches. American Journal of Epidemiology. 1999;149:689–692. [PubMed]
8. Haneuse S, Saegusa T, Lumley T. osDesign: an R package for the analysis, evaluation, and design of two-phase and case-control studies. Journal of Statistical Software. 2011;43:11. [PMC free article] [PubMed]
9. Gilbert PB, Grove D, Gabriel E, Huang Y, Gray G, Hammer S, Buchbinder S, Kublin J, Corey L, Self SG. A sequential Phase 2b trial design for evaluating vaccine efficacy and immune correlates for multiple HIV vaccine regimens. Statistical Communications in Infectious Diseases. 2011;3(Issue 1) Article 4. [PMC free article] [PubMed]
10. Holcroft CA, Spiegelman D. Design of validation studies for estimating the odds ratio of exposure-disease relationships when exposure is misclassified. Biometrics. 1999;55:1193–1201. [PubMed]
11. Tosteson TD, Buzas JS, Demidenko E, Karagas M. Power and sample size calculations for generalized regression models and covariate measurement error. Statistics in Medicine. 2003;22:1069–1082. [PubMed]
12. Chen LM, Ibrahim JG, Chu H. Sample size and power determination in joint modeling of longitudinal and survival data. Statistics in Medicine. 2011;30:2295–2309. [PMC free article] [PubMed]
13. Horvitz DG, Thompson DJ. A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association. 1952;47(260):663–685.
14. Armstrong BG, Whittemore AS, Howe GR. Analysis of case-control data with covariate measurement error: application to diet and colon cancer. Statistics in Medicine. 1989;8:1151–1163. [PubMed]
15. Rosner B, Willett WC, Spiegelman D. Correction of logistic regression relative risk estimates and confidence intervals for systematic within person measure error. Statistics in Medicine. 1989;8:1051–1070. [PubMed]
16. McKeown-Eyssen GE, Tibshirani R. Implications of measurement error in exposure for the sample sizes of case-control studies. American Journal of Epidemiology. 1994;139(4):415–421. [PubMed]
17. Halloran ME, Longini IM, Struchiner CJ. Design and Analysis of Vaccine Studies. New York: Springer; 2010.
18. Gilbert PB, Hudgens MG. Evaluating candidate principal surrogate endpoints. Biometrics. 2008;64:1146–1154. [PMC free article] [PubMed]
19. Qin L, Gilbert PB, Corey L, McElrath MJ, Self SG. A framework for assessing an immunological correlate of protection in vaccine trials. The Journal of Infectious Diseases. 2007;196:1304–1312. [PubMed]
20. Plotkin S, Gilbert PB. Nomenclature for immune correlates of protection after vaccination. Clinical Infectious Diseases. 2012;54:1615–1617. [PMC free article] [PubMed]
21. Gabriel E, Gilbert PB. Evaluating principle surrogate endpoints with time-to-event data accounting for time-varying treatment efficacy. Biostatistics. 2014;15:251–265. [PMC free article] [PubMed]
22. Breslow NE, Holubkov R. Maximum likelihood estimation of logistic regression parameters under two-phase, outcome-dependent sampling. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 1997;59:447–461.
23. Rerks-Ngarm S, Pitisuttithum P, Nitayaphan S, Kim JH. for the MOPH-TAVEG Investigators. Vaccination with ALVAC and AIDSVAX to prevent HIV-1 infection in Thailand. New England Journal of Medicine. 2009;361:2209–2220. [PubMed]
24. Rolland M, Gilbert PB. Evaluating immune correlates in HIV Type 1 vaccine efficacy trials: what RV144 may provide. AIDS Research and Human Retroviruses. 2012;28:400–404. [PMC free article] [PubMed]
25. Haynes BF, Gilbert PB, McElrath MJ, Zolla-Pazner S, Tomaras GD, Alam SM, Evans DT, Montefiori DC, Karnasuta C, Sutthent R, Liao HX, DeVico AL, Lewis GK, Williams C, Pinter A, Fong Y, Janes H, deCamp A, Huang Y, Rao M, Billings E, Karasavvas N, Robb ML, Ngauy V, de Souza MS, Paris R, Ferrari G, Bailer RT, Soderberg KA, Andrews C, Berman PW, Frahm N, De Rosa SC, Alpert MD, Yates NL, Shen X, Koup RA, Pitisuttithum P, Kaewkungwal J, Nitayaphan S, Rerks-Ngarm S, Michael NL, Kim JH. Immune correlates analysis of the ALVAC-AIDSVAX HIV-1 vaccine efficacy trial. New England Journal of Medicine. 2012;366:1275–1286. [PMC free article] [PubMed]
26. Borgan L, Langholz B, Samuelson S, Pogoda J. Exposure stratified case-cohort designs. Lifetime Data Analysis. 2000;6:39–58. [PubMed]
27. Frangakis C, Rubin D. Principal stratification in causal inference. Biometrics. 2002;58:21–29. [PMC free article] [PubMed]