Package {rsmart}

rsmart: Sequential Multiple Assignment Randomized Trials Design and Analyses

Description

Implements the interim augmented inverse probability weighted estimator (IAIPWE) for estimating the value of treatment regimes embedded in sequential multiple assignment randomized trials (SMARTs). The IAIPWE subsumes both the inverse probability weighted estimator (IPWE) and the augmented inverse probability weighted estimator (AIPWE), enabling inference at interim and final analyses. The package provides functions for value estimation, sandwich variance computation, group sequential stopping boundaries, and sample size determination for multi-stage SMARTs with up to two treatment options at each stage. See Manschot, Laber, and Davidian (2023) doi:10.1111/biom.13854 for additional details.

Author(s)

Maintainer: Cole Manschot cole.manschot@msd.com

Other contributors:

• Tabitha Peter [contributor]

• Gabriela Piasecki [contributor]

• Laura Frederick [contributor]

See Also

Useful links:

• https://MSDLLCpapers.github.io/rsmart/

• https://github.com/MSDLLCpapers/rsmart

• https://msdllcpapers.github.io/rsmart/

• Report bugs at https://github.com/MSDLLCpapers/rsmart/issues

Create a blocked randomization function

Description

Returns a randomization function compatible with sim_treatment that implements permuted block randomization. Within each block, every treatment appears exactly block_rep times, and the order within blocks is randomly permuted. This ensures approximate balance across treatments at any point during enrollment.

Usage

block_rand(block_rep = 2)

Arguments

Value

A function with signature function(n, n_treatments, prob) suitable for use as rand_prob_fn in sim_treatment. The prob argument is accepted but ignored, since blocked randomization enforces equal allocation within each block.

Examples

# Use blocked randomization with sim_treatment
set.seed(1)
a <- sim_treatment(n = 100, n_treatments = 2,
                   rand_prob_fn = block_rand(block_rep = 2))
table(a)

# Three treatments with block size 6 (block_rep = 2)
set.seed(1)
a <- sim_treatment(n = 300, n_treatments = 3,
                   rand_prob_fn = block_rand(block_rep = 2))
table(a)

# Stage 2 blocked randomization within prior treatment groups
set.seed(1)
df <- data.frame(a1 = c(rep(0, 50), rep(1, 50)))
df <- sim_treatment(n_treatments = 2, dat = df, stage = 2,
                    rand_prob_fn = block_rand(block_rep = 3))
table(df$a1, df$a2)

Estimate values for all treatment regimes

Description

Wraps the outcome regression and value term functions to estimate the value of all regimes in the provided list. For each regime, Q-functions are fitted (if specified), augmentation and IPW terms are computed, and the regime value is estimated.

Usage

estimate_values(df, q_list, regime_all, feasible_sets_indicator, pis, nus)

Arguments

Value

A list with the following components:

value

A numeric vector of estimated regime values.

df

A list of matrices of value term components, augmentation terms 1 through 2K then the IPW term for each regime.

q_all

A list of fitted Q-function objects for each regime, or an empty list if q_list is NULL.

Generate sample data for a two-stage SMART where responders are not re-randomized

Description

This function generates data for a two-stage Sequential Multiple Assignment Randomized Trial (SMART) where responders are not re-randomized in the second stage. The function allows for different value patterns and treatment assignments.

Usage

gen_no_trt_resp(n, s2 = 100, block_rep = 2, r2p = 0.5)

Arguments

Value

A data frame with n rows and the following columns:

t1

Study day of enrollment.

t2

Study day of second-stage randomization.

t3

Study day of outcome assessment.

a1

First-stage treatment assignment (0 or 1).

r2

Response indicator at stage 2 (0 = non-responder, 1 = responder).

a2

Second-stage treatment assignment (0 or 1 for non-responders, 0 for responders).

x11

Baseline covariate (continuous, range 25–75).

x12

Baseline covariate (binary, 0 or 1).

x21

