--- title: "RMSTpowerBoost: Sample Size and Power Calculations for RMST-based Clinical Trials" name: RMSTpowerBoost-Main css: css/Main.css output: rmarkdown::html_vignette: toc: true fig_caption: true code_folding: hide df_print: paged highlight: null self_contained: true bibliography: references.bib biblio-style: apalike link-citations: yes vignette: > %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{RMSTpowerBoost: Sample Size and Power Calculations for RMST-based Clinical Trials} %\VignetteDepends{RMSTpowerBoost, survival, dplyr, tidyr, knitr, ggplot2, mgcv, kableExtra} --- ```{r setup, include=FALSE, purl=FALSE} cache_root <- Sys.getenv("RMSTPOWERBOOST_VIGNETTE_CACHE", unset = "") if (!nzchar(cache_root)) { # R CMD build renders vignettes in a temporary copy; local builds need a stable cache path. cache_root <- if (nzchar(Sys.getenv("_R_CHECK_PACKAGE_NAME_"))) { "cache" } else { file.path(tools::R_user_dir("RMSTpowerBoost", "cache"), "vignettes") } } cache_path <- file.path(cache_root, "RMSTpowerBoost-Main", "html") dir.create(cache_path, recursive = TRUE, showWarnings = FALSE) cache_path <- paste0(normalizePath(cache_path, winslash = "/", mustWork = FALSE), "/") pkg_root <- if (file.exists("../DESCRIPTION") && dir.exists("../R")) ".." else "." pkg_version <- read.dcf(file.path(pkg_root, "DESCRIPTION"))[1, "Version"] r_files <- list.files(file.path(pkg_root, "R"), pattern = "[.]R$", recursive = TRUE, full.names = TRUE) r_hash <- paste(unname(tools::md5sum(r_files)), collapse = "") cache_extra <- paste( getRversion(), pkg_version, r_hash, sep = "-" ) object_cache_path <- file.path( cache_root, "RMSTpowerBoost-Main", "objects", paste(pkg_version, substr(r_hash, 1, 16), sep = "-") ) dir.create(object_cache_path, recursive = TRUE, showWarnings = FALSE) knitr::opts_chunk$set( message = FALSE, cache = TRUE, cache.path = cache_path, cache.extra = cache_extra, cache.lazy = FALSE, purl = FALSE, warning = FALSE, comment = NA, fig.align = 'center', fig.width = 7, fig.height = 5, table.width = "100%", table.align = "center", collapse = TRUE ) packages = c("survival", "dplyr", "tidyr", "knitr", "ggplot2", "mgcv", "kableExtra") lapply(packages, require, character.only = TRUE) library(RMSTpowerBoost) run_full_vignettes <- identical(Sys.getenv("RMSTPOWERBOOST_FULL_VIGNETTES"), "true") vignette_power_result <- function(results_data, x_col, title, x_label, summary = NULL, color = "#0072B2") { plot <- ggplot2::ggplot(results_data, ggplot2::aes_string(x = x_col, y = "Power")) + ggplot2::geom_line(color = color, linewidth = 1) + ggplot2::geom_point(color = color, size = 3) + ggplot2::labs(title = title, x = x_label, y = "Estimated Power") + ggplot2::ylim(0, 1) + ggplot2::theme_minimal() list( results_data = results_data, results_plot = plot, results_summary = summary, model_output = list() ) } vignette_ss_result <- function(results_data, plot_data, x_col, target_power, required_n, title, x_label, summary = NULL, color = "#009E73") { plot <- ggplot2::ggplot(plot_data, ggplot2::aes_string(x = x_col, y = "Power")) + ggplot2::geom_line(color = color, linewidth = 1) + ggplot2::geom_point(color = color, size = 3) + ggplot2::geom_hline(yintercept = target_power, linetype = "dashed", color = "red") + ggplot2::geom_vline(xintercept = required_n, linetype = "dotted", color = "blue") + ggplot2::labs(title = title, x = x_label, y = "Calculated Power") + ggplot2::theme_minimal() list( results_data = results_data, results_plot = plot, results_summary = summary, model_output = list() ) } vignette_summary <- function(statistic, value) { data.frame(Statistic = statistic, Value = value, row.names = NULL) } precomputed_vignette_result <- function(label) { switch( label, veteran_power_calc = vignette_power_result( data.frame(N_per_Arm = c(100, 150, 200, 250), Power = c(0.1266, 0.1687, 0.2106, 0.2521)), "N_per_Arm", "Power Curve: Linear IPCW RMST Model", "Sample Size Per Arm", vignette_summary("Assumed RMST Difference (from pilot)", -9.7429), "#D55E00" ), veteran_ss_calc = vignette_ss_result( data.frame(Target_Power = 0.40, Required_N_per_Arm = 4500), data.frame(N_per_Arm = seq(1000, 4500, by = 250), Power = c(0.1281, 0.1496, 0.1711, 0.1924, 0.2137, 0.2348, 0.2559, 0.2768, 0.2975, 0.3180, 0.3383, 0.3583, 0.3780, 0.3975, 0.4166)), "N_per_Arm", 0.40, 4500, "Sample Size Search Path: Linear IPCW RMST Model", "Sample Size Per Arm", vignette_summary("Assumed RMST Difference (from pilot)", -3.9666) ), linear_boot_example = vignette_power_result( data.frame(N_per_Arm = c(150, 200, 250), Power = c(0.08, 0.10, 0.12)), "N_per_Arm", "Power Curve: Linear IPCW RMST Model", "Sample Size Per Arm", vignette_summary( c("Mean RMST Difference", "Mean Standard Error", "95% CI Lower", "95% CI Upper"), c(-2.7125, 9.2260, -23.4963, 18.0712) ), "#D55E00" ), linear_boot_ss_example = vignette_ss_result( data.frame(Target_Power = 0.50, Required_N_per_Arm = 150), data.frame(N_per_Arm = c(50, 75, 100, 125, 150), Power = c(0.14, 0.27, 0.28, 0.42, 0.54)), "N_per_Arm", 0.50, 150, "Sample Size Search Path: Linear IPCW RMST Model", "Sample Size Per Arm", vignette_summary( c("Mean RMST Difference", "Mean Standard Error", "95% CI Lower", "95% CI Upper"), c(-12.6216, 6.3372, -24.6968, -0.5464) ) ), colon_ss_calc = vignette_ss_result( data.frame(Target_Power = 0.60, Required_N_per_Stratum = 2200), data.frame(N_per_Stratum = seq(100, 2200, by = 100), Power = c(0.0691, 0.0994, 0.1285, 0.1573, 0.1859, 0.2144, 0.2427, 0.2707, 0.2985, 0.3259, 0.3528, 0.3793, 0.4052, 0.4306, 0.4553, 0.4794, 0.5028, 0.5256, 0.5476, 0.5689, 0.5895, 0.6094)), "N_per_Stratum", 0.60, 2200, "Sample Size Search Path: Additive Stratified RMST Model", "Sample Size Per Stratum", vignette_summary("Assumed RMST Difference (from pilot)", -36.7735) ), additive_power_calc = vignette_power_result( data.frame(N_per_Stratum = c(1000, 3000, 5000), Power = c(0.3259, 0.7432, 0.9213)), "N_per_Stratum", "Power Curve for Additive Stratified Model", "Sample Size Per Stratum", NULL ), ms_power_analytical_example = vignette_power_result( data.frame(N_per_Stratum = c(300, 400, 500), Power = c(0.5062, 0.6259, 0.7225)), "N_per_Stratum", "Power Curve: Multiplicative Stratified RMST Model", "Sample Size Per Stratum", NULL, "#E69F00" ), ms_ss_analytical_example = vignette_ss_result( data.frame(Target_Power = 0.60, Required_N_per_Stratum = 400), data.frame(N_per_Stratum = seq(50, 400, by = 25), Power = c(0.1244, 0.1655, 0.2063, 0.2467, 0.2867, 0.3259, 0.3643, 0.4016, 0.4378, 0.4726, 0.5062, 0.5383, 0.5690, 0.5982, 0.6259)), "N_per_Stratum", 0.60, 400, "Sample Size Search Path: Multiplicative Stratified RMST Model", "Sample Size Per Stratum", vignette_summary("Assumed log(RMST Ratio) (from pilot)", -0.0898) ), ms_power_boot_example = vignette_power_result( data.frame(N_per_Stratum = c(100, 300, 500), Power = c(0.2333, 0.6333, 0.9667)), "N_per_Stratum", "Power Curve: Multiplicative Stratified RMST Model", "Sample Size Per Stratum", vignette_summary( c("Mean RMST Ratio", "95% CI Lower", "95% CI Upper"), c(1.0001, 0.9579, 1.0501) ), "#E69F00" ), ms_ss_boot_example = vignette_ss_result( data.frame(Target_Power = 0.75, Required_N_per_Stratum = 400), data.frame(N_per_Stratum = seq(100, 400, by = 50), Power = c(0.4667, 0.4667, 0.6333, 0.6667, 0.6667, 0.6667, 0.8000)), "N_per_Stratum", 0.75, 400, "Sample Size Search Path: Multiplicative Stratified RMST Model", "Sample Size Per Stratum", vignette_summary( c("Mean RMST Ratio", "95% CI Lower", "95% CI Upper"), c(1.0081, 0.9538, 1.0684) ) ), gbsg_power_calc = vignette_power_result( data.frame(N_per_Group = c(50, 200, 400), Power = c(1, 1, 1)), "N_per_Group", "Power Curve: Additive GAM RMST Model", "Sample Size Per Arm", vignette_summary( c("Mean RMST Difference", "Mean Standard Error", "95% CI Lower", "95% CI Upper"), c(861.2639, 30.3843, 786.0199, 936.5079) ) ), gbsg_ss_calc = vignette_ss_result( data.frame(Target_Power = 0.95, Required_N_per_Group = 50), data.frame(N_per_Group = 50, Power = 1), "N_per_Group", 0.95, 50, "Sample Size Search Path: Additive GAM RMST Model", "Sample Size Per Arm", vignette_summary( c("Mean RMST Difference", "Mean Standard Error", "95% CI Lower", "95% CI Upper"), c(90.6345, 0.3314, 89.7268, 91.5421) ) ), dc_power_calc = vignette_power_result( data.frame(N_per_Arm = c(100, 250, 500), Power = c(0.4760, 0.8517, 0.9889)), "N_per_Arm", "Power Curve: Dependent Censoring RMST Model", "Sample Size Per Arm", vignette_summary("Assumed RMST Difference (from pilot)", -10.6135) ), mgus2_ss_calc = vignette_ss_result( data.frame(Target_Power = 0.80, Required_N_per_Arm = 250), data.frame(N_per_Arm = c(100, 150, 200, 250), Power = c(0.4760, 0.6431, 0.7663, 0.8517)), "N_per_Arm", 0.80, 250, "Sample Size Search Path: Dependent Censoring RMST Model", "Sample Size Per Arm", vignette_summary("Assumed RMST Difference (from pilot)", -10.6135) ), NULL ) } vignette_cache <- function(label, value) { path <- file.path(object_cache_path, paste0(label, ".rds")) if (file.exists(path)) { readRDS(path) } else if (!run_full_vignettes) { result <- precomputed_vignette_result(label) if (is.null(result)) force(value) else result } else { result <- force(value) saveRDS(result, path) result } } ``` # Introduction Clinical trials with time-to-event endpoints often rely on hazard-based methods such as the proportional hazards (PH) model and the hazard ratio (HR). The HR can be hard to interpret when treatment effects vary over time, and the proportional hazards assumption is frequently violated in practice. The **Restricted Mean Survival Time (RMST)** is the expected event-free time up to a pre-specified follow-up point, $L$ [@royston2013; @uno2014]. It yields a treatment contrast on the time scale, such as an average gain of several months over a clinically relevant horizon. Recent work models RMST directly as a function of baseline covariates instead of estimating it only from a survival curve or hazard model. Methods based on Inverse Probability of Censoring Weighting (IPCW) [@tian2014] now cover stratified studies [@wang2019; @zhang2024] and covariate-dependent censoring [@wang2018]. Most available software emphasizes estimation from existing data rather than study design. As a result, trial statisticians often need custom code for sample size and power calculations. `RMSTpowerBoost` implements direct RMST methods for **power and sample size calculations**: * **Linear IPCW Models**: A direct regression model for RMST with IPCW-based estimation [@tian2014]. * **Stratified Models**: Efficient methods for studies with a large number of strata (e.g., clinical centers), including both **additive** [@zhang2024] and **multiplicative** [@wang2019] models. * **Dependent Censoring Models**: Methods for handling covariate-dependent censoring [@wang2018]. * **Flexible Non-Linear Models**: Bootstrap-based functions using Generalized Additive Models (GAMs) to capture non-linear covariate effects. * **Analytic vs. Bootstrap Methods**: For most models, the package offers a choice between a fast `analytical` calculation and a simulation-based `boot` method. This vignette outlines the main model families and shows how to use them. To keep CRAN checks fast, this vignette uses precomputed outputs for the long-running examples by default. Set `RMSTPOWERBOOST_FULL_VIGNETTES=true` before knitting to recompute every analytical and bootstrap result. ------------------------------------------------------------------------ # Core Concepts of `RMSTpowerBoost` Package {.tabset} The functions in this package are grounded in a regression-based formulation of the **restricted mean survival time (RMST)**. For a given subject \( i \) with event time \(T_i\), covariate vector \( Z_i \) and treatment indicator \(\mathrm{Trt}_i\), the conditional RMST is modeled as $$ \mathbb{E}[\min(T_i, L)\mid Z_i] \;=\; \beta_0 + \beta_{\text{effect}}\,\mathrm{Trt}_i + \beta_2^{\top} Z_i, $$ where \(L\) is the restriction time, and \( \beta_{\text{effect}} \) represents the modeled treatment contrast on the RMST scale, that is, the expected difference in event-free time between treatment arms after adjustment for the included covariates. The quantity \( \beta_{\text{effect}} \) therefore defines the **effect size** used throughout the analytical power and sample size functions in `RMSTpowerBoost`. ## The Analytic Method The analytical functions estimate power from a closed-form expression, so they are useful for rapid scenario exploration. The process is: 1. **One-Time Estimation**: The function first analyzes the provided reference data to estimate two key parameters: * The **treatment effect size** (e.g., the difference in RMST or the log-RMST ratio). * The **asymptotic variance** of that effect estimator, which measures its uncertainty. 2. **Power Formula**: It then plugs these fixed estimates into a standard power formula. For a given total sample size `N`, the power is calculated as: $$ \text{Power} = \Phi\left( \frac{|\beta_{\text{effect}}|}{\sigma_N} - z_{1-\alpha/2} \right) $$ where: * $\Phi$ is the cumulative distribution function (CDF) of the standard normal distribution. * $\beta_{\text{effect}}$ is the treatment effect. * $\sigma_N = \frac{\sigma_1}{\sqrt{N}}$ is the standard error of the effect for the target sample size `N`, which is scaled from the reference data's variance. * $z_{1-\alpha/2}$ is the critical value from the standard normal distribution (e.g., 1.96 for an alpha of 0.05). ## The Bootstrap Method The bootstrap functions estimate power empirically and therefore require more computation. The process is: 1. **Resample**: The function simulates a "future trial" of a given `sample_size` by resampling with replacement from the reference data. 2. **Fit Model**: On this new bootstrap sample, it performs the full analysis (e.g., calculating weights or pseudo-observations and fitting the specified model). 3. **Get P-Value**: It extracts the p-value for the treatment effect from the fitted model. 4. **Repeat**: This process is repeated thousands of times (`n_sim`). 5. **Calculate Power**: The final estimated power is the proportion of simulations where the p-value was less than the significance level `alpha`. $$ \text{Power} = \frac{\text{Number of simulations with } p < \alpha}{n_{\text{sim}}} $$ ## The Sample Size Search Algorithm `rmst.ss()` uses an iterative search algorithm to find the `N` required to achieve a `target_power`: 1. **Start**: The search begins with a sample size of `n_start`. 2. **Calculate Power**: It calculates the power for the `current_n` using either the **analytic formula** or a **full bootstrap simulation**. 3. **Check Condition**: * If `calculated_power >= target_power`, the search succeeds and returns `current_n`. * If not, it increments the sample size (`current_n = current_n + n_step`) and repeats the process. 4. **Stopping Rules**: The search terminates if the sample size exceeds `max_n` or, for bootstrap methods, if the power fails to improve for a set number of `patience` steps. ## The Unified Interface `RMSTpowerBoost` exposes all models through two top-level functions that use a familiar `Surv()`-based formula interface: | Function | Purpose | | :--- | :--- | | `rmst.power(Surv(time, status) ~ covariates, data, arm, sample_sizes, L, ...)` | Power curve over a vector of sample sizes | | `rmst.ss(Surv(time, status) ~ covariates, data, arm, target_power, L, ...)` | Minimum sample size to reach a power target | Key arguments that control model selection: | Argument | Values | Effect | | :--- | :--- | :--- | | `type` | `"analytical"` (default) / `"boot"` | Analytic vs. bootstrap method | | `strata` | column name / `~col` / `NULL` | Activates stratified model | | `strata_type` | `"additive"` / `"multiplicative"` | Stratified model type | | `dep_cens` | `TRUE` / `FALSE` | Dependent-censoring model | | Smooth terms | `s(var)` in formula | Activates GAM (forces `type = "boot"`) | The underlying model-specific functions remain exported for direct use when needed. ## Selecting an Appropriate Model Model choice depends on the study design and the assumptions you are willing to make. The table below summarizes typical settings for each approach: | Model | Key Assumption / Scenario | Use when | | :---------------------------- | :--------------------------------------------------------------------------- | :------------------------------------------------------------------------------------------------------------------------------------------------------------- | | **Linear IPCW** | Assumes a linear relationship between covariates and RMST. | there is no strong evidence of non-linear effects or complex stratification. | | **Additive Stratified** | Assumes the treatment adds a constant amount of survival time across strata. | treatment effects are expected to be comparable across centers or other strata. | | **Multiplicative Stratified** | Assumes the treatment multiplies survival time proportionally across strata. | treatment effects are better expressed as relative changes in RMST across strata. | | **Semiparametric GAM** | Allows for non-linear covariate effects on RMST. | covariates such as age or biomarkers are likely to have non-linear associations with the outcome. | | **Dependent Censoring** | Accounts for covariate-dependent censoring under a single censoring mechanism. | censoring depends on measured covariates and competing risks are not being modeled explicitly. | # Linear IPCW Models{.tabset} These functions implement the foundational direct linear regression model for the RMST. This model is appropriate when a linear relationship between covariates and the RMST is assumed, and when censoring is independent of the event of interest. ## Theory and Model Based on the methods of [@tian2014], these functions model the conditional RMST as a linear function of covariates: $$\mathbb{E}[\min(T_i, L) | Z_i] = \beta_0 + \beta_{\text{effect}} \text{Treatment}_i + \beta_2 \text{Covariate}_{i}$$ In this model, the expected RMST up to a pre-specified time **L** for subject *i* is modeled as a linear combination of their treatment arm and other variables $Z_i$. To handle right-censoring, the method uses **Inverse Probability of Censoring Weighting (IPCW)**. This is achieved through the following steps: 1. A survival curve for the **censoring distribution** is estimated using the Kaplan-Meier method (where "failure" is being censored). 2. For each subject who experienced the primary event, a weight is calculated. This weight is the inverse of the probability of *not* being censored up to their event time. 3. A standard weighted linear model (`lm()`) is then fitted using these weights. The model only includes subjects who experienced the event. ## Analytical Methods The analytical functions use a formula based on the asymptotic variance of the regression coefficients to calculate power or sample size, making them fast to evaluate. **Scenario**: We use the `veteran` dataset to estimate power for a trial comparing standard vs. test chemotherapy (`trt`), adjusting for the Karnofsky performance score (`karno`). ### Power Calculation First, let's inspect the prepared `veteran` dataset. ```{r veteran_data_prep, echo=FALSE} vet <- veteran %>% mutate( arm = ifelse(trt == 1, 0, 1), status = status ) head(vet) ``` Now, we calculate the power for a range of sample sizes using a truncation time of 9 months (270 days). ```{r veteran_power_calc} power_results_vet <- vignette_cache("veteran_power_calc", { rmst.power( Surv(time, status) ~ karno, data = vet, arm = "arm", sample_sizes = c(100, 150, 200, 250), L = 270 ) }) ``` The results are returned as a data frame and a `ggplot` object. ```{r veteran_table_plot, echo=FALSE} kbl(power_results_vet$results_data , caption = "Power Analysis for Veteran Dataset") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ``` ### Sample Size Calculation We can also use the analytical method to find the required sample size to achieve a target power for a truncation time of one year (365 days). ```{r veteran_ss_calc} ss_results_vet <- vignette_cache("veteran_ss_calc", { rmst.ss( Surv(time, status) ~ karno, data = vet, arm = "arm", target_power = 0.40, L = 365, n_start = 1000, n_step = 250, max_n = 5000 ) }) ``` ```{r veteran_ss_table, echo=FALSE} kbl(ss_results_vet$results_summary, caption = "Estimated Effect from Reference Data") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ss_results_vet$results_plot + theme_bw(base_size = 14) ``` ## Bootstrap Methods Passing `type = "boot"` to `rmst.power()` or `rmst.ss()` switches to a simulation-based approach. This method repeatedly resamples from the reference data, fits the model on each sample, and calculates power as the proportion of simulations where the treatment effect is significant. It relies less on closed-form large-sample approximations at the cost of greater computation time. ### Power and Sample Size Calculation (bootstrap) Here is how you would call the bootstrap method for power for the linear model. The following examples use the same `veteran` dataset, but with a smaller number of simulations for demonstration purposes. In practice, a larger number of simulations (e.g., 1,000 or more) is recommended to ensure stable results. First we calculate the power for a range of sample sizes. ```{r linear_boot_example} power_boot_vet <- vignette_cache("linear_boot_example", { rmst.power( Surv(time, status) ~ karno, data = vet, arm = "arm", sample_sizes = c(150, 200, 250), L = 365, type = "boot", n_sim = 50 ) }) ``` ```{r echo=FALSE} power_boot_vet$results_plot ``` Here is how you would call the bootstrap method for sample size calculation, targeting 50% power and truncating at 180 days (6 months). ```{r} ss_boot_vet <- vignette_cache("linear_boot_ss_example", { rmst.ss( Surv(time, status) ~ karno, data = vet, arm = "arm", target_power = 0.5, L = 180, type = "boot", n_sim = 100, patience = 5 ) }) ``` ```{r echo=FALSE} ss_boot_vet$results_plot + theme_bw(base_size = 14) ``` *** # Additive Stratified Models {.tabset} In multi-center clinical trials, the analysis is often stratified by a categorical variable with many levels, such as clinical center or a discretized biomarker. Estimating a separate parameter for each stratum can be inefficient when the number of strata is large. The additive stratified model removes the stratum-specific effects through conditioning. ## Theory and Model The semiparametric additive model for RMST, as developed by [@zhang2024], is defined as: $$\mu_{ij} = \mu_{0j} + \beta'Z_i$$ This model assumes that the effect of the covariates $Z_i$ (which includes the treatment arm) is **additive** and constant across all strata $j$. Each stratum has its own baseline RMST, denoted by $\mu_{0j}$. The common treatment effect, $\beta$, is estimated with a **stratum-centering** approach applied to IPCW-weighted data. This avoids direct estimation of the many $\mu_{0j}$ parameters. ## Analytical Methods ### Sample Size Calculation **Scenario**: We use the `colon` dataset to design a trial stratified by the extent of local disease (`extent`), a factor with 4 levels. We want to find the sample size per stratum to achieve 60% power. Let's inspect the prepared `colon` dataset. ```{r colon_data_prep, echo=FALSE} colon_death <- colon %>% filter(etype == 2) %>% select(time, status, rx, extent) %>% na.omit() %>% mutate( arm = ifelse(rx == "Obs", 0, 1), status = status, strata = factor(extent) ) head(colon_death) ``` Now, we run the sample size search for 60% power, truncating at 5 years (1825 days). ```{r colon_ss_calc} ss_results_colon <- vignette_cache("colon_ss_calc", { rmst.ss( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "additive", target_power = 0.60, L = 1825, n_start = 100, n_step = 100, max_n = 10000 ) }) ``` ```{r colon_ss_table, echo=FALSE} kbl(ss_results_colon$results_summary , caption = "Estimated Effect from Reference Data") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") final_n_colon <- ss_results_colon$results_data$Required_N_per_Stratum power_at_final_n_colon <- ss_results_colon$results_plot$data %>% filter(N_per_Stratum == final_n_colon) %>% pull(Power) ss_results_colon$results_plot ``` ### Power Calculation This example calculates the power for a given set of sample sizes in a stratified additive model using the same `colon` dataset. ```{r additive_power_calc} power_results_colon <- vignette_cache("additive_power_calc", { rmst.power( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "additive", sample_sizes = c(1000, 3000, 5000), L = 1825 ) }) ``` ```{r additive_power_table_plot, echo=FALSE} kbl(power_results_colon$results_data, caption = "Power for Additive Stratified Colon Trial") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") power_results_colon$results_plot + geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") + labs(title = "Power Curve for Additive Stratified Model") + theme_bw(base_size = 14) ``` *** # Multiplicative Stratified Models {.tabset} Use the multiplicative model when treatment is expected to act on the RMST scale through a relative effect, such as a percentage increase or decrease in survival time. ## Theory and Model The multiplicative model, based on the work of [@wang2019], is defined as: $$\mu_{ij} = \mu_{0j} \exp(\beta'Z_i)$$ In this model, the covariates $Z_i$ have a **multiplicative** effect on the baseline stratum-specific RMST, $\mu_{0j}$. This structure is equivalent to a linear model on the log-RMST. Formal estimation of $\beta$ requires an iterative solver. This package instead fits a weighted log-linear model (`lm(log(Y_rmst) ~ ...)`) to approximate the log-RMST ratio and its variance. ## Analytical Methods ### Power Calculation This function calculates the power for various sample sizes using the analytical method for the multiplicative stratified model. ```{r ms_power_analytical_example} power_ms_analytical <- vignette_cache("ms_power_analytical_example", { rmst.power( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "multiplicative", sample_sizes = c(300, 400, 500), L = 1825 ) }) ``` ```{r echo=FALSE} kbl(power_ms_analytical$results_data, caption = "Power for Multiplicative Stratified Model") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ``` ### Sample Size Calculation The following example demonstrates the sample size calculation for the multiplicative model. ```{r ms_ss_analytical_example} ms_ss_results_colon <- vignette_cache("ms_ss_analytical_example", { rmst.ss( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "multiplicative", target_power = 0.6, L = 1825 ) }) ``` ```{r echo=FALSE} kbl(ms_ss_results_colon$results_summary, caption = "Sample Size for Multiplicative Stratified Model") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ms_ss_results_colon$results_plot + theme_bw(base_size = 14) ``` ## Bootstrap Methods The bootstrap approach provides a simulation-based analysis for the multiplicative model. Pass `type = "boot"` together with `strata_type = "multiplicative"`. ### Power Calculation (bootstrap) ```{r ms_power_boot_example} power_ms_boot <- vignette_cache("ms_power_boot_example", { rmst.power( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "multiplicative", sample_sizes = c(100, 300, 500), L = 1825, type = "boot", n_sim = 30, parallel.cores = 1 ) }) ``` ```{r echo=FALSE} kbl(power_ms_boot$results_summary, caption = "Power for Multiplicative Stratified Model (Bootstrap)") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") power_ms_boot$results_plot ``` ### Sample Size Calculation (bootstrap) ```{r ms_ss_boot_example} ss_ms_boot <- vignette_cache("ms_ss_boot_example", { rmst.ss( Surv(time, status) ~ 1, data = colon_death, arm = "arm", strata = "strata", strata_type = "multiplicative", target_power = 0.75, L = 1825, type = "boot", n_sim = 30, n_start = 100, n_step = 50, patience = 4, parallel.cores = 1 ) }) ``` ```{r echo=FALSE} kbl(ss_ms_boot$results_summary, caption = "Sample Size for Multiplicative Stratified Model (Bootstrap)") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ``` *** # Semiparametric GAM Models {.tabset} When a covariate is expected to have a non-linear effect on the outcome, standard linear models may be misspecified. Generalized Additive Models (GAMs) handle this by fitting smooth functions. ## Theory and Model These functions use a bootstrap simulation approach combined with a GAM. The method involves two main steps: 1. **Jackknife Pseudo-Observations**: The time-to-event outcome is first converted into **jackknife pseudo-observations** for the RMST. This technique [@andersen2004; @parner2010] creates a continuous, uncensored variable that represents each subject's contribution to the RMST. This makes the outcome suitable for use in a standard regression framework. 2. **GAM Fitting**: A GAM is then fitted to these pseudo-observations. The model has the form: $$\mathbb{E}[\text{pseudo}_i] = \beta_0 + \beta_{\text{effect}} \cdot \text{Treatment}_i + \sum_{k=1}^{q} f_k(\text{Covariate}_{ik})$$ Here, $f_k()$ are the non-linear **smooth functions** (splines) that the GAM estimates from the data. ### Power Calculation Formula Because this is a bootstrap method, power is not calculated from a direct formula but is instead estimated empirically from the simulations: $$\text{Power} = \frac{1}{B} \sum_{b=1}^{B} \mathbb{I}(p_b < \alpha)$$ Where: * $B$ is the total number of bootstrap simulations (`n_sim`). * $p_b$ is the p-value for the treatment effect in the $b$-th simulation. * $\mathbb{I}(\cdot)$ is the indicator function, which is 1 if the condition is true and 0 otherwise. ## Bootstrap Methods Wrap smooth covariates in `s()` in the formula — this automatically routes to the GAM bootstrap engine regardless of the `type` argument. ### Power Calculation **Scenario**: We use the `gbsg` (German Breast Cancer Study Group) dataset, suspecting that the progesterone receptor count (`pgr`) has a non-linear effect on recurrence-free survival. Here is a look at the prepared `gbsg` data. ```{r gbsg_data_prep, echo=FALSE} gbsg_prepared <- gbsg %>% mutate( arm = ifelse(hormon == "no", 0, 1) ) head(gbsg_prepared) ``` The following code shows how to calculate power. Wrapping `pgr` in `s()` signals a smooth non-linear effect. ```{r gbsg_power_calc} power_gam <- vignette_cache("gbsg_power_calc", { rmst.power( Surv(rfstime, status) ~ s(pgr), data = gbsg_prepared, arm = "arm", sample_sizes = c(50, 200, 400), L = 2825, n_sim = 50, parallel.cores = 1 ) }) print(power_gam$results_plot) ``` ### Sample Size Calculation **Scenario**: We want to find the sample size needed to achieve 95% power for detecting an effect of `pgr` on recurrence-free survival. ```{r} ss_gam <- vignette_cache("gbsg_ss_calc", { rmst.ss( Surv(rfstime, status) ~ s(pgr), data = gbsg_prepared, arm = "arm", target_power = 0.95, L = 182, n_sim = 50, patience = 5, parallel.cores = 1 ) }) ``` ```{r echo=FALSE} ss_gam$results_plot + theme_bw(base_size = 14) ``` *** # Covariate-Dependent Censoring Models In many observational or longitudinal studies, the probability of being censored may depend on measured subject characteristics such as age, disease stage, or biomarker level, but **not** on the unobserved event time after conditioning on those covariates. This setting is referred to as **covariate-dependent (or conditionally independent) censoring**. Inverse-probability-of-censoring weighting (IPCW) can account for this dependence under the stated model and support valid estimation of the restricted mean survival time (RMST). ## Theory and Model We assume a single censoring mechanism described by a Cox proportional-hazards model for the censoring time: \[ \text{Pr}(C \le t \mid X) = 1 - G(t \mid X), \qquad G(t \mid X) = \exp\{-\Lambda_C(t \mid X)\}, \] where \(X\) denotes the observed covariates (e.g., age, sex, treatment arm). Each subject receives an inverse-probability weight \[ W_i = \frac{1}{\widehat{G}(Y_i \mid X_i)}, \qquad Y_i = \min(T_i, L), \] with \(T_i\) the event time and \(L\) the truncation horizon for RMST. Weights are stabilized and truncated to prevent numerical instability. The weighted regression \[ E[Y_i \mid A_i, X_i] = \beta_0 + \beta_A A_i + X_i^\top\beta_X \] is then fitted using weighted least squares. The variance of \(\hat\beta_A\) is estimated with a sandwich-type estimator that treats the censoring model as known. ### Power Calculation Formula Analytic power is computed as \[ \text{Power} = \Phi\!\left( \frac{|\hat\beta_A|}{\sigma_N} - z_{1-\alpha/2} \right), \qquad \sigma_N^2 = \frac{\widehat{\mathrm{Var}}(\hat\beta_A)}{N}, \] where \(\widehat{\mathrm{Var}}(\hat\beta_A)\) is the asymptotic variance estimated from the reference data and \(N\) is the total sample size. Because the same IPCW model is used for all subjects, no cause-specific components appear in the variance. ## Analytical Methods We demonstrate this approach using the `mgus2` dataset, modeling censoring as a function of age and sex under a single censoring mechanism. ```{r mgus2_data_prep, echo=FALSE} mgus_prepared <- mgus2 %>% mutate( status = ifelse(pstat == 1, 1, 0), # 1 = event, 0 = censored arm = ifelse(sex == "M", 1, 0) # treatment indicator ) %>% rename(time = futime) head(mgus_prepared) ``` ### Power Calculation Setting `dep_cens = TRUE` routes the call to the dependent-censoring adjusted estimator. ```{r dc_power_calc} dc_power_results <- vignette_cache("dc_power_calc", { rmst.power( Surv(time, status) ~ age, data = mgus_prepared, arm = "arm", sample_sizes = c(100, 250, 500), L = 120, dep_cens = TRUE, alpha = 0.05 ) }) ``` ```{r echo=FALSE} kbl(dc_power_results$results_summary, caption = "Analytic Power Analysis under Covariate-Dependent Censoring") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") dc_power_results$results_plot ``` ### Sample-Size Calculation We next estimate the per-arm sample size required to achieve 80% power at the same truncation time (10 years). ```{r mgus2_ss_calc} ss_dc_mgus <- vignette_cache("mgus2_ss_calc", { rmst.ss( Surv(time, status) ~ age, data = mgus_prepared, arm = "arm", target_power = 0.80, L = 120, dep_cens = TRUE, alpha = 0.05, n_start = 100, n_step = 50, max_n = 5000 ) }) ``` ```{r mgus2_table, echo=FALSE} kbl(ss_dc_mgus$results_summary, caption = "Estimated RMST Difference from Reference Data") %>% kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center") ss_dc_mgus$results_plot + theme_bw(base_size = 14) ``` --- # Interactive Shiny Application `RMSTpowerBoost` also includes a Shiny web application for users who prefer a graphical interface. ### Accessing the Application There are two ways to access the application: 1. **Live Web Version**: Access the application directly in your browser without any installation. - [Launch Web Application](https://arnab96.shinyapps.io/uthsc-app/) 2. **Run Locally from the R Package**: If you have installed the `RMSTpowerBoost-Package`, you can run the application on your own machine with the following command: ```{r, eval=FALSE} RMSTpowerBoost::run_app() ``` ### App Features - Upload a reference dataset or generate data inside the app. - Map input columns to the variables required for the analysis. - Select the RMST model, calculation method, and analysis parameters. - Review data summaries, Kaplan-Meier plots, result tables, and power curves. - Export the analysis as PDF or HTML. # Conclusion `RMSTpowerBoost` centers its interface on `rmst.power()` and `rmst.ss()` for RMST-based power and sample size calculations under linear, stratified, semiparametric, and covariate-dependent censoring models. Both functions use the familiar `Surv(time, status) ~ covariates` syntax, with model choice controlled by a small set of arguments (`strata`, `strata_type`, `dep_cens`, `type`). In practice, these procedures depend on representative reference data because effect sizes and variance components are estimated from that dataset. Bootstrap-based methods can also require substantial computation. Future development could extend simulation-based procedures for covariate-dependent censoring and add model diagnostic tools for assessing reference-data assumptions. # References