Introduction to R Package CoRpower

The CoRpower package performs power calculations for correlates of risk, as described in Gilbert, Janes, and Huang (2016). Power/Sample Size Calculations for Assessing Correlates of Risk in Clinical Efficacy Trials. The power calculations accommodate three types of biomarkers:

  • trichotomous
  • dichotomous
  • continuous,

as well as two types of sampling design:

  • without replacement sampling for retrospective case-control designs
  • Bernoulli sampling for prospective case-cohort designs.

The vignette aims to illustrate distinct features of the functions in the package (with some mathematical background) by walking through a number of power calculation scenarios for different biomarker types, sampling designs, and input parameters.

The functions included in this package are:

  • computeN()
  • computePower()
  • plotPowerTri()
  • plotPowerCont()
  • plotRRgradVE()
  • plotVElatCont()

Set-up and notation | Without replacement sampling

  • Assume a randomized vaccine vs. placebo/control vaccine efficacy trial

  • Participants are followed for the first occurrence of the primary clinical study endpoint, with follow-up through time τmax

  • T is the time from randomization (or enrollment) to the primary endpoint

  • Y = I(T ≤ τmax) is the binary outcome of interest

  • Δ is the indicator that Y is observed; i.e., Δ = 0 if the participants drops out before τmax and before experiencing the primary endpoint and Δ = 1 otherwise

  • N1 is the number of vaccine recipients observed or expected to be at risk at τ (typically, τ is the biomarker sampling timepoint)

  • ncases, 1 (ncontrols, 1) is the number of observed or expected cases (controls) in vaccine recipients between τ and τmax, where cases have ΔY = 1 and controls have Δ(1 − Y) = 1

    • Note that both cases and controls have Δ = 1
    • ncases, 1 + ncontrols, 1 ≤ N1

  • ncases, 1S (ncontrols, 1S) is the number of observed cases (controls) in vaccine recipients between τ and τmax with measured biomarker response S(1) (or S*(1))

  • If calculations done at design stage, then N1, ncases, 1, ncontrols, 1, ncases, 1S, and ncontrols, 1S are expected counts