Stage 2 covariate (continuous, range 0–1).

y

Continuous outcome (higher is better).

id

Individual identifier (1 to n).

Examples

set.seed(1)
dat <- gen_no_trt_resp(n=400, s2=100, block_rep=2, r2p = 0.5)
head(dat)

Compute the An matrix for the sandwich variance estimator

Description

Computes the An matrix, defined as -\partial \Psi / \partial \theta, where \Psi is the stacked estimating equations vector. The estimating equations are ordered to match Bn: \pi_1, \ldots, \pi_K, \nu_1, \ldots, \nu_{K+1}, \beta^{(\ell)}_1, \ldots, \beta^{(\ell)}_K for each regime \ell, and the value estimating equations V_1, \ldots, V_L. Based on Section 7 of Boos and Stefanski for the robust sandwich matrix.

Usage

get_an(
  df,
  pis,
  p_fits,
  nus,
  q_all,
  values,
  regime_all,
  dfs,
  feasible_sets_indicator,
  q_list
)

Arguments

Value

A numeric matrix representing the An component of the sandwich variance estimator, with dimensions equal to the total number of estimated parameters.

Compute the Bn matrix for the sandwich variance estimator

Description

Computes the Bn matrix, defined as n^{-1} \sum_{i=1}^{n} \Psi_i \Psi_i^T, the empirical variance of the estimating equations. Each row of the individual-level estimating equation matrix corresponds to one subject, with columns ordered as \pi_1, \ldots, \pi_K, \nu_1, \ldots, \nu_{K+1}, \beta^{(\ell)}_1, \ldots, \beta^{(\ell)}_K for each regime \ell, and the value estimating equations V_1, \ldots, V_L. Based on Section 7 of Boos and Stefanski for the robust sandwich matrix.

Usage

get_bn(df, p_fits, nus, q_all, dfs, values, feasible_sets_indicator, q_list)

Arguments

Value

A numeric matrix representing the Bn component of the sandwich variance estimator, with dimensions equal to the total number of estimated parameters.

Compute all stopping boundaries for a two-analysis group sequential design

Description

Wrapper function that computes the stopping boundaries for both the first and second analyses of a group sequential design with multiple treatment regimes. Calls get_first_bound and get_next_bound sequentially.

Usage

get_bounds(
  alpha = 0.05,
  inf_frac = c(0.5, 1),
  spend_fn = "OF",
  corr = diag(x = 1, nrow = 1, ncol = 1),
  test_type = "one-sided",
  lambda = 0.1,
  tol = 1e-06,
  max_iter = 1000
)

Arguments

Value

A list with the following components:

bounds / bound

A numeric value (single analysis) or numeric vector (sequential) of stopping boundaries. The key is "bounds" for a single analysis and "bound" for multiple analyses.

spending

A numeric value or vector of cumulative alpha spent at each analysis.

convergence

A logical value or vector indicating whether the algorithm converged at each analysis.

iterations

An integer value or vector of the number of iterations used at each analysis.

alpha

The input type I error rate.

inf_frac

The input information fractions.

spend_fn

The input spending function name.

corr

The input correlation matrix.

test_type

The input test type.

Compute the first chi-squared stopping boundary for a group sequential design

Description

Determines the first analysis stopping boundary based on a chi-squared global test statistic. This is used when the comparison type is a global chi-squared test (testing whether any regime differs from the others) rather than individual regime comparisons against a fixed control.

Usage

get_bounds_chi(
  alpha = 0.05,
  inf_frac = c(0.5, 1),
  spend_fn = "OF",
  corr = diag(x = 1, nrow = 1, ncol = 1),
  mu = NULL,
  B = 1000001,
  seed = 1
)

Arguments

Details

The function uses Monte Carlo simulation to estimate the boundary that controls the familywise type I error rate at a specified level. A contrast matrix C is constructed to form L-1 linearly independent comparisons among L regimes, and the chi-squared statistic is computed as (CZ)' (C \Sigma C')^{-1} (CZ).

