Overview

In this draft vignette of the weightedsurv package (Leon 2024) we illustrate various weighted survival analyses with two breast cancer datasets available in the base R survival library (Therneau 2023) both evaluating hormon replacement therapy in comparison to chemotherapy. The gbsg trial (Schumacher et al. 1994) was randomized, whereas the rotterdam study (Royston and Altman 2013) was observational; Wherein for comparisons of the non-randomized therapy, in contrast to un-adjusted comparisons, propensity-score adjustment (Cole and Hernán 2004) appears to suggest a stronger impact consistent with the randomized trial.

Survival Analysis and KM Plotting

Survival analysis functions allowing for time-dependent and subject-specific (eg, propensity-scores) weighting.

Weighted estimation for: Cox model; Kaplan-Meier treatment survival curves as well as treatment difference along with point-wise and simultaneous confidence bands; and Restricted mean survival time (RMST, Uno et al. (2014)) comparisons where RMST estimates are evaluated across all potential truncation times (point-wise and simultaneous bands).

Summary: Weighted Survival Analysis Capabilities

Weighted Kaplan-Meier and Log-Rank Analyses
- The code provides functions to compute and plot weighted Kaplan-Meier survival curves for two or more groups, using custom weights. Weighted log-rank statistics can be calculated using flexible, time-dependent weighting schemes (e.g., Fleming-Harrington, Schemper, Xu & O’Quigley, Maggir-Burman, and custom schemes). Functions such as wlr_dhat_estimates, wlr_cumulative, and KM_estimates (in km_wlr_calculations_helpers.R) allow for the calculation of weighted risk sets, event counts, and survival probabilities.
Weighted Cox Proportional Hazards Models
- The code supports fitting Cox models with custom weights, including rho-gamma weights for weighted log-rank tests. Functions like cox_rhogamma and cox_score_rhogamma (in weightedCox.R) implement weighted Cox regression, including score calculations, confidence intervals, and resampling for inference. The code allows for root-finding and score-based estimation in the presence of weights.
Flexible Weighting Schemes
- The weighted_helpers.R file provides a framework for defining, validating, and applying a variety of weighting schemes to survival data. The get_weights and extract_and_calc_weights functions allow users to specify and compute weights for each subject and time point, supporting both standard and custom schemes. Visualization of weight functions over time is supported via plot_weight_schemes.
Counting Process Data Preparation
- The code can transform raw survival data into a counting process format, which is essential for weighted analyses. Functions in df_counting.R and df_counting_helpers.R prepare risk sets, event counts, and weighted summaries for each group, supporting both unweighted and weighted analyses.
Weighted RMST (Restricted Mean Survival Time)
- The code includes functions for calculating cumulative RMST and confidence bands using resampling, with support for weights (rmst_calculations.R).
Plotting and Visualization
- Weighted survival curves, confidence intervals, and risk tables can be plotted using functions in kmplotting_helpers.R. The plotting functions support the display of weighted KM curves, weighted risk tables, and annotation of statistical results.

Analysis functions

df_counting

Prepares a comprehensive dataset for survival analysis using the counting process approach. Computes risk sets, event counts, Kaplan-Meier estimates, log-rank (Fleming-Harrington (\(\rho,\gamma\)) weights) and Cox model results, and quantile estimates for two groups (e.g., treatment and control), optionally handling (subject-specific case) weights (such as inverse probability weights) and stratification.

Time-dependent weights for log-rank statistics are implemented via wt.rg.S function: Supporting a variety of commonly used and custom weighting schemes for weighted log-rank and related tests. The function is flexible and can be used for Fleming-Harrington, Schemper, Xu & O’Quigley (XO), Maggir-Burman (MB), custom time-based, and exponential variants of Fleming-Harrington weights.

Subject-specific (case-weights) are also included in order to implement propensity-score weighted analyses.

Key Features: - Calculates weights for use in weighted log-rank and related survival tests. - Supports standard Fleming-Harrington weights (scheme = "fh"). - Implements Schemper weights (scheme = "schemper"), XO weights (scheme = "XO"), MB weights (scheme = "MB"). - Allows for custom time-based weights (scheme = "custom_time"). - Supports exponential variants of Fleming-Harrington weights (scheme = "fh_exp1", scheme = "fh_exp2").

plot_weighted_km

Creates Kaplan-Meier survival plots for two groups (e.g., treatment vs. control), allowing for the use of weights (such as inverse probability weights) in the estimation. The function can display survival curves, confidence intervals, risk tables, and optionally subgroup curves or additional annotations.

KM_diff

KM_diff compares Kaplan-Meier survival curves between two groups (typically treatment vs. control) using time-to-event data. It calculates survival estimates, their differences, and provides both pointwise and simultaneous confidence intervals. The function can also perform resampling to construct simultaneous confidence bands.

Allows for arbitrary weights (non-negative) to implement (for example) propensity-score adjustment (IPW).