Algorithm for trichotomous biomarker S(1) | Without replacement sampling

  1. Specify true overall VE between τ and τmax
    • Protocol-specified design alternative or $\widehat{VE}$
  2. Specify risk0 = P(Y = 1 ∣ Z = 0, Yτ = 0) where Yτ = I(T ≤ τ)
    • Protocol-specified placebo-group endpoint rate or $\widehat{risk}_0$
  3. Select a grid of VE0lat values
    • E.g., ranging from VE (H0) to 0 (maximal H1 not allowing harm by vaccination)
  4. Select a grid of VE1lat ≥ VE0lat values
    • E.g., VE1lat = VE
  5. Specify P0lat and P2lat, then P1lat = 1 − P0lat − P2lat
    • Assuming risk0lat(x) = risk0, VE = VE0latP0lat + VE1latP1lat + VE2latP2lat yields VE2lat
    • If VE0lat varies from VE to 0, then VE2lat varies from VE to VE (P0lat + P2lat)/P2lat
  6. Specify P0 and P2, then P1 = 1 − P0 − P2
    • E.g., P0 = P0lat and P2 = P2lat
  7. Approach 1: Specify two of the three (Sens, FP0, FP1) and two of the three (Spec, FN2, FN1)
    • E.g., specify Sens and Spec and FP0 = FN2 = 0
    Approach 2: Specify σobs2 and ρ = σtr2/σobs2 and determine (Sens, Spec, FP0, FP1, FN1, FN2) (see below)

    • Typically, σobs2 = 1

    • Calculation of (Sens, Spec, FP0, FP1, FN1, FN2)

      1. Assuming the classical measurement error model, where X* ∼ N(0, ρσobs2), solve P0lat = P(X* ≤ θ0)  and  P2lat = P(X* > θ2) for θ0 and θ2

      2. Generate B realizations of X* and S* = X* + e, where e ∼ N(0, (1 − ρ)σobs2), and X* independent of e

        • B = 20, 000 by default

      3. Using θ0 and θ2 from the first step, define $$ \begin{align} Spec(\phi_0) &= P(S^{\ast} \leq \phi_0 \mid X^{\ast} \leq \theta_0)\\ FN^1(\phi_0) &= P(S^{\ast} \leq \phi_0 \mid X^{\ast} \in (\theta_0,\theta_2])\\ FN^2(\phi_0) &= P(S^{\ast} \leq \phi_0 \mid X^{\ast} > \theta_2)\\ Sens(\phi_2) &= P(S^{\ast} > \phi_2 \mid X^{\ast} > \theta_2)\\ FP^1(\phi_2) &= P(S^{\ast} > \phi_2 \mid X^{\ast} \in (\theta_0,\theta_2])\\ FP^0(\phi_2) &= P(S^{\ast} > \phi_2 \mid X^{\ast} \leq \theta_0) \end{align} $$

        Estimate Spec(ϕ0) by $$\widehat{Spec}(\phi_0) = \frac{\#\{S^{\ast}_b \leq \phi_0, X^{\ast}_b \leq \theta_0\}}{\#\{X^{\ast}_b \leq \theta_0\}}\,$$ etc.

      4. Find ϕ0 = ϕ0* and ϕ2 = ϕ2* that numerically solve $$ \begin{align} P_0 &= \widehat{Spec}(\phi_0)P^{lat}_0 + \widehat{FN}^1(\phi_0)P^{lat}_1 + \widehat{FN}^2(\phi_0)P^{lat}_2\\ P_2 &= \widehat{Sens}(\phi_2)P^{lat}_2 + \widehat{FP}^1(\phi_2)P^{lat}_1 + \widehat{FP}^0(\phi_2)P^{lat}_0 \end{align} $$ and compute $$ Spec = \widehat{Spec}(\phi^{\ast}_0),\; Sens = \widehat{Sens}(\phi^{\ast}_2),\; \textrm{etc.} $$

    • In Approach 2, plot $$ RR_t \quad \textrm{vs.} \quad \frac{RR^{lat}_2}{RR^{lat}_0}, $$ where RRt is the CoR effect size defined as $$ RR_t := \frac{risk_1(2)}{risk_1(0)} = \frac{\sum_{x=0}^2 RR^{lat}_x P(X=x|S(1)=2)}{\sum_{x=0}^2 RR^{lat}_x P(X=x|S(1)=0)} $$

    • If ρ < 1, then RRt is closer to 1 than RR2lat/RR0lat

      • Note that, under the assumption of homogeneous risk in the placebo group, i.e., risk0lat(x) = risk0, x = 0, 1, 2, the relative risk ratio RR2lat/RR0lat = risk1lat(2)/risk1lat(0)

      • Consequently, risk1(2)/risk1(0) > risk1lat(2)/risk1lat(0) because, if ρ < 1, then risk1(2) > risk1lat(2) and risk1(0) < risk1lat(0)

  8. Simulate M data sets under the true parameter values:
      Full data
    1. Nx = (ncontrols, 1 + ncases, 1)Pxlat
    2. (ncases, 1(0), ncases, 1(1), ncases, 1(2)) ∼ Mult(ncases, 1, (p0, p1, p2)), where pk = P(X = k|Y = 1, Yτ = 0, Z = 1)
    3. For each of the Nx subjects, generate Si(1) by a draw from Mult(1, (p0, p1, p2)), where pk = P(S(1) = k|X = x)
      4. Observed data
    4. Sample without replacement ncases, 1S and ncontrols, 1S = Kncases, 1S controls with measured S(1) (R = 1), i.e., the control:case ratio is not fixed within subgroup X = x
  9. For each observed data set, compute the 1-sided one-degree-of-freedom Wald test statistic for H0 ⇔ {H̃0 : βS(1) ≥ 0} from IPW logistic regression model that accounts for biomarker sampling design (function tps in R package osDesign)
    • Alternatively, use the generalized two-degree-of-freedom Wald test
  10. Compute power as proportion of data sets with 1-sided Wald test p ≤ α/2 for specified α
  11. Repeat power calculation varying control:case ratio, ncases, 1, ncases, 1S, (P0lat, P2lat, P0, P2), (Sens, Spec), ρ

Illustration: hypothetical randomized placebo-controlled VE trial

Trial design

  • Nrand = 4, 100 participants randomized to each of the vaccine and placebo group and followed for τmax= 24 months
  • Samples for measurement of S(1) at τ= 3.5 months stored in all vaccine recipients
  • Cumulative endpoint rate between τ and τmax in placebo group  = 3.4% ( = risk0)
  • VEτ − τmax = 75%,   VE0 − τ = VEτ − τmax/2
  • Cumulative dropout rate between 0 and τmax in both groups  = 10%

Illustration: calculation of input parameters with computeN()

Assumptions



  1. Failure time T and censoring time C are independent
  2. T ∣ Z = 0 ∼ Exp(λT) and C ∣ Z = 0 ∼ Exp(λC)
  3. $RR_{\tau-\tau_{\mathrm{max}}} := \frac{P(T \leq \tau_{\mathrm{max}} \mid T>\tau, Z=1)}{P(T \leq \tau_{\mathrm{max}} \mid T>\tau, Z=0)} = \frac{P(T \leq t\mid T>\tau, Z=1)}{P(T \leq t\mid T>\tau, Z=0)}$ for all t ∈ (τ, τmax]

Number of vaccine recipients observed to be at risk at τ

$$ \begin{align} N_1 &= N_{rand}\, P(T>\tau, C>\tau \mid Z=1)\\ &= N_{rand}\, P(T>\tau \mid Z=1)\, P(C>\tau \mid Z=1)\\ &= N_{rand}\, \{1 - RR_{0-\tau}\, P(T \leq \tau \mid Z=0)\}\, P(C>\tau \mid Z=1)\\ &\approx 4,023 \end{align} $$

Number of observed cases in vaccine recipients between τ and τmax

$$ \begin{align} n_{cases,1} &= N_1\, P(T\leq \tau_{\mathrm{max}}, T\leq C \mid T>\tau, C>\tau, Z=1)\\ &= N_1\, P(T\leq \min(C,\tau_{\mathrm{max}}) \mid T>\tau, C>\tau, Z=1)\\ &= N_1\, \frac{\int_{\tau}^{\infty}P(\tau < T \leq \min(c,\tau_{\mathrm{max}})\mid Z=1)f_C(c)\mathop{\mathrm{d}}\!c}{P(T>\tau, C>\tau \mid Z=1)}\\ &= N_1\, \frac{\bigg\{\int_{\tau}^{\tau_{\mathrm{max}}}P(\tau < T \leq c\mid Z=1)f_C(c)\mathop{\mathrm{d}}\!c + P(\tau < T \leq \tau_{\mathrm{max}}\mid Z=1) P(C>\tau_{\mathrm{max}})\bigg\}}{P(T>\tau, C>\tau \mid Z=1)}\\ &\approx 32 \end{align} $$

Number of observed controls in vaccine recipients between τ and τmax

$$ \begin{align} n_{controls,1} &= N_1 \, P(T > \tau_{\mathrm{max}}, C > \tau_{\mathrm{max}} \mid T>\tau, C>\tau, Z=1)\\ &\approx 3,654 \end{align} $$

Number of observed cases (controls) in vaccine recipients between τ and τmax with measured S(1)

$$ \begin{align} n^S_{cases,1} &= n_{cases,1}\\ n^S_{controls,1} &= K n^S_{cases,1} \end{align} $$

Compute N1, ncases, 1, ncontrols, 1, and ncases, 1S with computeN()

library(CoRpower)
computeN(Nrand = 4100,          # participants randomized to vaccine arm
         tau = 3.5,             # biomarker sampling timepoint
         taumax = 24,           # end of follow-up
         VEtauToTaumax = 0.75,  # VE between 'tau' and 'taumax'
         VE0toTau = 0.75/2,     # VE between 0 and 'tau'
         risk0 = 0.034,         # placebo-group endpoint risk between 'tau' and 'taumax'
         dropoutRisk = 0.1,     # dropout risk between 0 and 'taumax'
         propCasesWithS = 1)    # proportion of observed cases with measured S(1)
## $N
## [1] 4023
## 
## $nCases
## [1] 33
## 
## $nControls
## [1] 3653
## 
## $nCasesWithS
## [1] 33

Illustration: CoRpower for trichotomous  S(1) | Without replacement sampling

Approach 1 (Sens, Spec, FP0, FN2 specified)

  • Scenario 1: vary control:case ratio
  • Scenario 2: vary Sens, Spec
  • Scenario 3: vary P0lat, P2lat, P0, P2

Approach 2 (σobs2 and ρ specified)

  • Scenario 4: vary ρ
  • Scenario 5: vary P0lat, P2lat, P0, P2
  • Scenario 6: vary ncases, 1

Scenario 1: vary control:case ratio (Approach 1) | Trichotomous S(1), without replacement sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,                 
                    nControlsTx = 3654,            
                    nCasesTxWithS = 32,           
                    controlCaseRatio = c(5, 3, 1), # n^S_controls : n^S_cases ratio         
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.6,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.6,                     # probability of high biomarker response
                    sens = 0.8, spec = 0.8, FP0 = 0, FN2 = 0,
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default

Plot power curves with plotPowerTri()

Basic plotting functions are included in the package to aid with visualizing results. plotPowerTri plots the power curve against the CoR relative risk, RRt, for trichotomous or binary biomarkers.

Output from computePower() can be saved as an object and assigned to the outComputePower input parameter.

plotPowerTri(outComputePower = pwr,  # 'computePower' output list of lists
             legendText = paste0("Control:Case = ", c("5:1", "3:1", "1:1")))

Alternatively, output from computePower() can be saved in RData files. In this case, the outComputePower input parameter should be the name(s) of the output file(s), and the outDir input parameter should be the name(s) of the file location(s). For more information, visit the plotPowerTri() help page.

computePower(..., saveDir = "myDir", saveFile = c("myFile1.RData", "myFile2.RData", "myFile3.RData"))
plotPowerTri(outComputePower = c("myFile1.RData", "myFile2.RData", "myFile3.RData"), # 'computePower' output files
             outDir = rep("~/myDir", 3),                           # path to each myFilex.RData
             legendText = paste0("Control:Case = ", c("5:1", "3:1", "1:1")))

Scenario 2: vary  Sens and  Spec (Approach 1) | Trichotomous  S(1), without replacement sampling

pwr <- computePower(nCasesTx = 32,                 
                    nControlsTx = 3654,            
                    nCasesTxWithS = 32,           
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.6,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.6,                     # probability of high biomarker response
                    sens = c(1, 0.9, 0.8, 0.7), spec = c(1, 0.9, 0.8, 0.7), 
                    FP0 = c(0, 0, 0, 0), FN2 = c(0, 0, 0, 0),
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default
plotPowerTri(outComputePower = pwr,    
             legendText = paste0("Sens = Spec = ", c(1, 0.9, 0.8, 0.7)))

Scenario 3: vary  P0lat, P2lat, P0, P2 (Approach 1) | Trichotomous  S(1), without replacement sampling

pwr <- computePower(nCasesTx = 32,                 
                    nControlsTx = 3654,            
                    nCasesTxWithS = 32,           
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = c(0.05, 0.1, 0.15, 0.2),          
                    Plat2 = c(0.15, 0.3, 0.45, 0.6),          
                    P0 = c(0.05, 0.1, 0.15, 0.2),            
                    P2 = c(0.15, 0.3, 0.45, 0.6),          
                    sens = 0.8, spec = 0.8, FP0 = 0, FN2 = 0,
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default
plotPowerTri(outComputePower = pwr, 
             legendText = c("Plat0=0.05, Plat2=0.15", 
                            "Plat0=0.1, Plat2=0.3", 
                            "Plat0=0.15, Plat2=0.45", 
                            "Plat0=0.2, Plat2=0.6"))

Scenario 4: vary  ρ (Approach 2) | Trichotomous  S(1), without replacement sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,                 
                    nControlsTx = 3654,            
                    nCasesTxWithS = 32,           
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax 
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.6,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.6,                     # probability of high biomarker response
                    sigma2obs = 1,                # variance of observed biomarker S(1)
                    rho = c(1, 0.9, 0.7, 0.5),    # protection-relevant fraction of variance of S(1)
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default

Plot power curves with plotPowerTri()

plotPowerTri(outComputePower = pwr,   
             legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))

Plot RRt vs. RR2lat/RR0lat with plotRRgradVE()

plotRRgradVE() plots the ratio of relative risks for the higher and lower latent subgroups RR2lat/RR0lat against the CoR relative risk effect size RRt = risk1(2)/risk1(0).

Output from computePower() can be saved as an object and assigned to the outComputePower input parameter.

plotRRgradVE(outComputePower = pwr,  # 'computePower' output list of lists
             legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))

Alternatively, output from computePower() can be saved in RData files. In this case, the outComputePower input parameter should be the name(s) of the output file(s), and the outDir input parameter should be the name(s) of the file location(s). For more information, visit the plotRRgradVE() help page.

computePower(..., saveDir = "myDir", saveFile = "myFile.RData")
plotRRgradVE(outComputePower = paste0("myFile_rho_", c(1, 0.9, 0.7, 0.5), ".RData"),    # files with 'computePower' output
             outDir = "~/myDir",            # path to myFile.RData
             legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))

