--- 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)$**
  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$
  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$ 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*
  4. 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()` ```{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
  1. Specify overall $VE$ between $\tau$ and $\tau_{\mathrm{max}}$
  2. Specify $risk_0$
  3. Specify $P^{lat}_{lowestVE}$, $\rho$, and a grid of $VE_{lowest}$ values (e.g., ranging from $VE$ to 0)
  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*
    4. 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()` ```{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))) ```