Value

A list with the following components:

bound

A numeric vector of chi-squared stopping boundaries, one per analysis.

spending

A numeric vector of observed cumulative alpha spent at each analysis.

typeI

The overall achieved type I error rate across all analyses.

convergence

A logical value; always TRUE for Monte Carlo based computation.

iters

The number of Monte Carlo samples used (B).

dfchi

The degrees of freedom of the chi-squared statistic, equal to the rank of C \Sigma C'.

alpha

The input type I error rate.

inf_frac

The input information fractions.

spend_fn

The input spending function name.

corr

The input correlation matrix.

mu

The input mean vector.

B

The input number of Monte Carlo samples.

seed

The input random seed.

Compute the first stopping boundary for a group sequential design

Description

Determines the first analysis stopping boundary that controls the familywise type I error rate at a specified level, adjusting for the multiplicity of multiple treatment regimes. Uses an iterative search to find the boundary value. This function works only for superiority (one-sided) testing.

Usage

get_first_bound(
  alpha = 0.05,
  inf_frac = c(0.5, 1),
  spend_fn = "OF",
  corr = diag(x = 1, nrow = 1, ncol = 1),
  test_type = "one-sided",
  lambda = 0.1,
  tol = 1e-06,
  max_iter = 1000
)

Arguments

Value

A list with the following components:

bound

The computed stopping boundary for the first analysis.

typeI

The achieved type I error rate at the boundary.

convergence

A logical value indicating whether the algorithm converged.

iters

The number of iterations used.

Compute the kappa (stage reached) for each individual

Description

Determines the total number of stages each individual has reached by a specified analysis time t_s, based on their arrival times t1, ..., t_{K+1}.

Usage

get_kappa(df, t_s, K)

Arguments

Value

An integer vector of length nrow(df) indicating the number of stages each individual has reached by time t_s. Values range from 0 (not yet enrolled) to K+1 (outcome observed).

Compute subsequent stopping boundaries for a group sequential design

Description

Given the boundary from the first analysis, determines the stopping boundary for the s-th analysis that controls the overall type I error rate. Uses an iterative search with the joint distribution of test statistics across analyses.

Usage

get_next_bound(
  alpha = 0.05,
  inf_frac = c(0.5, 1),
  spend_fn = "OF",
  corr = diag(x = 1, nrow = 1, ncol = 1),
  test_type = "one-sided",
  lambda = 0.1,
  tol = 1e-06,
  prev_bound = NULL,
  s = 2,
  max_iter = 1000
)

Arguments

Value

A list with the following components:

bound

The computed stopping boundary for analysis s.

typeI

The achieved cumulative type I error rate through analysis s.

convergence

A logical value indicating whether the algorithm converged.

iters

The number of iterations used.

Estimate stage arrival probabilities (nu)

Description

Estimates the probability that an individual has reached each stage given they are enrolled in the trial. Requires the data frame to have a kappa column indicating the stage reached. The returned list uses indexing such that nu[[k]] corresponds to stage k, and nu[[K+1]] is the probability that an individual has their final outcome observed.

Usage

get_nu(df, K)

Arguments

Value

A list with the following components:

nu

A list of length K+1 where nu[[k]] is the estimated probability that an individual has reached stage k given enrollment.

ns

The total number of individuals enrolled in the trial (kappa > 0).

nd

The proportion of individuals who have reached their last treatment stage and have their outcome observed.

Extract coefficients from all fitted Q-function models

Description

Extracts and concatenates the regression coefficients from all fitted Q-function models across all regimes and stages.

Usage

get_q_coefs(q_all)

Arguments

Value

A numeric vector of all Q-function regression coefficients, ordered by regime and then by stage.

Fit outcome regression (Q-function) models across all stages for a single regime

Description