Plot ROC curves with plotROCcurveTri()

plotROCcurveTri() plots the receiver operating characteristic (ROC) curve displaying sensitivity and specificity for a range of values for P2, P0, Plat2, and rho. For more information, visit the plotROCcurveTri() help page.

plotROCcurveTri(Plat0 = 0.2, 
                Plat2 = c(0.2, 0.3, 0.4, 0.5), 
                P0 = seq(0.90, 0.10, len=25), 
                P2 = seq(0.10, 0.90, len=25), 
                rho = c(1, 0.9, 0.7, 0.5))

Scenario 5: vary  P0lat, P0, P2lat, P2 (Approach 2) | Trichotomous  S(1), without replacement sampling

pwr <- computePower(nCasesTx = 32,                 
                    nControlsTx = 3654,            
                    nCasesTxWithS = 32,           
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax 
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = c(0.05, 0.1, 0.15, 0.2),          
                    Plat2 = c(0.15, 0.3, 0.45, 0.6),          
                    P0 = c(0.05, 0.1, 0.15, 0.2),            
                    P2 = c(0.15, 0.3, 0.45, 0.6), 
                    sigma2obs = 1,                # variance of observed biomarker S(1)
                    rho = 0.9,                    # protection-relevant fraction of variance of S(1)
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default
plotPowerTri(outComputePower = pwr,  
             legendText = c("Plat0=0.05, Plat2=0.15", 
                            "Plat0=0.1, Plat2=0.3", 
                            "Plat0=0.15, Plat2=0.45", 
                            "Plat0=0.2, Plat2=0.6"))

Scenario 6: vary  ncases, 1 (Approach 2) | Trichotomous  S(1), without replacement sampling

pwr <- computePower(nCasesTx = c(25, 32, 35, 40),             
                    nControlsTx = c(3661, 3654, 3651, 3646),  
                    nCasesTxWithS = c(25, 32, 35, 40),       
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk fom tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.6,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.6,                     # probability of high biomarker response
                    sigma2obs = 1,                # variance of observed biomarker S(1)
                    rho = 0.9,                    # protection-relevant fraction of variance of S(1)
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default
plotPowerTri(outComputePower = pwr,   
             legendText = paste0("nCasesTx = ", c(25, 32, 35, 40)))

Illustration: CoRpower for binary  S(1) | Without replacement sampling

Achieved by selecting P0lat, P2lat, P0, P2 such that $$ \begin{align} P_0^{lat} + P_2^{lat} &= 1\\ P_0 + P_2 &= 1 \end{align} $$

Approach 2 (σobs2 and ρ specified)

  • Scenario 7: vary ncases, 1

Scenario 7: vary  ncases, 1 (Approach 2) | Binary  S(1), without replacement sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = c(25, 32, 35, 40),             
                    nControlsTx = c(3661, 3654, 3651, 3646),  
                    nCasesTxWithS = c(25, 32, 35, 40),       
                    controlCaseRatio = 5,         # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.8,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.8,                     # probability of high biomarker response
                    sigma2obs = 1,                # variance of observed biomarker S(1)
                    rho = 0.9,                    # protection-relevant fraction of variance of S(1)
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "binary")          # "continuous" by default

Plot power curves with plotPowerTri()

plotPowerTri(outComputePower = pwr,   
             legendText = paste0("nCasesTx = ", c(25, 32, 35, 40)))

Algorithm for continuous biomarker  S*(1) | Without replacement sampling

  1. Specify overall VE between τ and τmax
    • Protocol-specified design alternative or $\widehat{VE}$
  2. Specify risk0
    • Protocol-specified placebo-group endpoint rate or $\widehat{risk}_0$
  3. Specify PlowestVElat, ρ, and a grid of VElowest values (e.g., ranging from VE to 0)
    • Fixed (VE, risk0, PlowestVElat, VElowest, ρ) and $$ \begin{align} risk^{lat}_1(x^{\ast}) &= (1 - VE_{lowest}) risk_0,\quad x^{\ast} \leq \nu\\ \mathop{\mathrm{logit}}\{risk^{lat}_1(x^{\ast})\} &= \alpha^{lat}+\beta^{lat}x^{\ast},\quad x^{\ast} \geq \nu\\ VE &= 1 - \frac{\int risk^{lat}_1(x^{\ast})\phi(x^{\ast}/\sqrt{\rho} \sigma_{obs})\mathop{\mathrm{d}}\!x^{\ast}}{risk_0} \end{align} $$ yield αlat and βlat
    • Plot VEx*lat vs. x* and calculate the pertaining CoR effect size exp (βlat) for each level of VElowest
  4. Simulate M data sets under the true parameter values:
      Full data
    1. Sample X* for ncases, 1 cases from fX*(x*|Y = 1, Yτ = 0, Z = 1) using Bayes rule
    2. Sample X* for ncontrols, 1 controls from fX*(x*|Y = 0, Yτ = 0, Z = 1) using Bayes rule
      - How? Use a fine grid of x̃* values and then
        sample(x̃*, prob=fX*(x̃*|Y = ⋅, Yτ = 0, Z = 1), replace=TRUE)
    3. Sample S*(1) following S*(1) = X* + e
      4. Observed data
    4. Sample without replacement ncases, 1S and ncontrols, 1S = Kncases, 1S controls with measured S*(1) (R = 1)
  5. For each observed data set, compute the 1-sided one-degree-of-freedom Wald test statistic for H0 ⇔ {H̃0 : βS*(1) ≥ 0} from IPW logistic regression model that accounts for biomarker sampling design (function tps in R package osDesign)
  6. Compute power as proportion of data sets with 1-sided Wald test p ≤ α/2 for specified α
  7. Repeat power calculation varying control:case ratio, ncases, 1, ncases, 1S, PlowestVElat, ρ

Illustration: CoRpower for continuous  S*(1) | Without replacement sampling

  • Scenario 8: vary ρ
  • Scenario 9: vary PlowestVElat

Scenario 8: vary  ρ | Continuous  S*(1), without replacement sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,             
                    nControlsTx = 3654,         
                    nCasesTxWithS = 32,         
                    controlCaseRatio = 5,        # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,            # overall VE
                    risk0 = 0.034,               # placebo-group endpoint risk from tau - taumax
                    PlatVElowest = 0.2,          # prevalence of VE_lowest
                    VElowest = seq(0, VEoverall, len=100), # lowest VE for true biomarker X*<=nu
                    sigma2obs = 1,               # variance of observed biomarker S
                    rho = c(1, 0.9, 0.7, 0.5)    # protection-relevant fraction of variance of S
                    M = 1000,                    # number of simulated clinical trials
                    alpha = 0.05,                # two-sided Wald test Type 1 error rate
                    biomType = "continuous")     # "continuous" by default

