# 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 $$\tau_{\mathrm{max}}$$

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

• $$Y = I(T \leq \tau_{\mathrm{max}})$$ is the binary outcome of interest

• $$\Delta$$ is the indicator that $$Y$$ is observed; i.e., $$\Delta = 0$$ if the participants drops out before $$\tau_{\mathrm{max}}$$ and before experiencing the primary endpoint and $$\Delta = 1$$ otherwise

• $$N_1$$ is the number of vaccine recipients observed or expected to be at risk at $$\tau$$ (typically, $$\tau$$ is the biomarker sampling timepoint)

• $$n_{cases,1}$$ ($$n_{controls,1}$$) is the number of observed or expected cases (controls) in vaccine recipients between $$\tau$$ and $$\tau_{\mathrm{max}}$$, where cases have $$\Delta Y = 1$$ and controls have $$\Delta (1-Y) = 1$$

• Note that both cases and controls have $$\Delta = 1$$
• $$n_{cases,1} + n_{controls,1} \leq N_1$$
• $$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

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

1. Specify true overall $$VE$$ between $$\tau$$ and $$\tau_{\mathrm{max}}$$
• Protocol-specified design alternative or $$\widehat{VE}$$
2. Specify $$risk_0 = P(Y=1 \mid Z=0, Y^{\tau}=0)$$ where $$Y^{\tau}=I(T \leq \tau)$$
• Protocol-specified placebo-group endpoint rate or $$\widehat{risk}_0$$
3. Select a grid of $$VE^{lat}_0$$ values
• E.g., ranging from $$VE$$ ($$H_0$$) to 0 (maximal $$H_1$$ not allowing harm by vaccination)
4. Select a grid of $$VE^{lat}_1 \geq VE^{lat}_0$$ values
• E.g., $$VE^{lat}_1 = VE$$
5. Specify $$P^{lat}_0$$ and $$P^{lat}_2$$, then $$P^{lat}_1 = 1-P^{lat}_0-P^{lat}_2$$
• Assuming $$risk^{lat}_0(x) = risk_0$$, $VE = VE^{lat}_0 P^{lat}_0 + VE^{lat}_1 P^{lat}_1 + VE^{lat}_2 P^{lat}_2$ yields $$VE^{lat}_2$$
• If $$VE^{lat}_0$$ varies from $$VE$$ to 0, then $$VE^{lat}_2$$ varies from VE to $$VE\, (P^{lat}_0 + P^{lat}_2)/P^{lat}_2$$
6. Specify $$P_0$$ and $$P_2$$, then $$P_1 = 1-P_0-P_2$$
• E.g., $$P_0 = P^{lat}_0$$ and $$P_2 = P^{lat}_2$$
7. Approach 1: Specify two of the three $$(Sens, FP^0, FP^1)$$ and two of the three $$(Spec, FN^2, FN^1)$$
• E.g., specify $$Sens$$ and $$Spec$$ and $$FP^0 = FN^2 = 0$$
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)$$

1. Assuming the classical measurement error model, where $$X^{\ast} \sim \mathsf{N}(0,\rho\sigma^2_{obs})$$, solve $P^{lat}_0 = P(X^{\ast} \leq \theta_0) \quad \textrm{and} \quad P^{lat}_2 = P(X^{\ast} > \theta_2)$ for $$\theta_0$$ and $$\theta_2$$

2. Generate $$B$$ realizations of $$X^{\ast}$$ and $$S^{\ast} = X^{\ast} + e$$, where $$e \sim \mathsf{N}(0,(1-\rho)\sigma^2_{obs})$$, and $$X^{\ast}$$ independent of $$e$$

• $$B = 20,000$$ by default

3. Using $$\theta_0$$ and $$\theta_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(\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.

4. 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.}$

• In Approach 2, plot $RR_t \quad \textrm{vs.} \quad \frac{RR^{lat}_2}{RR^{lat}_0},$ where $$RR_t$$ 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 $$\rho < 1$$, then $$RR_t$$ is closer to 1 than $$RR^{lat}_2\big/ RR^{lat}_0$$

