Processing math: 51%

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:

as well as two types of sampling design:

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:

Set-up and notation | Without replacement sampling

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

  1. Specify true overall VE between τ and τmax
    • Protocol-specified design alternative or ^VE
  2. Specify risk0=P(Y=1Z=0,Yτ=0) where Yτ=I(Tτ)
    • Protocol-specified placebo-group endpoint rate or ^risk0
  3. Select a grid of VElat0 values
    • E.g., ranging from VE (H0) to 0 (maximal H1 not allowing harm by vaccination)
  4. Select a grid of VElat1VElat0 values
    • E.g., VElat1=VE
  5. Specify Plat0 and Plat2, then Plat1=1Plat0Plat2
    • Assuming risklat0(x)=risk0, VE=VElat0Plat0+VElat1Plat1+VElat2Plat2 yields VElat2
    • If VElat0 varies from VE to 0, then VElat2 varies from VE to VE(Plat0+Plat2)/Plat2
  6. Specify P0 and P2, then P1=1P0P2
    • E.g., P0=Plat0 and P2=Plat2
  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 σ2obs and ρ=σ2tr/σ2obs and determine (Sens,Spec,FP0,FP1,FN1,FN2) (see below)
    • Typically, σ2obs=1

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

      1. Assuming the classical measurement error model, where XN(0,ρσ2obs), solve Plat0=P(Xθ0)andPlat2=P(X>θ2) for θ0 and θ2

      2. Generate B realizations of X and S=X+e, where eN(0,(1ρ)σ2obs), and X independent of e

        • B=20,000 by default

      3. Using θ0 and θ2 from the first step, define Spec(ϕ0)=P(Sϕ0Xθ0)FN1(ϕ0)=P(Sϕ0X(θ0,θ2])FN2(ϕ0)=P(Sϕ0X>θ2)Sens(ϕ2)=P(S>ϕ2X>θ2)FP1(ϕ2)=P(S>ϕ2X(θ0,θ2])FP0(ϕ2)=P(S>ϕ2Xθ0)

        Estimate Spec(ϕ0) by ^Spec(ϕ0)=#{Sbϕ0,Xbθ0}#{Xbθ0} etc.

      4. Find ϕ0=ϕ0 and ϕ2=ϕ2 that numerically solve P0=^Spec(ϕ0)Plat0+^FN1(ϕ0)Plat1+^FN2(ϕ0)Plat2P2=^Sens(ϕ2)Plat2+^FP1(ϕ2)Plat1+^FP0(ϕ2)Plat0 and compute Spec=^Spec(ϕ0),Sens=^Sens(ϕ2),etc.

    • In Approach 2, plot RRtvs.RRlat2RRlat0, where RRt is the CoR effect size defined as RRt:=risk1(2)risk1(0)=2x=0RRlatxP(X=x|S(1)=2)2x=0RRlatxP(X=x|S(1)=0)

    • If ρ<1, then RRt is closer to 1 than RRlat2/RRlat0

      • Note that, under the assumption of homogeneous risk in the placebo group, i.e., risklat0(x)=risk0, x=0,1,2, the relative risk ratio RRlat2/RRlat0=risklat1(2)/risklat1(0)

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

  8. Simulate M data sets under the true parameter values:
      Full data
    1. Nx=(ncontrols,1+ncases,1)Platx
    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
    5. Sample without replacement nScases,1 and nScontrols,1=KnScases,1 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{˜H0:β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, nScases,1, (Plat0,Plat2,P0,P2), (Sens,Spec), ρ

Illustration: hypothetical randomized placebo-controlled VE trial

Trial design

Illustration: calculation of input parameters with computeN()

Assumptions



  1. Failure time T and censoring time C are independent
  2. TZ=0Exp(λT) and CZ=0Exp(λC)
  3. RRττmax:=P(TτmaxT>τ,Z=1)P(TτmaxT>τ,Z=0)=P(TtT>τ,Z=1)P(TtT>τ,Z=0) for all t(τ,τmax]

Number of vaccine recipients observed to be at risk at τ

N1=NrandP(T>τ,C>τZ=1)=NrandP(T>τZ=1)P(C>τZ=1)=Nrand{1RR0τP(TτZ=0)}P(C>τZ=1)4,023

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

ncases,1=N1P(Tτmax,TCT>τ,C>τ,Z=1)=N1P(Tmin

Number of observed controls in vaccine recipients between \tau and \tau_{\mathrm{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 \tau and \tau_{\mathrm{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 N_1, n_{cases,1}, n_{controls,1}, and n^S_{cases,1} 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, FP^0, FN^2 specified)

Approach 2 (\sigma^2_{obs} and \rho specified)

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, RR_t, 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 \, P^{lat}_0, P^{lat}_2, P_0, P_2 (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 \, \rho (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 RR_t vs. RR_2^{lat}/RR_0^{lat} with plotRRgradVE()

plotRRgradVE() plots the ratio of relative risks for the higher and lower latent subgroups RR_2^{lat}/RR_0^{lat} against the CoR relative risk effect size RR_t = risk_1(2)/risk_1(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 \, P^{lat}_0, P_0, P^{lat}_2, P_2 (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 \, n_{cases,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 P_0^{lat}, P_2^{lat}, P_0, P_2 such that \begin{align} P_0^{lat} + P_2^{lat} &= 1\\ P_0 + P_2 &= 1 \end{align}

Approach 2 (\sigma^2_{obs} and \rho specified)

Scenario 7: vary \, n_{cases,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^{\ast}(1) | Without replacement sampling

  1. Specify overall VE between \tau and \tau_{\mathrm{max}}
    • Protocol-specified design alternative or \widehat{VE}
  2. Specify risk_0
    • Protocol-specified placebo-group endpoint rate or \widehat{risk}_0
  3. Specify P^{lat}_{lowestVE}, \rho, and a grid of VE_{lowest} values (e.g., ranging from VE to 0)
    • Fixed (VE, risk_0, P^{lat}_{lowestVE}, VE_{lowest}, \rho) 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 \alpha^{lat} and \beta^{lat}
    • Plot VE^{lat}_{x^{\ast}} vs. x^{\ast} and calculate the pertaining CoR effect size \exp(\beta^{lat}) for each level of VE_{lowest}
  4. Simulate M data sets under the true parameter values:
      Full data
    1. Sample X^{\ast} for n_{cases,1} cases from f_{X^{\ast}}(x^{\ast}|Y=1,Y^{\tau}=0,Z=1) using Bayes rule
    2. Sample X^{\ast} for n_{controls,1} controls from f_{X^{\ast}}(x^{\ast}|Y=0,Y^{\tau}=0,Z=1) using Bayes rule
      - How? Use a fine grid of \widetilde{x}^{\ast} values and then
      \;\;sample(\widetilde{x}^{\ast}, prob=f_{X^{\ast}}(\widetilde{x}^{\ast}|Y=\cdot,Y^{\tau}=0,Z=1), replace=TRUE)
    3. Sample S^{\ast}(1) following S^{\ast}(1) = X^{\ast} + e
    4. Observed data
    5. Sample without replacement n^S_{cases,1} and n^S_{controls,1} = K n^S_{cases,1} controls with measured S^{\ast}(1) (R=1)
  5. For each observed data set, compute the 1-sided one-degree-of-freedom Wald test statistic for H_0 \Leftrightarrow \{\widetilde{H}_0: \beta_{S^{\ast}(1)} \geq 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 \leq \alpha/2 for specified \alpha
  7. Repeat power calculation varying control:case ratio, n_{cases,1}, n^S_{cases,1}, P^{lat}_{lowestVE}, \rho

Illustration: CoRpower for continuous \, S^{\ast}(1) | Without replacement sampling

Scenario 8: vary \, \rho | Continuous \, S^{\ast}(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, RR_c, 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 VE^{lat}_{x^{\ast}} 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 VE^{lat}_0, VE^{lat}_1, 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 \, P_{lowestVE}^{lat} | Continuous \, S^{\ast}(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^{\ast}(1))

Illustration: CoRpower for trichotomous \, S(1) and continuous \, S^{\ast}(1) | Bernoulli sampling

Trichotomous S(1) (Approach 1)

Continuous S^{\ast}(1)

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^{\ast}(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)))