| Title: | Power and Sample Size for Restricted Mean Survival Time Based Clinical Trials |
|---|---|
| Description: | Tools for Restricted Mean Survival Time based study design and analysis planning. Provides power and sample size calculations for two-arm studies using direct modeling approaches from the literature, including stratified semiparametric additive models from Zhang and Schaubel (2024) <doi:10.1002/bimj.202200371>, linear Inverse Probability Weighting based models from Tian, Zhao and Wei (2014) <doi:10.1093/biostatistics/kxt050>, multiplicative stratified models from Wang (2019) <doi:10.1002/sim.8356>, and covariate-dependent censoring methods from Wang (2018) <doi:10.1007/s10985-017-9391-6>. |
| Authors: | Arnab Aich [aut, cre] (ORCID: <https://orcid.org/0009-0005-7801-6701>), Yuan Zhang [aut] (ORCID: <https://orcid.org/0000-0001-7626-9514>) |
| Maintainer: | Arnab Aich <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.3 |
| Built: | 2026-10-29 17:10:53 UTC |
| Source: | https://github.com/uthsc-zhang/rmstpowerboost-package |
Performs power analysis for a stratified, additive RMST model using the analytic variance estimator based on the method of Zhang & Schaubel (2024).
additive.power.analytical( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE )additive.power.analytical( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
strata_var |
A character string for the stratification variable. |
sample_sizes |
A numeric vector of sample sizes per stratum to calculate power for. |
linear_terms |
An optional character vector of other covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
verbose |
Logical; if |
This function computes power for the semiparametric additive RMST model
, where i indexes subjects and j
indexes strata.
The method uses Inverse Probability of Censoring Weighting (IPCW), with
weights derived from a stratified Cox model for
the censoring times fit on the original time scale. Here
if the event occurs before or follow-up reaches
. The regression coefficient is estimated by
centering the covariates and RMST values within each stratum and then solving
the resulting estimating equations in closed form.
The censoring model uses Breslow ties and is evaluated at each subject's
truncated outcome. The weights are used as fitted, without capping or
truncation: Equations (11)-(12) of Zhang & Schaubel (2024) derive the
estimator and its simplified sandwich variance under the fitted weights, so
capping would break the correspondence with the published formulas, and it
compresses exactly the upper tail that IPCW exists to up-weight. The weight
distribution is reported in
model_output$censoring_weights$raw_summary.
Power is obtained from the asymptotic sandwich variance of .
This implementation uses the robust variance estimator
, where and are empirical
estimates of the variance components.
A list containing:
results_data |
A |
results_plot |
A |
Zhang, Y. and Schaubel, D. E. (2024). Semiparametric Additive Modeling of the Restricted Mean Survival Time. Biometrical Journal, 66:e202200371. doi:10.1002/bimj.202200371
set.seed(123) pilot_df_strat <- data.frame( time = rexp(150, 0.1), status = rbinom(150, 1, 0.8), arm = rep(0:1, each = 75), region = factor(rep(c("A", "B", "C"), each = 50)), age = rnorm(150, 60, 10) ) # Introduce an additive treatment effect pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] + 1.5 power_results <- additive.power.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(100, 150, 200), linear_terms = "age", L = 12 ) print(power_results$results_data) print(power_results$results_plot)set.seed(123) pilot_df_strat <- data.frame( time = rexp(150, 0.1), status = rbinom(150, 1, 0.8), arm = rep(0:1, each = 75), region = factor(rep(c("A", "B", "C"), each = 50)), age = rnorm(150, 60, 10) ) # Introduce an additive treatment effect pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] + 1.5 power_results <- additive.power.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(100, 150, 200), linear_terms = "age", L = 12 ) print(power_results$results_data) print(power_results$results_plot)
Calculates the required sample size for a target power using the analytic method for a stratified, additive RMST model.
additive.ss.analytical( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )additive.ss.analytical( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
strata_var |
A character string for the stratification variable. |
target_power |
A single numeric value for the desired power. |
linear_terms |
An optional character vector of other covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
n_start |
The starting sample size per stratum for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per stratum to search up to. |
verbose |
Logical; if |
This function estimates the additive treatment effect and its asymptotic
variance once from the pilot data, then increases the per-stratum sample
size until the target power is reached or the search limit is hit. It uses
the same stratum-centering framework as additive.power.analytical,
including the uncapped IPCW weights described there.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
Zhang, Y. and Schaubel, D. E. (2024). Semiparametric Additive Modeling of the Restricted Mean Survival Time. Biometrical Journal, 66:e202200371. doi:10.1002/bimj.202200371
set.seed(123) pilot_df_strat <- data.frame( time = rexp(150, 0.1), status = rbinom(150, 1, 0.8), arm = rep(0:1, each = 75), region = factor(rep(c("A", "B", "C"), each = 50)), age = rnorm(150, 60, 10) ) # Introduce an additive treatment effect pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] + 1.5 # Find the required sample size per stratum for 80% power ss_results <- additive.ss.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.50, L = 18, # n_start = 200, n_step = 50, max_n_per_arm = 1000 ) print(ss_results$results_data) print(ss_results$results_plot)set.seed(123) pilot_df_strat <- data.frame( time = rexp(150, 0.1), status = rbinom(150, 1, 0.8), arm = rep(0:1, each = 75), region = factor(rep(c("A", "B", "C"), each = 50)), age = rnorm(150, 60, 10) ) # Introduce an additive treatment effect pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] + 1.5 # Find the required sample size per stratum for 80% power ss_results <- additive.ss.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.50, L = 18, # n_start = 200, n_step = 50, max_n_per_arm = 1000 ) print(ss_results$results_data) print(ss_results$results_plot)
Columns:
time: observed time
status: event indicator (1 = event, 0 = censored)
arm: treatment arm (0 = control, 1 = treatment)
age: baseline age
gender: factor with two levels
A data frame with 150 rows and 5 variables.
data(aft_lognormal_L12_n150) table(aft_lognormal_L12_n150$arm)data(aft_lognormal_L12_n150) table(aft_lognormal_L12_n150$arm)
Columns: time, status, arm, age, gender.
A data frame with 200 rows and 5 variables.
data(aft_weibull_L24_n200) mean(aft_weibull_L24_n200$status == 0)data(aft_weibull_L24_n200) mean(aft_weibull_L24_n200$status == 0)
Builds a binary (Bernoulli) covariate definition list.
covar_binary(name, p = 0.5)covar_binary(name, p = 0.5)
name |
Column name in the generated data. |
p |
Probability of value 1. Default |
A named list suitable as one element of the covariates argument.
covar_continuous, covar_categorical,
rmst.sim
covar_binary("female", p = 0.52)covar_binary("female", p = 0.52)
Builds a multi-level categorical covariate definition list.
covar_categorical(name, probs, labels = NULL)covar_categorical(name, probs, labels = NULL)
name |
Column name in the generated data. |
probs |
Numeric vector of category probabilities (must sum to 1). |
labels |
Optional character vector of level labels. Length must equal
|
A named list suitable as one element of the covariates argument.
covar_continuous, covar_binary,
rmst.sim
covar_categorical("stage", probs = c(0.4, 0.35, 0.25), labels = c("I", "II", "III"))covar_categorical("stage", probs = c(0.4, 0.35, 0.25), labels = c("I", "II", "III"))
Builds a covariate definition list for use in rmst.sim,
recipe_quick_aft, or recipe_quick_ph.
covar_continuous(name, dist = "normal", ...)covar_continuous(name, dist = "normal", ...)
name |
Column name in the generated data. |
dist |
Distribution name. One of |
... |
Named distribution parameters passed as the |
A named list suitable as one element of the covariates argument.
covar_binary, covar_categorical,
rmst.sim
covar_continuous("age", dist = "normal", mean = 50, sd = 10) covar_continuous("bmi", dist = "lognormal", meanlog = 3.2, sdlog = 0.2)covar_continuous("age", dist = "normal", mean = 50, sd = 10) covar_continuous("bmi", dist = "lognormal", meanlog = 3.2, sdlog = 0.2)
Performs power analysis for an RMST model when the censoring mechanism depends on observed covariates.
DC.power.analytical( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("wls", "sandwich") )DC.power.analytical( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("wls", "sandwich") )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
sample_sizes |
A numeric vector of sample sizes per arm to calculate power for. |
linear_terms |
Optional character vector of additional covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
verbose |
Logical; if |
test |
Which test the reported power refers to. |
This function assumes a single censoring process whose hazard depends on covariates (but not on treatment by default). It fits a Cox model for censoring
using Surv(time, status==0) ~ linear\_terms, then forms inverse-probability-of-censoring
weights (IPCW) evaluated at
, where if the event is observed before
or follow-up reaches . The censoring model uses Breslow ties, is fit
on the original follow-up, and is evaluated at each subject's truncated outcome,
with the uncentered baseline hazard paired with an uncentered linear predictor.
The weights are used as fitted, without capping, truncation, or a floor on
, so that the estimator and its variance reproduce the published
formulas exactly; their distribution is reported in
model_output$censoring_weights$raw_summary. With no observed censoring the
weights are one and no censoring model is fit.
The RMST regression is then fit by weighted least squares,
and the sampling variance of the treatment effect is the sandwich
with the fitted weights treated as fixed, so it ignores
uncertainty from estimating .
The test argument selects which hypothesis test the reported power refers to,
and this matters because the two candidates do not use the same variance. The
sandwich above is consistent for the sampling variability of the estimator. The
ordinary weighted-least-squares standard error reported by summary(lm(...)) is
not, because IPCW weights are sampling weights rather than inverse-variance
weights, so it typically runs above the true sampling standard deviation of the
estimate, which pushes the rejection threshold out and costs power.
With test = "wls" (the default) the power refers to the weighted least-squares
t-test that a standard lm(..., weights = w) analysis reports, as in
linear.power.analytical(). The sampling distribution of the estimate is taken
from the sandwich variance and the rejection threshold from the
WLS standard error , giving
and model_output$coefficient_table is the t-table that summary()
reports for the weighted fit.
With test = "sandwich" the power refers to the sandwich Wald test, in which one
variance sets both the sampling distribution and the threshold, giving the usual
. That is the more efficient test,
and model_output$coefficient_table then uses the robust variance. Both
standard errors are always returned in model_output$variance_components.
Note: This implementation models a single censoring process and does not handle competing risks.
A list with:
results_data |
A data.frame with |
results_plot |
A ggplot object of the power curve. |
results_summary |
A data.frame summarizing the pilot arm effect. |
Iteratively finds required per-arm sample size for a target power,
using the same IPCW-based analytic variance as DC.power.analytical.
DC.ss.analytical( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("wls", "sandwich") )DC.ss.analytical( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("wls", "sandwich") )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
target_power |
A single numeric value for the desired power. |
linear_terms |
Optional character vector of additional covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
n_start |
The starting sample size per arm for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per arm to search. |
verbose |
Logical; if |
test |
Which test the required sample size refers to. |
Uses a single censoring Cox model Surv(time, status==0) ~ linear_terms to form
IPCW weights and fits a weighted RMST regression
among subjects whose truncated RMST outcome is observed. Treatment is excluded from
the censoring model by default. Competing risks are not modeled. Variance ignores
uncertainty in . It uses the same uncapped IPCW weights, sandwich
sampling variance, and test argument as DC.power.analytical.
Note: This implementation models a single censoring process and does not handle competing risks.
A list with:
results_data |
data.frame with |
results_plot |
ggplot showing the search path. |
results_summary |
data.frame summarizing the pilot arm effect. |
Given one element returned by load_recipe_sets() (a list with data and
a rich meta list), this builds tidy tables describing the data-generating
process: model family, baseline parameters, linear predictor coefficients
(intercept / treatment / covariates), treatment assignment, and censoring.
describe_generation(set)describe_generation(set)
set |
A single element from |
A list of data frames:
header: n, L (analysis horizon; stored as attr "tau"), model,
event_rate, achieved_censoring.
baseline: flattened baseline parameters.
effects: coefficients for intercept/treatment and covariates (and formula terms if used).
treatment: assignment mechanism and knobs.
censoring: censoring mode/target/admin time.
covariates: each generated covariate with its distribution and parameters.
files: paths to on-disk files (csv/rds/rdata) when available.
covs <- list(covar_continuous("x", mean = 0, sd = 1)) rec <- recipe_quick_aft( n = 20, model = "aft_lognormal", baseline = list(mu = 2.7, sigma = 0.6), treat_effect = -0.2, covariates = covs, target_censoring = 0.2, seed = 123 ) out <- file.path(tempdir(), "rmst_describe_generation") generate_recipe_sets(rec, out_dir = out, formats = "rds", n_reps = 1) sets <- load_recipe_sets(file.path(out, "manifest.rds")) spec <- describe_generation(sets[[1]]) spec$header spec$baseline spec$effects spec$treatment spec$censoring spec$covariates spec$filescovs <- list(covar_continuous("x", mean = 0, sd = 1)) rec <- recipe_quick_aft( n = 20, model = "aft_lognormal", baseline = list(mu = 2.7, sigma = 0.6), treat_effect = -0.2, covariates = covs, target_censoring = 0.2, seed = 123 ) out <- file.path(tempdir(), "rmst_describe_generation") generate_recipe_sets(rec, out_dir = out, formats = "rds", n_reps = 1) sets <- load_recipe_sets(file.path(out, "manifest.rds")) spec <- describe_generation(sets[[1]]) spec$header spec$baseline spec$effects spec$treatment spec$censoring spec$covariates spec$files
Performs a power analysis for given sample sizes using a flexible, semiparametric additive model for the RMST based on pseudo-observations.
GAM.power.boot( pilot_data, time_var, status_var, arm_var, strata_var = NULL, sample_sizes, linear_terms = NULL, smooth_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, verbose = FALSE )GAM.power.boot( pilot_data, time_var, status_var, arm_var, strata_var = NULL, sample_sizes, linear_terms = NULL, smooth_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
strata_var |
An optional character string for a stratification variable. |
sample_sizes |
A numeric vector of sample sizes per arm/stratum. |
linear_terms |
Optional character vector of covariates with a linear effect. |
smooth_terms |
Optional character vector of covariates with a non-linear effect. |
L |
The numeric truncation time for RMST. |
n_sim |
Number of bootstrap simulations. |
alpha |
The significance level. |
parallel.cores |
Number of cores for parallel processing. |
verbose |
Logical; if |
This function estimates power by bootstrap resampling the pilot data,
computing RMST pseudo-observations, fitting a semiparametric additive model,
and recording whether the treatment effect is significant. If strata_var
is supplied, resampling is carried out within strata and the model includes
stratum-specific intercepts with a common treatment effect (the same
single-effect hypothesis tested by the analytical stratified models). The
stratification variable is treated as a factor, so a numeric stratum code
still gives one intercept per stratum rather than a single slope. The model
formula can include linear terms and non-linear smooth terms (s()).
The jackknife pseudo-observations are computed from the Kaplan-Meier curve
within each arm (within each stratum when strata_var is supplied), with the
area under the curve up to L computed exactly as a sum of rectangles.
Power is estimated as the proportion of simulations with a treatment-effect
p-value below alpha. The treatment effect is the coefficient of arm_var
coded as a 0/1 factor, selected by its exact name, so a covariate whose name
begins with the arm variable's name is never mistaken for it. This
pseudo-observation approach models RMST directly without using the IPCW
estimating equations used elsewhere in the package.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
status_var should be 1 for an event, 0 for censored. arm_var
should be 1 for treatment, 0 for control.
pilot_df <- data.frame( time = rexp(100, 0.08), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 60, 10) ) # Add a treatment effect pilot_df$time[pilot_df$arm == 1] <- pilot_df$time[pilot_df$arm == 1] * 1.3 power_results <- GAM.power.boot( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", sample_sizes = c(100, 150), linear_terms = "age", L = 15, n_sim = 100, # Use more sims in practice, e.g., 1000 parallel.cores = 2 ) print(power_results$results_data) print(power_results$results_plot)pilot_df <- data.frame( time = rexp(100, 0.08), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 60, 10) ) # Add a treatment effect pilot_df$time[pilot_df$arm == 1] <- pilot_df$time[pilot_df$arm == 1] * 1.3 power_results <- GAM.power.boot( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", sample_sizes = c(100, 150), linear_terms = "age", L = 15, n_sim = 100, # Use more sims in practice, e.g., 1000 parallel.cores = 2 ) print(power_results$results_data) print(power_results$results_plot)
Performs an iterative sample size search to achieve a target power using a flexible, semiparametric additive model for the RMST.
GAM.ss.boot( pilot_data, time_var, status_var, arm_var, strata_var = NULL, target_power, linear_terms = NULL, smooth_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )GAM.ss.boot( pilot_data, time_var, status_var, arm_var, strata_var = NULL, target_power, linear_terms = NULL, smooth_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
strata_var |
An optional string for a stratification variable. |
target_power |
A single numeric value for the target power. |
linear_terms |
Optional character vector of covariates with a linear effect. |
smooth_terms |
Optional character vector of covariates with a non-linear effect. |
L |
The numeric truncation time for RMST. |
n_sim |
Number of bootstrap simulations per search step. |
alpha |
The significance level. |
parallel.cores |
Number of cores for parallel processing. |
patience |
Number of consecutive non-improving steps in the search before terminating. |
n_start |
The starting sample size per arm/stratum for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per arm/stratum to search up to. |
verbose |
Logical; if |
This function increases the sample size step by step until the estimated
power reaches target_power or a stopping rule is met. At each step it runs
the same bootstrap procedure described in GAM.power.boot, including the
factor treatment of strata_var, the exact pseudo-observations, and the
exact selection of the treatment coefficient.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
This function's methodology is bootstrap-based.
pilot_df_effect <- data.frame( time = c(stats::rexp(50, 0.1), stats::rexp(50, 0.04)), # Effect status = stats::rbinom(100, 1, 0.9), arm = rep(0:1, each = 50) ) ss_results <- GAM.ss.boot( pilot_data = pilot_df_effect, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 15, n_sim = 100, # Low n_sim for example n_start = 100, n_step = 50, patience = 2, parallel.cores = 2 ) print(ss_results$results_data) print(ss_results$results_plot)pilot_df_effect <- data.frame( time = c(stats::rexp(50, 0.1), stats::rexp(50, 0.04)), # Effect status = stats::rbinom(100, 1, 0.9), arm = rep(0:1, each = 50) ) ss_results <- GAM.ss.boot( pilot_data = pilot_df_effect, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 15, n_sim = 100, # Low n_sim for example n_start = 100, n_step = 50, patience = 2, parallel.cores = 2 ) print(ss_results$results_data) print(ss_results$results_plot)
Generate covariate matrix/data frame from a recipe
gen_covariates(n, covariates)gen_covariates(n, covariates)
n |
sample size |
covariates |
list(defs = list(...)) |
data.frame of covariates
defs <- list( list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1)), list(name="z", type="categorical", dist="categorical", params=list(prob=c(0.3,0.7), labels=c("A","B"))) ) X <- gen_covariates(10, list(defs = defs))defs <- list( list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1)), list(name="z", type="categorical", dist="categorical", params=list(prob=c(0.3,0.7), labels=c("A","B"))) ) X <- gen_covariates(10, list(defs = defs))
Builds a grid of scenarios from a base recipe and a named list of variations,
simulates one or more replicates per scenario, and writes the datasets to files
in a target folder as .txt (tab), .csv, .rds, and/or .RData.
A manifest is written as manifest.rds. The recipe specification does not
include an analysis horizon; specify L later in downstream analysis functions.
generate_recipe_sets( base_recipe, vary = list(), out_dir, formats = c("txt", "csv", "rds", "rdata"), n_reps = 1L, seed_base = NULL, filename_template = "sc{scenario_id}_r{rep}" )generate_recipe_sets( base_recipe, vary = list(), out_dir, formats = c("txt", "csv", "rds", "rdata"), n_reps = 1L, seed_base = NULL, filename_template = "sc{scenario_id}_r{rep}" )
base_recipe |
A recipe list (use |
vary |
Named list; keys are dotted paths inside the recipe (e.g., |
out_dir |
Directory to write datasets and |
formats |
Character vector subset of |
n_reps |
Integer; number of replicates per scenario. |
seed_base |
Optional integer retained in the manifest as per-replicate
provenance metadata computed as |
filename_template |
Base filename (no extension) with placeholders:
|
Invisibly returns the manifest data.frame and writes manifest.rds.
covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) rec <- recipe_quick_aft(60, "aft_lognormal", baseline=list(mu=2.7, sigma=0.6), treat_effect=-0.2, covariates=covs, target_censoring=0.2) out <- file.path(tempdir(), "rmst_checks") dir.create(out, showWarnings = FALSE, recursive = TRUE) man <- generate_recipe_sets(rec, vary=list(n=c(60,80)), out_dir=out, formats=c("csv","rds"), n_reps=1, seed_base=123)covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) rec <- recipe_quick_aft(60, "aft_lognormal", baseline=list(mu=2.7, sigma=0.6), treat_effect=-0.2, covariates=covs, target_censoring=0.2) out <- file.path(tempdir(), "rmst_checks") dir.create(out, showWarnings = FALSE, recursive = TRUE) man <- generate_recipe_sets(rec, vary=list(n=c(60,80)), out_dir=out, formats=c("csv","rds"), n_reps=1, seed_base=123)
Performs power analysis using a direct formula based on the asymptotic variance estimator for the linear RMST model.
linear.power.analytical( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("wls", "sandwich") )linear.power.analytical( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("wls", "sandwich") )
pilot_data |
A |
time_var |
A character string specifying the name of the time-to-event variable. |
status_var |
A character string specifying the name of the event status variable (1=event, 0=censored). |
arm_var |
A character string specifying the name of the treatment arm variable (1=treatment, 0=control). |
sample_sizes |
A numeric vector of sample sizes per arm to calculate power for. |
linear_terms |
An optional character vector of other covariate names to include in the model. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level for the power calculation (Type I error rate). |
verbose |
Logical; if |
test |
Which test the reported power refers to. |
This function implements the analytic power calculation for the direct linear regression model of the Restricted Mean Survival Time (RMST) proposed by Tian et al. (2014). The core of the method is a weighted linear model of the form
where is the event time truncated at , is the
treatment indicator, and is the treatment effect of interest.
To handle right-censoring, the method uses Inverse Probability of Censoring
Weighting (IPCW). Let and
if the event occurs before or follow-up reaches
. The weight is , where
estimates censoring survival conditional on treatment and
linear_terms, using a Cox proportional hazards model with Breslow
ties fit on the original time scale and evaluated at each subject's
truncated outcome. The weights are used as fitted, without capping or
truncation, so that the estimator and its variance reproduce the published
formulas exactly; their distribution is reported in
model_output$censoring_weights$raw_summary.
Power is calculated analytically based on the asymptotic properties of the
coefficient estimators. The variance of the treatment effect estimator, , is derived from a
robust sandwich variance estimator of the form . This is the
simplified sandwich variance of Zhang and Schaubel (2024, Equations 11-12), which
treats the fitted IPCW weights as fixed: A is the weighted information matrix
, with the weights entering linearly, and B is the empirical second
moment of the weighted score contributions,
with .
The intercept-plus-slopes parameterization used here is, for a single stratum,
algebraically equivalent to the weighted-centered-covariate form without intercept
used in the paper. As in the paper's simulations and data analysis, the additional
influence-function terms for estimation of the Cox censoring model given in their
Theorem 3.1 are not included.
The test argument selects which hypothesis test the reported power refers to, and
this matters because the two candidates do not use the same variance. The sandwich
above is consistent for the sampling variability of the estimator. The ordinary
weighted-least-squares standard error reported by summary(lm(...)), which is the
test that linear.power.boot() simulates, is not, because IPCW weights are sampling
weights rather than inverse-variance weights. In a representative pilot the WLS
standard error runs about 20 percent above the true sampling standard deviation of
the estimate, which pushes the rejection threshold out and costs power.
With test = "wls" (the default) the power refers to the weighted least-squares
t-test. The sampling distribution of the estimate is taken from the sandwich
variance and the rejection threshold from the WLS standard error
, giving
This is the analysis that linear.power.boot() performs, so the analytic and
bootstrap engines then describe the same test.
With test = "sandwich" the power refers to the sandwich Wald test, in which one
variance sets both the sampling distribution and the threshold, giving the usual
. That is the more efficient test and
the one matching the paper's inference, but it is not what the bootstrap simulates.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
pilot_df <- data.frame( time = rexp(100, 0.1), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 55, 10) ) power_results <- linear.power.analytical( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", linear_terms = "age", sample_sizes = c(100, 200, 300), L = 10 ) print(power_results$results_data) print(power_results$results_plot)pilot_df <- data.frame( time = rexp(100, 0.1), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 55, 10) ) power_results <- linear.power.analytical( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", linear_terms = "age", sample_sizes = c(100, 200, 300), L = 10 ) print(power_results$results_data) print(power_results$results_plot)
Performs a power analysis for given sample sizes based on the direct linear regression model for RMST, using a bootstrap simulation approach.
linear.power.boot( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, verbose = FALSE )linear.power.boot( pilot_data, time_var, status_var, arm_var, sample_sizes, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
sample_sizes |
A numeric vector of sample sizes per arm to calculate power for. |
linear_terms |
Optional character vector of other covariates for the linear model. |
L |
The numeric truncation time for RMST. |
n_sim |
The number of bootstrap simulations to run for each sample size. |
alpha |
The significance level (Type I error rate). |
verbose |
Logical; if |
This function estimates power by generating a number of bootstrap
samples (n_sim) from the provided pilot data by resampling with replacement.
For each bootstrap sample, it performs the following steps:
Estimates conditional censoring survival using a Cox model with treatment
and linear_terms as predictors (survival::coxph, Breslow ties).
Calculates Inverse Probability of Censoring Weights (IPCW) for each
observation using , where
if the event occurs before or follow-up
reaches . Here is conditional on those predictors.
The weights are used as fitted, without capping, with the censoring model
refitted in each bootstrap sample.
Fits a weighted linear model (stats::lm) to the truncated RMST outcome
among subjects with observed .
Extracts the p-value for the treatment arm_var coefficient.
The final power for a given sample size is the proportion of the n_sim simulations
where this p-value is less than the significance level alpha. This simulation-based
approach can be useful when analytic approximations are less reliable, but it can
be computationally intensive.
The p-values come from the weighted least-squares fit with model-based
standard errors, so the estimated power corresponds to a trial analyzed the
same way. Because model-based standard errors treat IPCW weights as
precision weights, this test is somewhat conservative relative to the robust
sandwich test underlying linear.power.analytical; bootstrap power can
therefore be lower than analytic power for the same design.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
status_var should be 1 for an event, 0 for censored. arm_var
should be 1 for treatment, 0 for control.
pilot_df <- data.frame( time = rexp(100, 0.1), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 60, 8) ) # Introduce a treatment effect for a more interesting example pilot_df$time[pilot_df$arm == 1] <- pilot_df$time[pilot_df$arm == 1] * 1.5 power_results <- linear.power.boot( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", linear_terms = "age", sample_sizes = c(100, 150, 200), L = 10, n_sim = 200 # Use more simulations in practice (e.g., 1000) ) print(power_results$results_data) print(power_results$results_plot)pilot_df <- data.frame( time = rexp(100, 0.1), status = rbinom(100, 1, 0.7), arm = rep(0:1, each = 50), age = rnorm(100, 60, 8) ) # Introduce a treatment effect for a more interesting example pilot_df$time[pilot_df$arm == 1] <- pilot_df$time[pilot_df$arm == 1] * 1.5 power_results <- linear.power.boot( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", linear_terms = "age", sample_sizes = c(100, 150, 200), L = 10, n_sim = 200 # Use more simulations in practice (e.g., 1000) ) print(power_results$results_data) print(power_results$results_plot)
Calculates the required sample size for a target power using an analytic formula based on the methods of Tian et al. (2014).
linear.ss.analytical( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("wls", "sandwich") )linear.ss.analytical( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("wls", "sandwich") )
pilot_data |
A |
time_var |
A character string specifying the name of the time-to-event variable. |
status_var |
A character string specifying the name of the event status variable (1=event, 0=censored). |
arm_var |
A character string specifying the name of the treatment arm variable (1=treatment, 0=control). |
target_power |
A single numeric value for the desired power (e.g., 0.80 or 0.90). |
linear_terms |
An optional character vector of other covariate names to include in the model. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
n_start |
The starting sample size per arm for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per arm to search up to. |
verbose |
Logical; if |
test |
Which test the required sample size refers to. |
This function performs an iterative search to find the sample size needed to
achieve a specified target_power. It uses the same underlying theory as
linear.power.analytical, including the uncapped IPCW weights, the
simplified sandwich variance, and the test argument described there. First, it estimates the treatment effect size and its
asymptotic variance from the pilot data. Then, it iteratively calculates the
power for increasing sample sizes using the analytic formula until the
target power is achieved.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
pilot_df <- data.frame( time = c(rexp(50, 0.1), rexp(50, 0.07)), # Introduce an effect status = rbinom(100, 1, 0.8), arm = rep(0:1, each = 50), age = rnorm(100, 55, 10) ) ss_results <- linear.ss.analytical( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 10 ) print(ss_results$results_data) print(ss_results$results_plot)pilot_df <- data.frame( time = c(rexp(50, 0.1), rexp(50, 0.07)), # Introduce an effect status = rbinom(100, 1, 0.8), arm = rep(0:1, each = 50), age = rnorm(100, 55, 10) ) ss_results <- linear.ss.analytical( pilot_data = pilot_df, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 10 ) print(ss_results$results_data) print(ss_results$results_plot)
Performs an iterative sample size search to achieve a target power based on the direct linear regression model for RMST, using bootstrap simulation.
linear.ss.boot( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )linear.ss.boot( pilot_data, time_var, status_var, arm_var, target_power, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
target_power |
A single numeric value for the target power (e.g., 0.80). |
linear_terms |
Optional character vector of other covariates for the linear model. |
L |
The numeric truncation time for RMST. |
n_sim |
The number of bootstrap simulations per search step. |
alpha |
The significance level. |
patience |
The number of consecutive non-improving steps in the search before terminating. |
n_start |
The starting sample size per arm for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per arm to search up to. |
verbose |
Logical; if |
This function iteratively searches for the required sample size to achieve
a specified target_power. At each step of the search, it runs a full bootstrap
simulation (n_sim iterations), as described in linear.power.boot, to
estimate the power for the current sample size. The search stops when the
target power is achieved or other stopping criteria (e.g., patience) are met.
Due to the nested simulation structure, this function can be very time-consuming.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
status_var should be 1 for an event, 0 for censored. arm_var
should be 1 for treatment, 0 for control.
pilot_df_effect <- data.frame( time = c(rexp(50, 0.1), rexp(50, 0.05)), # Effect present status = rbinom(100, 1, 0.8), arm = rep(0:1, each = 50) ) ss_results <- linear.ss.boot( pilot_data = pilot_df_effect, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 10, n_sim = 200, # Low n_sim for example patience = 2, n_start = 100, n_step = 50, max_n_per_arm = 500 ) print(ss_results$results_data) print(ss_results$results_plot)pilot_df_effect <- data.frame( time = c(rexp(50, 0.1), rexp(50, 0.05)), # Effect present status = rbinom(100, 1, 0.8), arm = rep(0:1, each = 50) ) ss_results <- linear.ss.boot( pilot_data = pilot_df_effect, time_var = "time", status_var = "status", arm_var = "arm", target_power = 0.80, L = 10, n_sim = 200, # Low n_sim for example patience = 2, n_start = 100, n_step = 50, max_n_per_arm = 500 ) print(ss_results$results_data) print(ss_results$results_plot)
Reads a manifest.rds created by generate_recipe_sets(), loads one
dataset per row (preferring rds to rdata to csv to txt),
restores attribute "achieved_censoring", and returns a named list
of list(data = <data.frame>, meta = <list>).
load_recipe_sets(manifest_path)load_recipe_sets(manifest_path)
manifest_path |
Path to |
A named list where each element is list(data=..., meta=...).
covs <- list(covar_continuous("x", mean = 0, sd = 1)) rec <- recipe_quick_aft( n = 20, model = "aft_lognormal", baseline = list(mu = 2.7, sigma = 0.6), treat_effect = -0.2, covariates = covs, target_censoring = 0.2, seed = 123 ) out <- file.path(tempdir(), "rmst_load_recipe_sets") generate_recipe_sets(rec, out_dir = out, formats = "rds", n_reps = 1) sets <- load_recipe_sets(file.path(out, "manifest.rds")) names(sets) str(sets[[1]]$meta)covs <- list(covar_continuous("x", mean = 0, sd = 1)) rec <- recipe_quick_aft( n = 20, model = "aft_lognormal", baseline = list(mu = 2.7, sigma = 0.6), treat_effect = -0.2, covariates = covs, target_censoring = 0.2, seed = 123 ) out <- file.path(tempdir(), "rmst_load_recipe_sets") generate_recipe_sets(rec, out_dir = out, formats = "rds", n_reps = 1) sets <- load_recipe_sets(file.path(out, "manifest.rds")) names(sets) str(sets[[1]]$meta)
Performs power analysis for a multiplicative, stratified RMST model using an analytic method based on the work of Wang et al. (2019).
MS.power.analytical( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("sandwich", "wls") )MS.power.analytical( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, alpha = 0.05, verbose = FALSE, test = c("sandwich", "wls") )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
strata_var |
A character string for the stratification variable. |
sample_sizes |
A numeric vector of sample sizes per stratum to calculate power for. |
linear_terms |
An optional character vector of other covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
verbose |
Logical; if |
test |
Which test the reported power refers to. |
This function is based on the method for evaluating center (stratum) effects
using a multiplicative model for RMST: .
The method uses IPCW with a stratified Cox model for the censoring distribution
fit on the original time scale and weights of the form
, where if the event occurs
before or follow-up reaches . The stratified Cox model uses
Breslow ties and is evaluated at each subject's truncated outcome. The
weights are used as fitted, without capping or truncation, so that the
estimator and its variance reproduce the published formulas exactly; their
distribution is reported in
model_output$censoring_weights$raw_summary.
Formal estimation of requires an iterative solver for the estimating
equation given in Equation (8) of Wang et al. (2019). Because this is computationally
intensive, this implementation uses a log-linear working model fitted by weighted
least squares to pseudo-observations (lm(log(Y_rmst) ~ ...)) as a tractable
approximation. The approximation is consistent when the log-linear mean structure
is well-specified but may differ from the formal estimator under strong misspecification.
The power calculation relies on the asymptotic variance of the log-RMST ratio
estimator, . The variance is the robust sandwich estimator
of the weighted least-squares estimating
equations, with and the empirical second moment
of the weighted residuals, mirroring the sandwich form described in Theorem 1
of Wang et al. (2019).
The test argument selects which hypothesis test the reported power refers to,
and this matters because the two candidates do not use the same variance. The
sandwich above is consistent for the sampling variability of the estimator. The
ordinary weighted-least-squares standard error reported by summary(lm(...)) is
not, because IPCW weights are sampling weights rather than inverse-variance
weights, so it typically runs above the true sampling standard deviation of the
estimate, which pushes the rejection threshold out and costs power.
With test = "sandwich" (the default) the power refers to the sandwich Wald test,
in which one variance sets both the sampling distribution and the threshold,
giving the usual . This is the
inference of Wang et al. (2019), and coefficient standard errors in
model_output$coefficient_table then use the robust variance.
With test = "wls" the power refers to the weighted least-squares t-test: the
sampling distribution of the estimate is taken from the sandwich variance
and the rejection threshold from the WLS standard error
, giving
and model_output$coefficient_table becomes the t-table that
summary() reports for the weighted fit. Both standard errors are always
returned in model_output$variance_components.
Note that the default differs from linear.power.analytical(), which defaults to
test = "wls" because linear.power.boot() simulates exactly that test. No
bootstrap engine in this package simulates the IPCW-weighted t-test for a
stratified model: MS.power.boot() fits an unweighted model to jackknife
pseudo-observations. "sandwich" is therefore the default here.
A list containing:
results_data |
A |
results_plot |
A |
set.seed(123) pilot_df_strat <- data.frame( time = rexp(120, 0.15), status = rbinom(120, 1, 0.6), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] * 1.5 power_results <- MS.power.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(50, 75, 100), L = 10, alpha = 0.05 ) print(power_results$results_data)set.seed(123) pilot_df_strat <- data.frame( time = rexp(120, 0.15), status = rbinom(120, 1, 0.6), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] * 1.5 power_results <- MS.power.analytical( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(50, 75, 100), L = 10, alpha = 0.05 ) print(power_results$results_data)
Performs power analysis based on a multiplicative model for RMST for stratified trials, using a bootstrap simulation approach.
MS.power.boot( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, verbose = FALSE )MS.power.boot( pilot_data, time_var, status_var, arm_var, strata_var, sample_sizes, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
strata_var |
A character string for the stratification variable. |
sample_sizes |
A numeric vector of sample sizes per stratum to calculate power for. |
linear_terms |
Optional character vector of covariates for the model. |
L |
The numeric truncation time for RMST. |
n_sim |
Number of bootstrap simulations. |
alpha |
The significance level. |
parallel.cores |
Number of cores for parallel processing. |
verbose |
Logical; if |
This function estimates power through bootstrap simulation, resampling
within each stratum defined by strata_var. In each of the n_sim iterations:
Jackknife pseudo-observations for the RMST are calculated.
A log-linear model (stats::lm) is fitted to the log(pseudo_obs).
This models the multiplicative relationship on the original RMST scale,
i.e., . The model formula
includes stratum-specific intercepts and a common treatment effect,
matching the single-effect hypothesis tested by MS.power.analytical.
The p-value for the treatment effect is extracted from the model summary.
Note that pseudo-observations that are non-positive are dropped before
taking logs; when many pseudo-values are non-positive (e.g., heavy
truncation or very small L), this selection can bias the estimate.
Power is determined as the proportion of simulations where the p-value is
less than alpha.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
status_var should be 1/0. arm_var should be 1/0. strata_var
is a mandatory argument.
pilot_df_strat <- data.frame( time = rexp(120, 0.15), status = rbinom(120, 1, 0.6), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] * 1.4 power_results <- MS.power.boot( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(50, 75), L = 10, n_sim = 100 # Low n_sim for example ) print(power_results$results_data)pilot_df_strat <- data.frame( time = rexp(120, 0.15), status = rbinom(120, 1, 0.6), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) pilot_df_strat$time[pilot_df_strat$arm == 1] <- pilot_df_strat$time[pilot_df_strat$arm == 1] * 1.4 power_results <- MS.power.boot( pilot_data = pilot_df_strat, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", sample_sizes = c(50, 75), L = 10, n_sim = 100 # Low n_sim for example ) print(power_results$results_data)
Calculates the required sample size for a target power using the analytic (approximate) method from Wang et al. (2019).
MS.ss.analytical( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("sandwich", "wls") )MS.ss.analytical( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, alpha = 0.05, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE, test = c("sandwich", "wls") )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable (1=event, 0=censored). |
arm_var |
A character string for the treatment arm variable (1=treatment, 0=control). |
strata_var |
A character string for the stratification variable. |
target_power |
A single numeric value for the desired power. |
linear_terms |
An optional character vector of other covariate names. |
L |
The numeric value for the RMST truncation time. |
alpha |
The significance level (Type I error rate). |
n_start |
The starting sample size per stratum for the search. |
n_step |
The increment in sample size at each step of the search. |
max_n_per_arm |
The maximum sample size per stratum to search up to. |
verbose |
Logical; if |
test |
Which test the required sample size refers to. |
This function estimates the log-RMST ratio and its asymptotic variance once
from the pilot data, then increases the per-stratum sample size until the
target power is reached or the search limit is hit. It uses the same
log-linear approximation and robust sandwich variance as
MS.power.analytical, including the uncapped IPCW weights and the test
argument described there.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
set.seed(456) pilot_df_strat_effect <- data.frame( time = c(rexp(60, 0.15), rexp(60, 0.08)), # Effect status = rbinom(120, 1, 0.7), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B"), each = 60)) ) ss_results <- MS.ss.analytical( pilot_data = pilot_df_strat_effect, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.80, L = 10, n_start = 100, n_step = 50 ) print(ss_results$results_data)set.seed(456) pilot_df_strat_effect <- data.frame( time = c(rexp(60, 0.15), rexp(60, 0.08)), # Effect status = rbinom(120, 1, 0.7), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B"), each = 60)) ) ss_results <- MS.ss.analytical( pilot_data = pilot_df_strat_effect, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.80, L = 10, n_start = 100, n_step = 50 ) print(ss_results$results_data)
Performs sample size estimation based on a multiplicative model for RMST for stratified trials, using iterative bootstrap simulations.
MS.ss.boot( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )MS.ss.boot( pilot_data, time_var, status_var, arm_var, strata_var, target_power, linear_terms = NULL, L, n_sim = 1000, alpha = 0.05, parallel.cores = 1, patience = 5, n_start = 50, n_step = 25, max_n_per_arm = 2000, verbose = FALSE )
pilot_data |
A |
time_var |
A character string for the time-to-event variable. |
status_var |
A character string for the event status variable. |
arm_var |
A character string for the treatment arm variable. |
strata_var |
A character string for the stratification variable. |
target_power |
A single numeric value for the target power (e.g., 0.80). |
linear_terms |
Optional vector of covariates for the model. |
L |
The numeric truncation time for RMST. |
n_sim |
Number of bootstrap simulations per search step. |
alpha |
The significance level. |
parallel.cores |
Number of cores for parallel processing. |
patience |
Number of consecutive non-improving steps in the search before terminating. |
n_start |
Starting sample size per stratum for the search. |
n_step |
Increment for the sample size search. |
max_n_per_arm |
Maximum sample size per stratum to try. |
verbose |
Logical; if |
This function iteratively searches for the sample size required to
achieve a target_power. At each step of the search, it runs a full bootstrap simulation
(as described in MS.power.boot) to estimate the power for the
current sample size. The search proceeds until the target power is met or
other stopping criteria are satisfied. This process can be very computationally
intensive.
A list containing:
results_data |
A |
results_plot |
A |
results_summary |
A |
status_var should be 1/0. arm_var should be 1/0. strata_var
is a mandatory argument.
pilot_df_strat_effect <- data.frame( time = c(rexp(60, 0.15), rexp(60, 0.08)), # Effect status = rbinom(120, 1, 0.7), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) ss_results <- MS.ss.boot( pilot_data = pilot_df_strat_effect, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.80, L = 10, n_sim = 100, # Low n_sim for example n_start = 100, n_step = 50, patience = 2 ) print(ss_results$results_data)pilot_df_strat_effect <- data.frame( time = c(rexp(60, 0.15), rexp(60, 0.08)), # Effect status = rbinom(120, 1, 0.7), arm = rep(0:1, each = 60), region = factor(rep(c("A", "B", "C"), each = 40)) ) ss_results <- MS.ss.boot( pilot_data = pilot_df_strat_effect, time_var = "time", status_var = "status", arm_var = "arm", strata_var = "region", target_power = 0.80, L = 10, n_sim = 100, # Low n_sim for example n_start = 100, n_step = 50, patience = 2 ) print(ss_results$results_data)
Columns: time, status, arm, age, gender.
A data frame with 250 rows and 5 variables.
data(ph_pwexp_L18_n250) summary(ph_pwexp_L18_n250$time)data(ph_pwexp_L18_n250) summary(ph_pwexp_L18_n250$time)
Columns: time, status, arm, age, gender.
A data frame with 300 rows and 5 variables.
data(ph_weibull_L24_n300) prop.table(table(ph_weibull_L24_n300$status))data(ph_weibull_L24_n300) prop.table(table(ph_weibull_L24_n300$status))
Prints the model metadata and power table returned by rmst.power().
Returns a summary object for printing and further inspection.
Prints and returns the stored ggplot2 object.
## S3 method for class 'rmst_power' print(x, ...) ## S3 method for class 'rmst_power' summary(object, ...) ## S3 method for class 'rmst_power' plot(x, ...)## S3 method for class 'rmst_power' print(x, ...) ## S3 method for class 'rmst_power' summary(object, ...) ## S3 method for class 'rmst_power' plot(x, ...)
x |
An object returned by |
... |
Additional arguments passed to or from other methods. |
object |
An object returned by |
The input object x, invisibly.
An object of class "summary.rmst_power".
A ggplot2 object, invisibly.
Prints the model metadata and sample-size table returned by rmst.ss().
Returns a summary object for printing and further inspection.
Prints and returns the stored ggplot2 object.
## S3 method for class 'rmst_ss' print(x, ...) ## S3 method for class 'rmst_ss' summary(object, ...) ## S3 method for class 'rmst_ss' plot(x, ...)## S3 method for class 'rmst_ss' print(x, ...) ## S3 method for class 'rmst_ss' summary(object, ...) ## S3 method for class 'rmst_ss' plot(x, ...)
x |
An object returned by |
... |
Additional arguments passed to or from other methods. |
object |
An object returned by |
The input object x, invisibly.
An object of class "summary.rmst_ss".
A ggplot2 object, invisibly.
Reads datasets already written in out_dir and reconstructs a
manifest.rds with rich metadata (model, baseline, effects, etc.).
rebuild_manifest( base_recipe, vary, out_dir, filename_template = "sc{scenario_id}_r{rep}" )rebuild_manifest( base_recipe, vary, out_dir, filename_template = "sc{scenario_id}_r{rep}" )
base_recipe |
A validated recipe list (use |
vary |
Named list used originally in |
out_dir |
Directory that already contains the datasets. |
filename_template |
The same template you used when writing files
(default |
The rebuilt manifest (also writes manifest.rds in out_dir).
Expand a recipe over a grid of values (list-only)
recipe_grid(base, vary)recipe_grid(base, vary)
base |
A recipe list. |
vary |
Named list of vectors; keys are dotted paths. |
A list of recipe lists.
r <- recipe_quick_aft(60, "aft_lognormal", baseline=list(mu=2.7, sigma=0.6), treat_effect=-0.2, covariates=list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1)))) recipe_grid(r, list(n=c(60,80), "event_time.effects.treatment"=c(-0.2,-0.4)))r <- recipe_quick_aft(60, "aft_lognormal", baseline=list(mu=2.7, sigma=0.6), treat_effect=-0.2, covariates=list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1)))) recipe_grid(r, list(n=c(60,80), "event_time.effects.treatment"=c(-0.2,-0.4)))
Quick AFT recipe builder for list-based simulation recipes
recipe_quick_aft( n, model = c("aft_lognormal", "aft_weibull"), baseline, treat_effect, covariates, target_censoring = 0.25, allocation = "1:1", seed = NULL )recipe_quick_aft( n, model = c("aft_lognormal", "aft_weibull"), baseline, treat_effect, covariates, target_censoring = 0.25, allocation = "1:1", seed = NULL )
n |
Sample size. |
model |
One of |
baseline |
Baseline parameter list (see model). |
treat_effect |
Numeric treatment coefficient (on log-time scale). |
covariates |
Covariate definitions (list of defs). |
target_censoring |
Target overall censoring fraction (0-1). |
allocation |
Allocation ratio string (e.g., "1:1"). |
seed |
Optional seed retained as recipe/provenance metadata. It does not set the random-number generator. |
A recipe list suitable for simulate_from_recipe.
covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) r <- recipe_quick_aft(120, "aft_lognormal", baseline=list(mu=2.3, sigma=0.5), treat_effect=-0.2, covariates=covs, target_censoring=0.25, allocation="1:1") set.seed(1) dat <- simulate_from_recipe(r)covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) r <- recipe_quick_aft(120, "aft_lognormal", baseline=list(mu=2.3, sigma=0.5), treat_effect=-0.2, covariates=covs, target_censoring=0.25, allocation="1:1") set.seed(1) dat <- simulate_from_recipe(r)
Convenience wrapper that builds a recipe list for proportional-hazard (PH)
models, parallel to recipe_quick_aft for AFT models.
The returned recipe is validated and ready for simulate_from_recipe.
recipe_quick_ph( n, model = c("ph_exponential", "ph_weibull", "ph_gompertz", "cox_pwexp"), baseline, treat_effect, covariates = list(), target_censoring = 0.25, allocation = "1:1", seed = NULL )recipe_quick_ph( n, model = c("ph_exponential", "ph_weibull", "ph_gompertz", "cox_pwexp"), baseline, treat_effect, covariates = list(), target_censoring = 0.25, allocation = "1:1", seed = NULL )
n |
Total sample size (integer). |
model |
PH model. One of |
baseline |
Named list of baseline hazard parameters:
|
treat_effect |
Numeric log-hazard ratio for the treatment arm. |
covariates |
List of covariate definitions. Use |
target_censoring |
Target overall censoring fraction (0-1). Default |
allocation |
Treatment allocation ratio string, e.g. |
seed |
Optional seed retained as recipe/provenance metadata. It does not
set the random-number generator; call |
A recipe list suitable for simulate_from_recipe.
recipe_quick_aft, rmst.sim,
simulate_from_recipe
r <- recipe_quick_ph(100, "ph_weibull", baseline = list(shape = 1.5, scale = 10), treat_effect = -0.5, covariates = list(covar_continuous("age")), target_censoring = 0.30) set.seed(1) dat <- simulate_from_recipe(r)r <- recipe_quick_ph(100, "ph_weibull", baseline = list(shape = 1.5, scale = 10), treat_effect = -0.5, covariates = list(covar_continuous("age")), target_censoring = 0.30) set.seed(1) dat <- simulate_from_recipe(r)
Routes a formula-based call to the matching RMST power routine.
rmst.power( formula, data, arm, sample_sizes, L, strata = NULL, strata_type = c("additive", "multiplicative"), dep_cens = FALSE, type = c("analytical", "boot"), alpha = 0.05, n_sim = 1000L, parallel.cores = 1L, verbose = FALSE, test = NULL )rmst.power( formula, data, arm, sample_sizes, L, strata = NULL, strata_type = c("additive", "multiplicative"), dep_cens = FALSE, type = c("analytical", "boot"), alpha = 0.05, n_sim = 1000L, parallel.cores = 1L, verbose = FALSE, test = NULL )
formula |
A formula of the form |
data |
A |
arm |
Character string naming the treatment arm column (binary 0/1). |
sample_sizes |
Integer vector of per-arm (or per-stratum) sample sizes to evaluate. |
L |
Numeric truncation time for RMST. |
strata |
Character column name, one-sided formula ( |
strata_type |
One of |
dep_cens |
Logical; use dependent-censoring model? Default |
type |
One of |
alpha |
Significance level. Default |
n_sim |
Number of bootstrap replicates (boot methods only). Default |
parallel.cores |
Number of cores for parallel processing. Default |
verbose |
Logical; if |
test |
Which hypothesis test the reported power refers to: |
An object of class c("rmst_power", "list") with elements
results_data, results_plot, results_summary,
model_output, and .meta.
rmst.ss, print.rmst_power,
summary.rmst_power, plot.rmst_power
data(aft_lognormal_L12_n150, package = "RMSTpowerBoost") r <- rmst.power(Surv(time, status) ~ age, data = aft_lognormal_L12_n150, arm = "arm", sample_sizes = c(50, 100, 150), L = 12) print(r) s <- summary(r) plot(r)data(aft_lognormal_L12_n150, package = "RMSTpowerBoost") r <- rmst.power(Surv(time, status) ~ age, data = aft_lognormal_L12_n150, arm = "arm", sample_sizes = c(50, 100, 150), L = 12) print(r) s <- summary(r) plot(r)
A unified, single-call wrapper for generating survival data suitable for use
as reference/pilot data in rmst.power and rmst.ss.
Supports all seven built-in event-time models (AFT and PH families).
rmst.sim( n, model = "aft_lognormal", baseline, treat_effect = 0, covariates = list(), target_censoring = 0.25, allocation = "1:1", L = NULL, seed = NULL )rmst.sim( n, model = "aft_lognormal", baseline, treat_effect = 0, covariates = list(), target_censoring = 0.25, allocation = "1:1", L = NULL, seed = NULL )
n |
Total sample size (split by |
model |
Event-time model. One of:
|
baseline |
Named list of baseline parameters (model-specific; see
|
treat_effect |
Numeric treatment coefficient. Log-time scale for AFT
models; log-hazard ratio for PH models. Default |
covariates |
List of covariate definitions. Elements can be created with
|
target_censoring |
Target overall censoring fraction (0-1). Default |
allocation |
Treatment allocation ratio string, e.g. |
L |
Optional numeric truncation time. Stored as an attribute on the
returned object for downstream use by |
seed |
Optional seed retained as recipe/provenance metadata. It does not
set the random-number generator; call |
Internally routes to recipe_quick_aft for AFT models and
recipe_quick_ph for PH models, then calls
simulate_from_recipe.
A data.frame of class c("rmst_simdata", "data.frame")
with columns time, status, arm (when treatment is
present), and one column per covariate. Attributes:
The validated recipe list used for generation.
The truncation time if supplied, else NULL.
Actual censoring fraction achieved.
rmst.power, rmst.ss,
recipe_quick_ph, recipe_quick_aft,
covar_continuous
set.seed(42) df <- rmst.sim( n = 150, model = "aft_lognormal", baseline = list(mu = 2.2, sigma = 0.5), treat_effect = -0.3, covariates = list(covar_continuous("age"), covar_binary("female")), L = 12, seed = 42 ) print(df) s <- summary(df)set.seed(42) df <- rmst.sim( n = 150, model = "aft_lognormal", baseline = list(mu = 2.2, sigma = 0.5), treat_effect = -0.3, covariates = list(covar_continuous("age"), covar_binary("female")), L = 12, seed = 42 ) print(df) s <- summary(df)
Routes a formula-based call to the matching RMST sample-size routine.
rmst.ss( formula, data, arm, target_power, L, strata = NULL, strata_type = c("additive", "multiplicative"), dep_cens = FALSE, type = c("analytical", "boot"), alpha = 0.05, n_sim = 1000L, parallel.cores = 1L, n_start = 50L, n_step = 25L, max_n = 2000L, patience = 5L, verbose = FALSE, test = NULL )rmst.ss( formula, data, arm, target_power, L, strata = NULL, strata_type = c("additive", "multiplicative"), dep_cens = FALSE, type = c("analytical", "boot"), alpha = 0.05, n_sim = 1000L, parallel.cores = 1L, n_start = 50L, n_step = 25L, max_n = 2000L, patience = 5L, verbose = FALSE, test = NULL )
formula |
A formula of the form |
data |
A |
arm |
Character string naming the treatment arm column (binary 0/1). |
target_power |
Numeric target power (e.g., |
L |
Numeric truncation time for RMST. |
strata |
Character column name, one-sided formula ( |
strata_type |
One of |
dep_cens |
Logical; use dependent-censoring model? Default |
type |
One of |
alpha |
Significance level. Default |
n_sim |
Number of bootstrap replicates. Default |
parallel.cores |
Number of cores for parallel processing. Default |
n_start |
Starting sample size for the search. Default |
n_step |
Search increment. Default |
max_n |
Maximum sample size to try. Default |
patience |
Number of consecutive non-improving steps before stopping. Default |
verbose |
Logical; if |
test |
Which hypothesis test the required sample size refers to: |
An object of class c("rmst_ss", "list") with elements
results_data, results_plot, results_summary,
model_output, and .meta.
rmst.power, print.rmst_ss,
summary.rmst_ss, plot.rmst_ss
data(aft_lognormal_L12_n150, package = "RMSTpowerBoost") r <- rmst.ss(Surv(time, status) ~ age, data = aft_lognormal_L12_n150, arm = "arm", target_power = 0.80, L = 12) print(r)data(aft_lognormal_L12_n150, package = "RMSTpowerBoost") r <- rmst.ss(Surv(time, status) ~ age, data = aft_lognormal_L12_n150, arm = "arm", target_power = 0.80, L = 12) print(r)
Launches the Shiny application bundled with the package. App-specific dependencies are installed only when they are needed.
run_app(install_missing = TRUE, repos = getOption("repos"))run_app(install_missing = TRUE, repos = getOption("repos"))
install_missing |
Logical; if |
repos |
CRAN mirror(s) passed to |
Invisible return value from shiny::runApp().
RMSTpowerBoost::run_app()RMSTpowerBoost::run_app()
Simulate a dataset from a validated recipe (list-only)
simulate_from_recipe(recipe, seed = NULL)simulate_from_recipe(recipe, seed = NULL)
recipe |
A validated recipe list (use |
seed |
Optional integer retained as recipe/provenance metadata; it does
not set the random-number generator. Call |
A data.frame with columns time, status, arm (if treatment present),
plus covariates. Attribute "achieved_censoring" is attached.
covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) rec <- recipe_quick_aft(120, "aft_lognormal", baseline=list(mu=2.2, sigma=0.5), treat_effect=-0.2, covariates=covs, target_censoring=0.25) set.seed(11) dat <- simulate_from_recipe(rec)covs <- list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))) rec <- recipe_quick_aft(120, "aft_lognormal", baseline=list(mu=2.2, sigma=0.5), treat_effect=-0.2, covariates=covs, target_censoring=0.25) set.seed(11) dat <- simulate_from_recipe(rec)
Checks/normalizes the simulation recipe and fills reasonable defaults.
This function does not use or require an analysis horizon; specify L
later in downstream analysis functions.
validate_recipe(recipe)validate_recipe(recipe)
recipe |
A named list defining n, covariates, treatment (optional), event_time (model, baseline, effects, optional frailty), and censoring. |
A validated recipe list.
r <- recipe_quick_aft( n = 100, model = "aft_lognormal", baseline = list(mu = 2.2, sigma = 0.5), treat_effect = -0.2, covariates = list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))), target_censoring = 0.25, allocation = "1:1" ) r2 <- validate_recipe(r)r <- recipe_quick_aft( n = 100, model = "aft_lognormal", baseline = list(mu = 2.2, sigma = 0.5), treat_effect = -0.2, covariates = list(list(name="x", type="continuous", dist="normal", params=list(mean=0, sd=1))), target_censoring = 0.25, allocation = "1:1" ) r2 <- validate_recipe(r)