• Note that, under the assumption of homogeneous risk in the placebo group, i.e., $$risk_0^{lat}(x) = risk_0$$, $$x=0,1,2$$, the relative risk ratio $$RR^{lat}_2\big/ RR^{lat}_0 = risk_1^{lat}(2) \big/ risk_1^{lat}(0)$$

• Consequently, $$risk_1(2) \big/ risk_1(0) > risk_1^{lat}(2) \big/ risk_1^{lat}(0)$$ because, if $$\rho < 1$$, then $$risk_1(2) > risk_1^{lat}(2)$$ and $$risk_1(0) < risk_1^{lat}(0)$$

8. Simulate $$M$$ data sets under the true parameter values:
Full data
1. $$N_x = (n_{controls,1} + n_{cases,1}) P^{lat}_x$$
2. $$\left(n_{cases,1}(0),n_{cases,1}(1),n_{cases,1}(2)\right) \sim \mathsf{Mult}(n_{cases,1},(p_0,p_1,p_2))$$, where $$p_k=P(X=k|Y=1,Y^{\tau}=0,Z=1)$$
3. For each of the $$N_x$$ subjects, generate $$S_i(1)$$ by a draw from $$\mathsf{Mult}(1,(p_0,p_1,p_2))$$, where $$p_k=P(S(1)=k|X=x)$$
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(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 $$H_0 \Leftrightarrow \{\widetilde{H}_0: \beta_S(1) \geq 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 \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$$

## Illustration: hypothetical randomized placebo-controlled VE trial

### Trial design

• $$N_{rand} = 4,100$$ participants randomized to each of the vaccine and placebo group and followed for $$\tau_{\mathrm{max}}=$$ 24 months
• Samples for measurement of $$S(1)$$ at $$\tau=$$ 3.5 months stored in all vaccine recipients
• Cumulative endpoint rate between $$\tau$$ and $$\tau_{\mathrm{max}}$$ in placebo group $$=3.4\%$$ $$(=risk_0)$$
• $$VE_{\tau-\tau_{\mathrm{max}}}=75\%$$, $$\quad VE_{0-\tau}=VE_{\tau-\tau_{\mathrm{max}}}/2$$
• Cumulative dropout rate between 0 and $$\tau_{\mathrm{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\mid Z=0 \sim \mathsf{Exp}(\lambda_T)$$ and $$C\mid Z=0 \sim \mathsf{Exp}(\lambda_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 \in (\tau,\tau_{\mathrm{max}}]$$

### Number of vaccine recipients observed to be at risk at $$\tau$$

\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 $$\tau$$ and $$\tau_{\mathrm{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 $$\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$$)$$

• Scenario 1: vary control:case ratio
• Scenario 2: vary $$Sens$$, $$Spec$$
• Scenario 3: vary $$P_0^{lat}$$, $$P_2^{lat}$$, $$P_0$$, $$P_2$$

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

• Scenario 4: vary $$\rho$$
• Scenario 5: vary $$P_0^{lat}$$, $$P_2^{lat}$$, $$P_0$$, $$P_2$$
• Scenario 6: vary $$n_{cases,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, $$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}$$

## 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$$
• Scenario 9: vary $$P_{lowestVE}^{lat}$$

## 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)$$)

• Bernoulli sample at baseline with sampling probability $$p$$
• $$S(1)$$ (or $$S^{\ast}(1)$$) measured at $$\tau$$ in
• a subset of the sample with $$Y^{\tau}=0$$, and
• all cases with $$Y^{\tau}=0$$
• Implications:
• $$n_{cases,1} = n^S_{cases,1}$$
• design parameter $$n^S_{controls,1}$$ replaced by probability $$p$$ because $$n^S_{controls,1}$$ is a random variable in case-cohort designs

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

Trichotomous $$S(1)$$ (Approach 1)

• Scenario 10: vary $$p$$

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