Plot power curves with plotPowerCont()

plotPowerCont() plots the power curve against the CoR relative risk, RRc, for continuous biomarkers.

Output from computePower() can be saved as an object and assigned to the outComputePower input parameter. In this scenario, since computePower() was run multiple times to vary the controls:cases ratio, these multiple output objects can be read into the function as a list.

plotPowerCont(outComputePower = pwr,          # output list of lists from 'computePower'
              legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))

Alternatively, output from computePower() can be saved in RData files. In this case, the outComputePower input parameter should be the name(s) of the output file(s), and the outDir input parameter should be the name(s) of the file location(s). For more information, visit the plotPowerCont() help page.

computePower(..., saveDir = "myDir", saveFile = "myFile.RData")
plotPowerCont(outComputePower = paste0("myFile_rho_", c(1, 0.9, 0.7, 0.5), ".RData"),     # files with 'computePower' output
              outDir = "~/myDir",             # path to myFile.RData
              legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))

Plot VEx*lat curves with plotVElatCont()

plotVElatCont() plots the vaccine efficacy (VE) curve for the true biomarker X=x for eight different values of the true CoR relative risk, , in vaccine recipients and the lowest vaccine efficacy level for the true biomarker, .

outComputePower contains output from a single run of computePower() with no varying argument (i.e., no vectorized input parameters other than VE0lat, VE1lat, and ). This output can be in the form of an assigned object, which is a list of lists of length 1, or a character string specifying the file containing the output. Note that this is unlike plotPowerTri() and plotPowerCont(), which can take in output from computePower() in the form of a list of lists of length greater than 1 or a character vector. For more information, visit the plotVElatCont() help page.