plotKM.band_subgroups

Plots the difference (via KM_diff calculations) in Kaplan-Meier survival curves between two groups (e.g., treatment vs. control), optionally including simultaneous confidence bands and subgroup curves. The function also displays risk tables for the overall population and specified subgroups.

wlr_dhat_estimates

Computes the weighted log-rank test statistic, its variance, the difference in survival between two groups at a specified time point (tzero), the variance of this difference, their covariance, and correlation. It supports flexible time-dependent weighting schemes via the wt.rg.S function.

cumulative_rmst_bands

The cumulative_rmst_bands function calculates and plots cumulative Restricted Mean Survival Time (RMST) estimates and their confidence bands for survival curves, typically comparing two groups (e.g., treatment vs. control).

Key Features

Computes cumulative RMST over time for each group.
Uses resampling (bootstrap or similar) to estimate confidence bands for the RMST curves.
Supports weighted analyses (optional).
Plots the RMST curves and their confidence intervals along with simultaneous bands.

oldpar <- par(no.readonly = TRUE)
library(survival)
# to install weightedKMplots
# library(devtools)
# github_install("larry-leon/weightedsurv")
library(weightedsurv)
library(dplyr)

GBSG data analysis example

—- Data Preparation —- Prepare GBSG data

library(survival)
df_gbsg <- gbsg
df_gbsg$tte <- df_gbsg$rfstime / 30.4375
df_gbsg$event <- df_gbsg$status
df_gbsg$treat <- df_gbsg$hormon
df_gbsg$grade3 <- ifelse(df_gbsg$grade == "3", 1, 0)


# Note that these are default names
# For alternative names, for example if the 
# time-to-event outcome (tte) is time_months
# then tte.name <- "time_months"
tte.name <- "tte"
event.name <- "event"
treat.name <- "treat"
arms <- c("treat", "control")

—- Main Analyses —-

GBSG - ITT Analysis

# Returns standard log-rank (FH(0,0) via scheme="fh" rho,gamma = 0 in scheme_params)

dfcount_gbsg <- df_counting(df=df_gbsg, tte.name = tte.name, event.name = event.name, treat.name = treat.name, arms = arms, 
by.risk = 12, scheme = "fh", scheme_params = list(rho = 0, gamma =0), lr.digits = 4, cox.digits =3)

# test default names
#test <- df_counting(df=df_gbsg, by.risk = 12, scheme = "fh", scheme_params = list(rho = 0, gamma =0), lr.digits = 4, cox.digits =3)

# MB weighting
#test <- df_counting(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, scheme = "MB", scheme_params = list(mb_tstar = 12))
# 0/1 weighting (0 for time <= t.tau, 1 for time > t.tau)
#test <- df_counting(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, scheme = "custom_time", scheme_params = list(t.tau = 6, w0.tau = 0, w1.tau =1))

Show some weight functions

g <- plot_weight_schemes(dfcount_gbsg)
print(g)

—- Plotting —- ## GBSG KM plots

par(mfrow=c(1,2))
# Mine
# Set ymax a little above 1.0 to allow for logrank placement in topleft
# topleft is default; xmed.fraction = 0.65 positions the display of the median estimates below the HR estimates ("topright" legend [default])
# For details, please see summary of xmed.fraction in the Appendix below.

plot_weighted_km(dfcount=dfcount_gbsg, conf.int=FALSE, show.logrank = TRUE, 
                 put.legend.lr = "topleft", ymax = 1.05, xmed.fraction = 0.65)
title(main="GBSG (trial) data un-weighted")


