inst/doc/PSsurvival-tutorial.R

In PSsurvival: Propensity Score Methods for Survival Analysis

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

## ----load-package-------------------------------------------------------------
library(PSsurvival)

## ----data---------------------------------------------------------------------
data(simdata_bin)   # Binary treatment (A vs B)
data(simdata_multi) # Four treatment groups (A, B, C, D)

# Binary treatment data
head(simdata_bin)
table(simdata_bin$Z)

# Multiple treatment data
table(simdata_multi$Z)

## ----surveff-basic------------------------------------------------------------
result_bin <- surveff(
  data = simdata_bin,
  ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
  censoring_formula = survival::Surv(time, event) ~ X1 + B1,
  estimand = "overlap",
  censoring_method = "weibull"
)

print(result_bin, max.len = 3)

## ----surveff-multi------------------------------------------------------------
# Define pairwise comparisons
contrast_mat <- matrix(
  c(1, -1, 0, 0,   # A vs B
    1, 0, -1, 0,   # A vs C
    1, 0, 0, -1),  # A vs D
  nrow = 3, byrow = TRUE
)
colnames(contrast_mat) <- c("A", "B", "C", "D")  # Must match treatment levels
rownames(contrast_mat) <- c("A vs B", "A vs C", "A vs D") # row names are optional

result_multi <- surveff(
  data = simdata_multi,
  ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
  censoring_formula = survival::Surv(time, event) ~ X1 + B1,
  estimand = "ATE",
  contrast_matrix = contrast_mat,
  censoring_method = "cox",
  variance_method = "bootstrap",
  B = 50,
  seed = 123
)

print(result_multi, max.len = 3)

## ----estimands, eval=FALSE----------------------------------------------------
# # ATE (default)
# surveff(..., estimand = "ATE")
# 
# # ATT targeting group B
# surveff(..., estimand = "ATT", att_group = "B")
# 
# # Overlap weights
# surveff(..., estimand = "overlap")

## ----trimming, eval=FALSE-----------------------------------------------------
# # Symmetric trimming
# surveff(..., estimand = "ATE", trim = "symmetric", delta = 0.1)
# 
# # Asymmetric trimming
# surveff(..., estimand = "ATE", trim = "asymmetric", alpha = 0.05)

## ----censoring-methods, eval=FALSE--------------------------------------------
# # Weibull censoring model
# surveff(..., censoring_method = "weibull")
# 
# # Cox censoring model (requires bootstrap)
# surveff(..., censoring_method = "cox", variance_method = "bootstrap", B = 200)

## ----variance-methods, eval=FALSE---------------------------------------------
# # Analytical variance (binary + weibull only)
# surveff(..., censoring_method = "weibull", variance_method = "analytical")
# 
# # Bootstrap variance
# surveff(..., variance_method = "bootstrap", B = 200)

## ----boot-level, eval=FALSE---------------------------------------------------
# # Stratified bootstrap
# surveff(..., variance_method = "bootstrap", boot_level = "strata", B = 200)

## ----summary-surveff----------------------------------------------------------
summary(result_bin, conf_level = 0.95, max.len = 5)

## ----summary-returns----------------------------------------------------------
summ <- summary(result_bin, style = "returns")
names(summ)
head(summ$survival_summary$A)
head(summ$difference_summary$`B vs A`)

## ----plot-surv, fig.width=7, fig.height=5-------------------------------------
# Survival curves for all groups
plot(result_multi, type = "surv", curve_width = 0.5)

## ----plot-survdiff, fig.width=7, fig.height=5---------------------------------
# Treatment effect curves
plot(result_multi, type = "survdiff", curve_width = 0.5)

## ----plot-subset-surv, fig.width=7, fig.height=5------------------------------
# Only groups A and C
plot(result_multi, type = "surv", curve_width = 0.5,
     strata_to_plot = c("A", "C"))

## ----plot-subset-diff, fig.width=7, fig.height=5------------------------------
# Only specific contrasts (names must match contrast_matrix rownames exactly)
plot(result_multi, type = "survdiff", curve_width = 0.5,
     strata_to_plot = c("A vs B", "A vs D"))

## ----plot-custom, fig.width=7, fig.height=5-----------------------------------
plot(result_multi,
     type = "surv",
     curve_width = 0.5,
     strata_to_plot = c("A", "B"),
     strata_colors = c("steelblue", "coral"),
     max_time = 15,
     include_CI = TRUE,
     legend_position = "bottom",
     plot_title = "Survival by Treatment Group")

## ----marCoxph-basic-----------------------------------------------------------
hr_result <- marCoxph(
  data = simdata_bin,
  ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
  time_var = "time",
  event_var = "event",
  reference_level = "A",
  estimand = "overlap"
)

print(hr_result)

## ----marCoxph-multi-----------------------------------------------------------
hr_multi <- marCoxph(
  data = simdata_multi,
  ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
  time_var = "time",
  event_var = "event",
  reference_level = "A",
  estimand = "ATE",
  variance_method = "bootstrap",
  B = 50,
  seed = 456
)

print(hr_multi)

## ----marcoxph-estimands, eval=FALSE-------------------------------------------
# # ATE (default)
# marCoxph(..., estimand = "ATE")
# 
# # ATT targeting specific group
# marCoxph(..., estimand = "ATT", att_group = "treated")
# 
# # Overlap weights
# marCoxph(..., estimand = "overlap")

## ----marcoxph-trimming, eval=FALSE--------------------------------------------
# # Symmetric trimming
# marCoxph(..., estimand = "ATE", trim = "symmetric", delta = 0.1)
# 
# # Asymmetric trimming
# marCoxph(..., estimand = "ATE", trim = "asymmetric", alpha = 0.05)

## ----marCoxph-variance, eval=FALSE--------------------------------------------
# # Bootstrap variance (default)
# marCoxph(..., variance_method = "bootstrap", B = 200)
# 
# # Robust sandwich variance (no bootstrap)
# marCoxph(..., variance_method = "robust")

## ----summary-marcoxph---------------------------------------------------------
summary(hr_multi, conf_level = 0.95)

## ----summary-marcoxph-returns-------------------------------------------------
summ_hr <- summary(hr_multi, style = "returns")
names(summ_hr)
summ_hr$HR  # Exponentiated hazard ratios