Fits Q-function models backwards from the last stage to the first for a single treatment regime. At each stage, both regime-modified and unmodified predicted values are computed. Handles feasible sets (where responders may not be re-randomized) and interim analyses.

Usage

get_q_fits(df, q_list, regime, feasible_sets_indicator = FALSE)

Arguments

Value

A list with the following components:

q_fits

A list of fitted modelObj objects, one per stage. When feasible sets with multiple models are used, the element is a named list with "r0" and "r1".

mod_regime_vhats

A matrix of regime-modified predicted values with columns q1, ..., q_{K+1}.

unmod_regime_vhats

A matrix of unmodified predicted values with columns q1_nochange, ..., q_{K+1}_nochange.

Determine sample size for a group sequential SMART design

Description

For specified operating characteristics, iteratively increases the sample size until the desired power is achieved. Assumes the information fraction remains unchanged as the sample size increases, which is reasonable for small changes or at the design stage.

Usage

get_sample_size(
  variances,
  beta,
  delta,
  bounds,
  n_init = 100,
  corr = bounds$corr,
  inf_frac = bounds$inf_frac,
  n_split = NULL
)

Arguments

Value

A list with the following components:

N

A numeric vector of sample sizes at each analysis.

power

The achieved power at the final sample size.

prop_rej

A numeric vector of cumulative rejection probabilities at each analysis.

variances

The input variances.

beta

The input type II error rate.

delta

The input alternative differences.

bounds

The input stopping boundaries.

n_init

The input initial sample size.

corr

The input correlation matrix.

inf_frac

The input information fractions.

n_split

The input sample size split proportions.

Determine sample size for a chi-squared global test in a group sequential SMART design

Description

For specified operating characteristics, iteratively increases the sample size until the desired power is achieved using a chi-squared global test statistic. The chi-squared test assesses whether any treatment regime differs from the others, rather than comparing individual regimes against a fixed control.

Usage

get_sample_size_chi(
  variances,
  beta,
  delta,
  bounds,
  n_init = 100,
  corr = bounds$corr,
  inf_frac = bounds$inf_frac,
  n_split = NULL,
  seed = 2,
  B = 10001,
  lambda = 20,
  max_iter = 100
)

Arguments

Details

The function uses Monte Carlo simulation to estimate power at each candidate sample size. A contrast matrix C is constructed to form L-1 linearly independent comparisons among L regimes, and the chi-squared statistic is computed as (CZ)' (C \Sigma C')^{-1} (CZ).

Value

A list with the following components:

N

A numeric vector of sample sizes at each analysis.

power

The achieved power at the final sample size.

prop_rej

A numeric vector of cumulative rejection probabilities at each analysis.

variances

The input variances.

beta

The input type II error rate.

delta

The input alternative differences.

bounds

The input stopping boundaries.

n_init

The input initial sample size.

corr

The input correlation matrix.

inf_frac

The input information fractions.

n_split

The input sample size split proportions.

seed

The input random seed.

B

The input number of Monte Carlo samples.

lambda

The input step size for the sample size search.

IAIPWE for K-stage SMARTs with up to 2 treatment options at each stage

Description

The IAIPWE was developed for arbitrary K stage SMART and subsumes both IPWE and AIPWE. To implement the IPWE, q_list should be NULL and the t_s should be set to the maximum available time of the outcome observed. To implement the AIPWE, t_s should be set to the maximum available time of the outcome observed.

Usage

iaipwe(df, pi_list, q_list, regime_all, feasible_sets_indicator, t_s, B = NULL)

Arguments

Details

If times are not recorded for the study, "dummy" times can be used with the analysis time set to be the maximum of the dummy times.

Value

A list with the following components:

values

A numeric vector of estimated regime values.

se

A numeric vector of standard errors for the regime values.

covariance

The estimated L x L covariance matrix of the regime value estimators.

params

A numeric vector of all estimated parameters.

nus

A list of estimated stage arrival probabilities, as returned by get_nu().

q_all

A list of fitted Q-function objects for each regime.