# Compare with survfit
plot_km(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name)
title(main="GBSG (trial) data un-weighted
      via survfit")

Include 95% CIs

plot_weighted_km(dfcount=dfcount_gbsg, conf.int=TRUE, show.logrank = TRUE, ymax = 1.05)
title(main="GBSG (trial) data un-weighted")

Rotterdam observational data analysis with weighting

—- Propensity Score Weighting (Rotterdam) —-

Prepare Rotterdam data

gbsg_validate <- within(rotterdam, {
  rfstime <- ifelse(recur == 1, rtime, dtime)
  t_months <- rfstime / 30.4375
  time_months <- t_months
  status <- pmax(recur, death)
  ignore <- (recur == 0 & death == 1 & rtime < dtime)
  status2 <- ifelse(recur == 1 | ignore, recur, death)
  rfstime2 <- ifelse(recur == 1 | ignore, rtime, dtime)
  time_months2 <- rfstime2 / 30.4375
  grade3 <- ifelse(grade == "3", 1, 0)
  treat <- hormon
  event <- status2
  tte <- time_months
  id <- as.numeric(1:nrow(rotterdam))
  SG0 <- ifelse(er <= 0, 0, 1)
})

Node positive only to correspond to GBSG population

df_rotterdam <- subset(gbsg_validate, nodes > 0)

Baseline demographics table

library(dplyr)
rotterdam2 <- df_rotterdam %>%
  mutate(
    meno = factor(meno, levels = c(0, 1), labels = c("Pre", "Post")),
    grade3 = factor(grade3, levels = c(0, 1), labels = c("Not Grade 3", "Grade 3")),
    chemo = factor(chemo, levels = c(0, 1), labels = c("No", "Yes")),
    er_negative = factor(ifelse(er == 0, 1, 0), labels = c("ER negative","ER positive"))
  )


baseline_table <- create_baseline_table(
  data = rotterdam2,
  treat_var = "treat",
  vars_continuous = c("age", "pgr", "er", "nodes"),
  vars_categorical = c("size","meno","chemo","er_negative"),
  show_pvalue = TRUE,
  show_smd = TRUE
)

baseline_table

Characteristic		Control (n=1207)	Treatment (n=339)	P-value¹	SMD²
Baseline Characteristics by Treatment Arm
age	Mean (SD)	54.1 (13.2)	62.5 (9.9)	<0.001	0.67
pgr	Mean (SD)	169.7 (320.6)	108.2 (200.3)	<0.001	0.21
er	Mean (SD)	160.7 (266.0)	180.6 (271.7)	0.232	0.07
nodes	Mean (SD)	5.1 (5.0)	5.7 (4.6)	0.030	0.13
size				0.581	0.03
	<=20	397 (32.9%)	104 (30.7%)
	20-50	611 (50.6%)	172 (50.7%)
	>50	199 (16.5%)	63 (18.6%)
meno				<0.001	0.31
	Pre	587 (48.6%)	41 (12.1%)
	Post	620 (51.4%)	298 (87.9%)
chemo				<0.001	0.32
	No	655 (54.3%)	311 (91.7%)
	Yes	552 (45.7%)	28 (8.3%)
er_negative				1.000	0.00
	ER negative	1058 (87.7%)	297 (87.6%)
	ER positive	149 (12.3%)	42 (12.4%)
¹ P-values: t-test for continuous, chi-square/Fisher's exact for categorical/binary variables
² SMD = Standardized mean difference

Propensity score model

fit.ps <- glm(treat ~ age + meno + size + grade3 + nodes + pgr + chemo + er, data=df_rotterdam, family="binomial")
pihat <- fit.ps$fitted
pihat.null <- glm(treat ~ 1, family="binomial", data=df_rotterdam)$fitted
wt.1 <- pihat.null / pihat
wt.0 <- (1 - pihat.null) / (1 - pihat)
df_rotterdam$sw.weights <- ifelse(df_rotterdam$treat == 1, wt.1, wt.0)
# truncated weights
df_rotterdam$sw.weights_trunc <- with(df_rotterdam, pmin(pmax(sw.weights, quantile(sw.weights, 0.05)), quantile(sw.weights, 0.95)))

Un-weighted and weighted analyses

dfcount_rotterdam_unwtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24)

dfcount_rotterdam_wtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24, 
                                  weight.name="sw.weights")

—- Plotting Weighted vs Unweighted —-

par(mfrow=c(1,2))
plot_weighted_km(dfcount=dfcount_rotterdam_unwtd, xmed.fraction = 0.65)
title(main="Rotterdam un-weighted K-M curves")
plot_weighted_km(dfcount=dfcount_rotterdam_wtd, xmed.fraction = 0.65)
title(main="Rotterdam weighted K-M curves")

## Compare with GBSG trial data

par(mfrow=c(1,2))
plot_weighted_km(dfcount=dfcount_gbsg, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.05)
title(main="GBSG (trial) data un-weighted K-M curves")
plot_weighted_km(dfcount=dfcount_rotterdam_wtd, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.05)
title(main="Rotterdam weighted K-M curves")

Simultaneous bands and point-wise CIs

Note that the first graph displays 20 (1st 20) resampled approximations to the centered survival difference

\(\Delta(t) = (\widehat{S_1} - \widehat{S_0})(t) - (S_{1} - S_{0})(t)),\) for timepoints \(t\) within “qtau” quartiles of the event times (i.e., for qtau = 2.5% we calculate \(\Delta\) for time points within the 2.5% and 97.5% quantiles).

—- Simultaneous CI bands —-

par(mfrow=c(1,2))
 
 temp <- plotKM.band_subgroups(
    df=df_rotterdam,
    xlabel="Months", ylabel="Survival difference",
    yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7,
    tau_add=42, by.risk = 12, risk_delta=0.075, risk.pad=0.03, ymax.pad = 0.125,
    tte.name = tte.name, treat.name = treat.name, event.name = event.name, weight.name = "sw.weights",
    draws.band = 1000, qtau = 0.025, show_resamples = TRUE
  )
  legend("topleft", c("Difference estimate", "95% pointwise CIs"),
         lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7)
  title(main="Rotterdam data
        propensity score weighted")

Compare Rotterdam observational data with propensity-score weighting and GBSG randomized data