---
title: "Introduction to R Package `CoRpower`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Introduction to R Package CoRpower}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
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)$**
<ol type="i">
<li> 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$
<li> 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$
<ul>
<li> $B = 20,000$ by default
</ul>
<li> 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.
<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>
+ 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$
<ul>
<li> 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)$
<li> 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)$
</ul>
8. Simulate $M$ data sets under the true parameter values:
<ol type="a">
*Full data*
<li> $N_x = (n_{controls,1} + n_{cases,1}) P^{lat}_x$
<li> $\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)$
<li> 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)$ </li>
*Observed data*
<li> 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$
</ol>
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
<br><br>
<ol type="i">
<li> Failure time $T$ and censoring time $C$ are independent
<li> $T\mid Z=0 \sim \mathsf{Exp}(\lambda_T)$ and $C\mid Z=0 \sim \mathsf{Exp}(\lambda_C)$
<li> $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}}]$
</ol>
### 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()`
```{r}
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)
```
## 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()`
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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
```{r, eval=FALSE}
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
```
```{r, eval=FALSE}
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
```{r, eval=FALSE}
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
```
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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
```{r, eval=FALSE}
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
```
```{r, eval=FALSE}
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
```{r, eval=FALSE}
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
```
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
plotPowerTri(outComputePower = pwr,
legendText = paste0("nCasesTx = ", c(25, 32, 35, 40)))
```
## Algorithm for continuous biomarker \(\, S^{\ast}(1)\) | Without replacement sampling
<ol>
<li> Specify overall $VE$ between $\tau$ and $\tau_{\mathrm{max}}$
<ul>
<li> Protocol-specified design alternative or $\widehat{VE}$
</ul>
</li>
<li> Specify $risk_0$
<ul>
<li> Protocol-specified placebo-group endpoint rate or $\widehat{risk}_0$
</ul>
</li>
<li> Specify $P^{lat}_{lowestVE}$, $\rho$, and a grid of $VE_{lowest}$ values (e.g., ranging from $VE$ to 0)
<ul>
<li> 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}$
<li> Plot $VE^{lat}_{x^{\ast}}$ vs. $x^{\ast}$ and calculate the pertaining CoR effect size $\exp(\beta^{lat})$ for each level of $VE_{lowest}$
</ul>
</li>
<li> Simulate $M$ data sets under the true parameter values:
<ol type="a">
*Full data*
<li> Sample $X^{\ast}$ for $n_{cases,1}$ cases from $f_{X^{\ast}}(x^{\ast}|Y=1,Y^{\tau}=0,Z=1)$ using Bayes rule
<li> Sample $X^{\ast}$ for $n_{controls,1}$ controls from $f_{X^{\ast}}(x^{\ast}|Y=0,Y^{\tau}=0,Z=1)$ using Bayes rule
<br>
- How? Use a fine grid of $\widetilde{x}^{\ast}$ values and then
<br>
$\;\;$`sample(`$\widetilde{x}^{\ast}$, `prob=`$f_{X^{\ast}}(\widetilde{x}^{\ast}|Y=\cdot,Y^{\tau}=0,Z=1)$, `replace=TRUE)`
<li> Sample $S^{\ast}(1)$ following $S^{\ast}(1) = X^{\ast} + e$ </li>
*Observed data*
<li> 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)$
</ol>
</li>
<li> 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`)
<li> Compute power as proportion of data sets with 1-sided Wald test $p \leq \alpha/2$ for specified $\alpha$
<li> Repeat power calculation varying control:case ratio, $n_{cases,1}$, $n^S_{cases,1}$, $P^{lat}_{lowestVE}$, $\rho$
</ol>
## 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()`
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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.
```{r, eval=FALSE}
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, \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`:
```{r, eval=FALSE}
plotVElatCont(outComputePower = pwr)
```
Using the function when `computePower()` output is saved in a file with name "myFile" and location "~/myDir":
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
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()`
```{r, eval=FALSE}
plotPowerCont(outComputePower = pwr, # 'computePower' output
legendText = paste0("Cohort p = ", c(0.01, 0.02, 0.03, 0.05)))
```