regime_all

The input regime list, returned for convenience.

dfs

A list of value term matrices for each regime.

variance_choice

A character string indicating the variance estimation method used ("Asymptotic" or "Bootstrap").

chi_square

A list with Statistic (the chi-squared test statistic for equality of regime values) and p_value.

pcsttrial: Simulated Clinical Trial Data for Pain Coping Skills Training

Description

A dataset containing simulated patient data to mimic the characteristics of the Pain Coping Skills Training trial presented in Manschot, Laber, and Davidian (2023).

Usage

pcsttrial

Format

A data frame with 284 rows and 14 variables:

pctchange

Percent change in outcome from baseline at final assessment

height

Patient height (cm)

weight

Patient weight (kg)

comorbidity

Comorbidity indicator (0 = no, 1 = yes)

painmed

Pain medication use indicator (0 = no, 1 = yes)

chemo

Chemotherapy indicator (0 = no, 1 = yes)

pctchangek2

Percent change in outcome at stage 2

adherence

Treatment adherence measure, from 0.5 to 1 in increments of 0.1

a1

First-stage treatment assignment

a2

Second-stage treatment assignment

r2

Response indicator at stage 2 (0 = non-responder, 1 = responder)

study_day_enroll

Study day of enrollment

study_day_rerand

Study day of re-randomization

study_day_outcome

Study day of outcome assessment

Details

Note: this dataset is for illustrative purposes and does not use the actual trial data.

Fit propensity score models for all stages

Description

Fits propensity score models at each stage of the SMART, using only individuals who have reached that stage. Returns estimated propensity scores and fitted model objects for all stages.

Usage

pi_fits(df, p_list)

Arguments

Value

A list with the following components:

ps

A data frame of estimated propensity scores with columns pi1, ..., piK. Individuals who have not reached a stage are assigned a value of 99.

p_fits

A list of fitted modelObj objects, one per stage.

Fit a propensity score model for a single stage

Description

Fits a propensity score model at a single stage of the SMART and returns the fitted probabilities. The probability that each individual received the treatment they actually received is computed.

Usage

pstep(pmodel, data, response, k)

Arguments

Value

A list with the following components:

pk

A numeric vector of predicted probabilities that treatment is 1 at stage k.

ps

A numeric vector of estimated propensity scores (probability that the individual received the treatment they actually received).

pfit

The fitted modelObj object.

Fit an outcome regression (Q-function) model for a single stage

Description

Fits a single Q-function model at one stage of the SMART. Predictions are obtained both for the unmodified data and for data where the treatment is set to the recommended regime.

Usage

qstep(qmodel, data, response, newdata, regime, txName)

Arguments

Value

A list with the following components:

hats_mod

A numeric vector of predicted values under the regime-consistent treatment.

hats_unmod

A numeric vector of predicted values under the actual (unmodified) treatment.

qfit

The fitted modelObj object.

Generate regime lists for a two-stage SMART where responders are not re-randomized

Description

Constructs the treatment assignment matrices and consistency indicator matrices for each embedded regime in a two-stage SMART. For stages where response status determines treatment deterministically, the regime matrix is updated to reflect the responder treatment.

Usage

regime_list_no_trt_resp(
  emb_regimes,
  dat,
  resp_trt = list(r2 = list(0, 0, 0, 0))
)

Arguments

Value

A list of length length(emb_regimes). Each element is a list with two components:

regime

An n x K matrix of treatment assignments each individual would receive if they followed that regime.

regime_ind

An n x K indicator matrix (0 or 1) of whether each individual's observed treatment was consistent with the regime at each stage.

Simulate treatment assignments for a single stage

Description

Generates treatment assignments for n individuals given a specified number of possible treatments and a randomization probability function. This is a utility function for constructing SMART simulations with flexible randomization schemes. If an input data frame is provided, the treatment assignments are appended as a new column a<stage>.

Usage

