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:
computeN()
computePower()
plotPowerTri()
plotPowerCont()
plotRRgradVE()
plotVElatCont()
$n^S_{cases,1}$ ($n^S_{controls,1}$) is the number of observed cases (controls) in vaccine recipients between $\tau$ and $\tau_{\mathrm{max}}$ with measured biomarker response $S(1)$ (or $S^{\ast}(1)$)
If calculations done at design stage, then $N_1$, $n_{cases,1}$, $n_{controls,1}$, $n^S_{cases,1}$, and $n^S_{controls,1}$ are expected counts
Approach 1: Specify two of the three $(Sens, FP^0, FP^1)$ and two of the three $(Spec, FN^2, FN^1)$
Approach 2: Specify $\sigma^2_{obs}$ and $\rho = \sigma^2_{tr} / \sigma^2_{obs}$ and determine $(Sens, Spec, FP^0, FP^1, FN^1, FN^2)$ (see below) + Typically, $\sigma^2_{obs} = 1$ + Calculation of $(Sens, Spec, FP^0, FP^1, FN^1, FN^2)$
Estimate $Spec(\phi_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. <li> Find $\phi_0 = \phi^{\ast}_0$ and $\phi_2 = \phi^{\ast}_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.} \] </ol>
Observed data
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 \leq \alpha/2$ for specified $\alpha$
11. Repeat power calculation varying control:case ratio, $n{cases,1}$, $n^S_{cases,1}$, $(P^{lat}_0, P^{lat}_2, P_0, P_2)$, $(Sens, Spec)$, $\rho$
computeN()
[ \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} ]
[ \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} ]
[ \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} ]
[ \begin{align} n^S_{cases,1} &= n_{cases,1}\ n^S_{controls,1} &= K n^S_{cases,1} \end{align} ]
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)
CoRpower
for trichotomous (\, S(1)) | Without replacement samplingApproach 1 $(Sens, Spec, FP^0, FN^2$ specified$)$
Approach 2 $(\sigma^2_{obs}$ and $\rho$ specified$)$
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
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")))
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)))
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"))
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
plotPowerTri()
plotPowerTri(outComputePower = pwr, legendText = paste0("rho = ", c(1, 0.9, 0.7, 0.5)))
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)))
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))
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"))
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)))
CoRpower
for binary (\, S(1)) | Without replacement samplingAchieved 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$)$
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
plotPowerTri()
plotPowerTri(outComputePower = pwr, legendText = paste0("nCasesTx = ", c(25, 32, 35, 40)))
sample(
$\widetilde{x}^{\ast}$, prob=
$f_{X^{\ast}}(\widetilde{x}^{\ast}|Y=\cdot,Y^{\tau}=0,Z=1)$, replace=TRUE)
tps
in R package osDesign
)
CoRpower
for continuous (\, S^{\ast}(1)) | Without replacement samplingcomputePower()
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
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)))
plotVElatCont()
plotVElatCont()
plots the vaccine efficacy (VE) curve for the true biomarker X=x for eight different values of the true CoR relative risk, \eqn{RR_c (\rho=1)}, in vaccine recipients and the lowest vaccine efficacy level for the true biomarker, \eqn{VE_lowest}.
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 \eqn{VE_lowest}). 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")
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
plotPowerCont()
plotPowerCont(outComputePower = pwr, # output list of lists from 'computePower' legendText = paste0("PlatVElowest = ", c(0.05, 0.1, 0.15, 0.2)))
CoRpower
for trichotomous (\, S(1)) and continuous (\, S^{\ast}(1)) | Bernoulli samplingTrichotomous $S(1)$ (Approach 1)
Continuous $S^{\ast}(1)$
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
plotPowerTri()
plotPowerTri(outComputePower = pwr, # 'computePower' output legendText = paste0("Cohort p = ", c(0.01, 0.02, 0.03, 0.05)))
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
plotPowerCont()
plotPowerCont(outComputePower = pwr, # 'computePower' output legendText = paste0("Cohort p = ", c(0.01, 0.02, 0.03, 0.05)))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.