Using the function when computePower() output is saved as list object pwr:

plotVElatCont(outComputePower = pwr)

Using the function when computePower() output is saved in a file with name “myFile” and location “~/myDir”:

computePower(..., saveDir = "myDir", saveFile = "myFile.RData")
plotVElatCont(outComputePower = "myFile.RData",   
              outDir = "~/myDir")

Scenario 9: vary  PlowestVElat | Continuous  S*(1), without replacement sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,               
                    nControlsTx = 3654,         
                    nCasesTxWithS = 32,         
                    controlCaseRatio = 5,        # n^S_controls : n^S_cases ratio
                    VEoverall = 0.75,            # overall VE
                    risk0 = 0.034,               # placebo-group endpoint risk from tau - taumax
                    PlatVElowest = c(0.05, 0.1, 0.15, 0.2),         
                    VElowest = seq(0, VEoverall, len=100), # lowest VE for true biomarker X*<=nu
                    sigma2obs = 1,               # variance of observed biomarker S(1)
                    rho = 0.9                    # protection-relevant fraction of variance of S(1)
                    M = 1000,                    # number of simulated clinical trials
                    alpha = 0.05,                # two-sided Wald test Type 1 error rate
                    biomType = "continuous")     # "continuous" by default

Plot power curves with plotPowerCont()