sim_treatment(
  n = NULL,
  n_treatments,
  rand_prob_fn = function(n, n_treatments, prob) {
     sample(0:(n_treatments - 1), size
    = n, replace = TRUE, prob = prob)
 },
  dat = NULL,
  stage = 1,
  randomize_response = NULL,
  prob = NULL
)

Arguments

Value

If dat is NULL (default), an integer vector of length n with treatment assignments coded as 0, 1, ..., n_treatments - 1. If dat is provided, the input data frame with a new integer column a<stage> appended.

Examples

# Default: Bernoulli randomization with 2 treatments (prob 0.5)
set.seed(1)
a <- sim_treatment(n = 100, n_treatments = 2)
table(a)

# Unequal randomization: 70% to treatment 0, 30% to treatment 1
set.seed(1)
a <- sim_treatment(n = 100, n_treatments = 2, prob = c(0.7, 0.3))
table(a)

# Three treatments with equal probability
set.seed(1)
a <- sim_treatment(n = 300, n_treatments = 3)
table(a)

# With an input data frame
set.seed(1)
df <- data.frame(x1 = rnorm(100), x2 = rbinom(100, 1, 0.5))
df <- sim_treatment(n_treatments = 2, dat = df)
head(df)

# With an input data frame at stage 2
set.seed(1)
df <- data.frame(x1 = rnorm(100), a1 = rbinom(100, 1, 0.5))
df <- sim_treatment(n_treatments = 2, dat = df, stage = 2)
head(df)


# Stage 2 with different randomization probabilities depending on a1:
# equal (50/50) if a1 == 0, unequal (70/30) if a1 == 1
# Using a named list for prob, where names match group levels
set.seed(1)
n <- 200
df <- data.frame(a1 = rbinom(n, 1, 0.5))
df <- sim_treatment(n_treatments = 2, dat = df, stage = 2,
                    prob = list("0" = c(0.5, 0.5), "1" = c(0.7, 0.3)))
table(df$a1, df$a2)

Compute stopping boundaries and sample size for a group sequential SMART

Description

Calculates the stopping boundaries and required sample size for a group sequential design with multiple treatment regimes embedded in a SMART. The boundaries are found first to control the type I error rate, then the sample size is determined to achieve the desired power. This function wraps get_bounds and get_sample_size.

Usage

smart_design(
  test_type = "one-sided",
  comp_type = "fixed.control",
  alpha = 0.05,
  beta = 0.1,
  delta = NULL,
  n_init = 400,
  inf_frac = 1,
  spend_fn = "OF",
  corr = NULL,
  variances = NULL,
  lambdaB = 0.1,
  lambdaC = 20,
  tol = 1e-06,
  max_iterB = 1000,
  max_iterC = 100,
  seed = 1,
  Bb = 100001,
  Bc = 10001,
  mu = NULL
)

Arguments

Value

A list with the following components:

boundaries

The output from get_bounds (when comp_type = "fixed.control") or get_bounds_chi (when comp_type = "global.chi.sq"), containing the stopping boundaries and associated metadata. See those functions for details.

sample.size

The output from get_sample_size or get_sample_size_chi, containing the required sample sizes, achieved power, cumulative rejection probabilities, and echoed input parameters. If delta is NULL (no alternative specified), the sample size is not computed and this element is NULL.

inputs

A list echoing all input parameters passed to smart_design.

Compute the individual-level value terms for a single regime

Description

Computes the 2K+1 coarsening-level value terms for each individual under a single treatment regime. These include augmentation terms for levels r = 1, \ldots, 2K and the IPW term for R = \infty (i.e., R = 2K+1). When outcome regression models (qs) are provided, the augmented terms are computed; otherwise only the IPW term is non-zero.

Usage

value_terms(df, regime_ind, pis, qs, nus)

Arguments

Value

A numeric matrix with nrow(df) rows and 2K+1 columns, where each column corresponds to a coarsening-level term in the value estimator.