Estimates the causal treatment effect parameter using g-estimation based on the log-rank test, Cox model, or parametric survival/accelerated failure time (AFT) model. The method uses counterfactual untreated survival times to estimate the causal parameter and derives the adjusted hazard ratio from the Cox model using counterfactual unswitched survival times.
rpsftm(
data,
id = "id",
stratum = "",
time = "time",
event = "event",
treat = "treat",
rx = "rx",
censor_time = "censor_time",
base_cov = "",
psi_test = "logrank",
aft_dist = "weibull",
strata_main_effect_only = TRUE,
low_psi = -2,
hi_psi = 2,
n_eval_z = 101,
treat_modifier = 1,
recensor = TRUE,
admin_recensor_only = TRUE,
autoswitch = TRUE,
gridsearch = TRUE,
root_finding = "brent",
alpha = 0.05,
ties = "efron",
tol = 1e-06,
boot = FALSE,
n_boot = 1000,
seed = NA
)A list with the following components:
psi: The estimated causal parameter.
psi_roots: Vector of psi values at which the Z-statistic
is zero, identified using grid search and linear interpolation.
psi_CI: The confidence interval for psi.
psi_CI_type: The type of confidence interval for psi,
i.e., "grid search", "root finding", or "bootstrap".
logrank_pvalue: The two-sided p-value of the log-rank test
for the ITT analysis.
cox_pvalue: The two-sided p-value for treatment effect based on
the Cox model applied to counterfactual unswitched survival times.
If boot is TRUE, this value represents the
bootstrap p-value.
hr: The estimated hazard ratio from the Cox model.
hr_CI: The confidence interval for hazard ratio.
hr_CI_type: The type of confidence interval for hazard ratio,
either "log-rank p-value" or "bootstrap".
event_summary: A data frame containing the count and percentage
of deaths and switches by treatment arm.
eval_z: A data frame containing the Z-statistics for treatment
effect evaluated at a sequence of psi values. Used to plot and
check if the range of psi values to search for the solution and
limits of confidence interval of psi need be modified.
Sstar: A data frame containing the counterfactual untreated
survival times and event indicators for each treatment group.
The variables include id, stratum,
"t_star", "d_star", "treated", base_cov,
and treat.
kmstar: A data frame containing the Kaplan-Meier estimates
based on the counterfactual untreated survival times by treatment arm.
data_outcome: The input data for the outcome Cox model of
counterfactual unswitched survival times.
The variables include id, stratum, "t_star",
"d_star", "treated", base_cov, and treat.
km_outcome: The Kaplan-Meier estimates of the survival
functions for the treatment and control groups based on the
counterfactual unswitched survival times.
lr_outcome: The log-rank test results for the treatment
effect based on the counterfactual unswitched survival times.
fit_outcome: The fitted outcome Cox model.
fail: Whether a model fails to converge.
psimissing: Whether the psi parameter cannot be estimated.
settings: A list containing the input parameter values.
fail_boots: The indicators for failed bootstrap samples
if boot is TRUE.
fail_boots_data: The data for failed bootstrap samples
if boot is TRUE.
hr_boots: The bootstrap hazard ratio estimates
if boot is TRUE.
psi_boots: The bootstrap psi estimates
if boot is TRUE.
The input data frame that contains the following variables:
id: The subject id.
stratum: The stratum.
time: The survival time for right censored data.
event: The event indicator, 1=event, 0=no event.
treat: The randomized treatment indicator, 1=treatment,
0=control.
rx: The proportion of time on active treatment.
censor_time: The administrative censoring time. It should
be provided for all subjects including those who had events.
base_cov: The baseline covariates (excluding treat).
The name of the id variable in the input data.
The name(s) of the stratum variable(s) in the input data.
The name of the time variable in the input data.
The name of the event variable in the input data.
The name of the treatment variable in the input data.
The name of the rx variable in the input data.
The name of the censor_time variable in the input data.
The names of baseline covariates (excluding
treat) in the input data for the outcome Cox model.
These covariates will also be used in the Cox model for estimating
psi when psi_test = "phreg" and in the parametric
survival regression/AFT model for
estimating psi when psi_test = "lifereg".
The survival function to calculate the Z-statistic, e.g., "logrank" (default), "phreg", or "lifereg".
The assumed distribution for time to event for the AFT
model when psi_test = "lifereg". Options include "exponential",
"weibull" (default), "loglogistic", and "lognormal".
Whether to only include the strata main
effects in the AFT model. Defaults to TRUE, otherwise all
possible strata combinations will be considered in the AFT model.
The lower limit of the causal parameter.
The upper limit of the causal parameter.
The number of points between low_psi and
hi_psi (inclusive) at which to evaluate the Z-statistics.
The optional sensitivity parameter for the constant treatment effect assumption.
Whether to apply recensoring to counterfactual
survival times. Defaults to TRUE.
Whether to apply recensoring to administrative
censoring times only. Defaults to TRUE. If FALSE,
recensoring will be applied to the actual censoring times for dropouts.
Whether to exclude recensoring for treatment arms
with no switching. Defaults to TRUE.
Whether to use grid search to estimate the causal
parameter psi. Defaults to TRUE, otherwise, a root
finding algorithm will be used.
Character string specifying the univariate
root-finding algorithm to use. Options are "brent" (default)
for Brent's method, or "bisection" for the bisection method.
The significance level to calculate confidence intervals.
The method for handling ties in the Cox model, either "breslow" or "efron" (default).
The desired accuracy (convergence tolerance) for psi
for the root finding algorithm.
Whether to use bootstrap to obtain the confidence
interval for hazard ratio. Defaults to FALSE, in which case,
the confidence interval will be constructed to match the log-rank
test p-value.
The number of bootstrap samples.
The seed to reproduce the bootstrap results. The default is
NA, in which case, the seed from the environment will be used.
Kaifeng Lu, kaifenglu@gmail.com
Assuming one-way switching from control to treatment, the hazard ratio and confidence interval under a no-switching scenario are obtained as follows:
Estimate the causal parameter \(\psi\) using g-estimation based on the log-rank test (default), Cox model, or parametric survival/AFT model, using counterfactual untreated survival times for both arms: $$U_{i,\psi} = T_{C_i} + e^{\psi}T_{E_i}$$
Compute counterfactual survival times for control patients using the estimated \(\psi\).
Fit a Cox model to the observed survival times for the treatment group and the counterfactual survival times for the control group to estimate the hazard ratio.
Obtain the confidence interval for the hazard ratio using either
the ITT log-rank test p-value or bootstrap. When bootstrapping,
the interval and p-value are derived from a t-distribution
with n_boot - 1 degrees of freedom.
If grid search is used to estimate \(\psi\), the estimated \(\psi\) is the one with the smallest absolute value among those at which the Z-statistic is zero based on linear interpolation. If root finding is used, the estimated \(\psi\) is the solution to the equation where the Z-statistic is zero.
James M. Robins and Anastasios A. Tsiatis. Correcting for non-compliance in randomized trials using rank preserving structural failure time models. Communications in Statistics. 1991;20(8):2609-2631.
Ian R. White, Adbel G. Babiker, Sarah Walker, and Janet H. Darbyshire. Randomization-based methods for correcting for treatment changes: Examples from the CONCORDE trial. Statistics in Medicine. 1999;18(19):2617-2634.
library(dplyr)
# Example 1: one-way treatment switching (control to active)
data <- immdef %>% mutate(rx = 1-xoyrs/progyrs)
fit1 <- rpsftm(
data, id = "id", time = "progyrs", event = "prog", treat = "imm",
rx = "rx", censor_time = "censyrs", boot = FALSE)
fit1
# Example 2: two-way treatment switching (illustration only)
# the eventual survival time
shilong1 <- shilong %>%
arrange(bras.f, id, tstop) %>%
group_by(bras.f, id) %>%
slice(n()) %>%
select(-c("ps", "ttc", "tran"))
shilong2 <- shilong1 %>%
mutate(rx = ifelse(co, ifelse(bras.f == "MTA", dco/ady,
1 - dco/ady),
ifelse(bras.f == "MTA", 1, 0)))
fit2 <- rpsftm(
shilong2, id = "id", time = "tstop", event = "event",
treat = "bras.f", rx = "rx", censor_time = "dcut",
base_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f"),
low_psi = -3, hi_psi = 3, boot = FALSE)
fit2
Run the code above in your browser using DataLab