plotPowerCont(outComputePower = pwr,          # output list of lists from 'computePower'
              legendText = paste0("PlatVElowest = ", c(0.05, 0.1, 0.15, 0.2)))

Bernoulli / case-cohort sampling of  S(1) (or  S*(1))

  • Bernoulli sample at baseline with sampling probability p
  • S(1) (or S*(1)) measured at τ in
    • a subset of the sample with Yτ = 0, and
    • all cases with Yτ = 0
  • Implications:
    • ncases, 1 = ncases, 1S
    • design parameter ncontrols, 1S replaced by probability p because ncontrols, 1S is a random variable in case-cohort designs

Illustration: CoRpower for trichotomous  S(1) and continuous  S*(1) | Bernoulli sampling

Trichotomous S(1) (Approach 1)

  • Scenario 10: vary p

Continuous S*(1)

  • Scenario 11: vary p

Scenario 10: vary  p (Approach 1) | Trichotomous  S(1), Bernoulli sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,             
                    nControlsTx = 3654,       
                    nCasesTxWithS = 32,       
                    cohort = TRUE,                # FALSE by default
                    p = c(0.01, 0.02, 0.03, 0.05),               
                    VEoverall = 0.75,             # overall VE
                    risk0 = 0.034,                # placebo-group endpoint risk from tau - taumax
                    VElat0 = seq(0, VEoverall, len=100), # grid of VE (V/P) among lower protected
                    VElat1 = rep(VEoverall, 100), # grid of VE (V/P) among medium protected
                    Plat0 = 0.2,                  # prevalence of lower protected
                    Plat2 = 0.6,                  # prevalence of higher protected
                    P0 = 0.2,                     # probability of low biomarker response
                    P2 = 0.6,                     # probability of high biomarker response
                    sens = 0.8, spec = 0.8, FP0 = 0, FN2 = 0,
                    M = 1000,                     # number of simulated clinical trials
                    alpha = 0.05,                 # two-sided Wald test Type 1 error rate
                    biomType = "trichotomous")    # "continuous" by default

Plot power curves with plotPowerTri()

plotPowerTri(outComputePower = pwr,  # 'computePower' output
             legendText = paste0("Cohort p = ", c(0.01, 0.02, 0.03, 0.05)))

Scenario 11: vary  p | Continuous  S*(1), Bernoulli sampling

Run simulations and compute power with computePower()

pwr <- computePower(nCasesTx = 32,             
                    nControlsTx = 3654,        
                    nCasesTxWithS = 32,       
                    cohort = TRUE,               # FALSE by default
                    p = c(0.01, 0.02, 0.03, 0.05),                  
                    VEoverall = 0.75,            # overall VE
                    risk0 = 0.034,               # placebo-group endpoint risk from tau - taumax
                    PlatVElowest = 0.2,          # prevalence of VE_lowest
                    VElowest = seq(0, VEoverall, len=100), # lowest VE for true biomarker X*<=nu
                    sigma2obs = 1,               # variance of observed biomarker S(1)
                    rho = 0.9                    # protection-relevant fraction of variance of S(1)
                    M = 1000,                    # number of simulated clinical trials
                    alpha = 0.05,                # two-sided Wald test Type 1 error rate
                    biomType = "continuous")     # "continuous" by default

Plot power curves with plotPowerCont()

plotPowerCont(outComputePower = pwr,  # 'computePower' output
              legendText = paste0("Cohort p = ", c(0.01, 0.02, 0.03, 0.05)))