PSsurvival: Propensity Score Methods for Survival Analysis

Introduction

PSsurvival implements propensity score methods for causal inference with time-to-event outcomes. The package provides two main functions:

Both functions support binary and multiple treatment groups, multiple weighting schemes, and propensity score trimming.

library(PSsurvival)

Data

The package includes two simulated datasets:

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

# Binary treatment data
head(simdata_bin)
#>        X1      X2      X3 B1 B2 Z    time event
#> 1  0.9820 -0.0517  1.0119  0  1 A 10.0035     1
#> 2  0.4687 -0.4994 -1.2844  0  0 A  0.4150     0
#> 3 -0.1080 -0.9025 -0.5866  1  0 B  7.6058     1
#> 4 -0.2129 -0.3720 -2.3708  0  0 B 10.0073     0
#> 5  1.1581 -0.1884 -0.3767  1  1 B 17.3095     1
#> 6  1.2924 -0.7366 -2.0261  0  1 A 14.5914     1
table(simdata_bin$Z)
#> 
#>   A   B 
#> 408 592

# Multiple treatment data
table(simdata_multi$Z)
#> 
#>   A   B   C   D 
#> 244 188 208 360

Both datasets contain:

Counterfactual Survival Analysis with surveff()

Key Arguments

Data and model specification:

Weighting and trimming:

Variance estimation:

Basic Usage: Binary Treatment

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)
#> 
#> === Survival Effect Estimation ===
#> 
#> Treatment groups: A, B 
#> Number of groups: 2 
#> Estimand: overlap 
#> Censoring method: weibull 
#> Variance method: analytical 
#> Evaluation times: 666 time points
#> 
#> Survival function estimates:
#>           A      B
#> [1,] 0.9979 1.0000
#> [2,] 0.9960 1.0000
#> [3,] 0.9960 0.9973
#>   ... (663 more rows not shown)
#> 
#> Treatment effect estimates:
#>      B vs A
#> [1,] 0.0021
#> [2,] 0.0040
#> [3,] 0.0013
#>   ... (663 more rows not shown)

The ps_formula specifies the propensity score model (treatment ~ covariates). The censoring_formula models the censoring distribution within each treatment group.

Multiple Treatment Groups

For multiple groups, survival differences require a contrast_matrix specifying which comparisons to make. Each row represents one contrast with exactly two non-zero elements: 1 and -1. Column names must match treatment levels exactly.

# 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)
#> 
#> === Survival Effect Estimation ===
#> 
#> Treatment groups: A, B, C, D 
#> Number of groups: 4 
#> Estimand: ATE 
#> Censoring method: cox 
#> Variance method: bootstrap 
#> Evaluation times: 667 time points
#> 
#> Survival function estimates:
#>           A      B C D
#> [1,] 1.0000 0.9945 1 1
#> [2,] 0.9961 0.9945 1 1
#> [3,] 0.9917 0.9945 1 1
#>   ... (664 more rows not shown)
#> 
#> Treatment effect estimates:
#>       A vs B  A vs C  A vs D
#> [1,]  0.0055  0.0000  0.0000
#> [2,]  0.0016 -0.0039 -0.0039
#> [3,] -0.0028 -0.0083 -0.0083
#>   ... (664 more rows not shown)

Without contrast_matrix, only group-specific survival functions are estimated (no differences in survivals).

Estimand Options

Three estimands are supported and estimated by corresponding weighting methods:

# ATE (default)
surveff(..., estimand = "ATE")

# ATT targeting group B
surveff(..., estimand = "ATT", att_group = "B")

# Overlap weights
surveff(..., estimand = "overlap")

Propensity Score Trimming

Trimming excludes observations with extreme propensity scores. Not available with overlap weights (which already downweight extremes).

# Symmetric trimming
surveff(..., estimand = "ATE", trim = "symmetric", delta = 0.1)

# Asymmetric trimming
surveff(..., estimand = "ATE", trim = "asymmetric", alpha = 0.05)

Censoring Methods

Two methods are available for estimating censoring distributions:

# Weibull censoring model
surveff(..., censoring_method = "weibull")

# Cox censoring model (requires bootstrap)
surveff(..., censoring_method = "cox", variance_method = "bootstrap", B = 200)

No censoring adjustment: To skip censoring adjustment entirely (e.g., when there is no censoring), use censoring_formula = Surv(time, event) ~ 0. In this case, all censoring scores are set to 1.

Variance Estimation

Two variance estimation methods are available:

# Analytical variance (binary + weibull only)
surveff(..., censoring_method = "weibull", variance_method = "analytical")

# Bootstrap variance
surveff(..., variance_method = "bootstrap", B = 200)

Bootstrap sampling level (boot_level): Controls how bootstrap samples are drawn.

# Stratified bootstrap
surveff(..., variance_method = "bootstrap", boot_level = "strata", B = 200)

Summary Output

The summary() method provides two output styles:

Printed output (style = "prints", default): Displays formatted tables with confidence intervals. Customize with:

Programmatic output (style = "returns"): Returns a list of matrices for further processing:

summary(result_bin, conf_level = 0.95, max.len = 5)
#> 
#> === Survival Effect Estimation Summary ===
#> 
#> Treatment groups: A, B 
#> Number of groups: 2 
#> Sample size: 1000 complete cases
#> Estimand: overlap 
#> Censoring method: weibull 
#> Variance method: analytical 
#> Evaluation times: 666 time points
#> Confidence level: 0.95 
#> 
#> --- Survival Function Estimates ---
#> 
#> Group A :
#>    Time Estimate     SE 95%CI.lower 95%CI.upper
#>  0.0492   0.9979 0.0021      0.9937           1
#>  0.0655   0.9960 0.0028      0.9905           1
#>  0.1875   0.9960 0.0028      0.9905           1
#>  0.3118   0.9960 0.0028      0.9905           1
#>  0.3199   0.9960 0.0028      0.9905           1
#>   ... (661 more rows not shown)
#> 
#> Group B :
#>    Time Estimate     SE 95%CI.lower 95%CI.upper
#>  0.0492   1.0000 0.0000      1.0000           1
#>  0.0655   1.0000 0.0000      1.0000           1
#>  0.1875   0.9973 0.0027      0.9920           1
#>  0.3118   0.9953 0.0033      0.9888           1
#>  0.3199   0.9938 0.0037      0.9866           1
#>   ... (661 more rows not shown)
#> 
#> --- Treatment Effect Estimates ---
#> 
#> B vs A :
#>    Time Estimate     SE 95%CI.lower 95%CI.upper
#>  0.0492   0.0021 0.0021     -0.0020      0.0063
#>  0.0655   0.0040 0.0028     -0.0015      0.0095
#>  0.1875   0.0013 0.0039     -0.0064      0.0089
#>  0.3118  -0.0007 0.0044     -0.0092      0.0079
#>  0.3199  -0.0022 0.0046     -0.0113      0.0068
#>   ... (661 more rows not shown)

For programmatic access:

summ <- summary(result_bin, style = "returns")
names(summ)
#> [1] "survival_summary"   "difference_summary"
head(summ$survival_summary$A)
#>     Time  Estimate          SE 95%CI.lower 95%CI.upper
#> 1 0.0492 0.9978679 0.002126095   0.9937008           1
#> 2 0.0655 0.9960009 0.002822286   0.9904693           1
#> 3 0.1875 0.9960009 0.002822286   0.9904693           1
#> 4 0.3118 0.9960009 0.002822286   0.9904693           1
#> 5 0.3199 0.9960009 0.002822286   0.9904693           1
#> 6 0.3958 0.9934268 0.003811892   0.9859556           1
head(summ$difference_summary$`B vs A`)
#>     Time      Estimate          SE  95%CI.lower 95%CI.upper
#> 1 0.0492  0.0021321189 0.002126095 -0.002034951 0.006299189
#> 2 0.0655  0.0039991123 0.002822286 -0.001532467 0.009530692
#> 3 0.1875  0.0012872620 0.003897586 -0.006351866 0.008926390
#> 4 0.3118 -0.0006676973 0.004358730 -0.009210651 0.007875257
#> 5 0.3199 -0.0022366401 0.004632447 -0.011316070 0.006842790
#> 6 0.3958  0.0003374480 0.005295230 -0.010041012 0.010715908

Plotting Results

The plot() method produces survival curves (type = "surv") or treatment effect curves (type = "survdiff").

# Survival curves for all groups
plot(result_multi, type = "surv", curve_width = 0.5)

# Treatment effect curves
plot(result_multi, type = "survdiff", curve_width = 0.5)

Selecting specific strata: Use strata_to_plot to display a subset of groups or contrasts.

# Only groups A and C
plot(result_multi, type = "surv", curve_width = 0.5,
     strata_to_plot = c("A", "C"))

# 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"))

Customization options:

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")

Marginal Hazard Ratios with marCoxph()

marCoxph() fits a weighted Cox model to estimate marginal hazard ratios.

Key Arguments

Data and model specification:

Weighting and trimming:

Variance estimation:

Basic Usage: Binary Treatment

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)
#> 
#> === Marginal Cox Proportional Hazards Model ===
#> 
#> Treatment variable: Z 
#> Treatment levels: A, B 
#> Reference level: A 
#> Number of groups: 2 
#> Sample size: 1000 complete cases
#> Estimand: overlap 
#> Variance method: bootstrap-full 
#> Bootstrap iterations: 100 
#> Trimming: None
#> 
#> Sample sizes used in Cox model:
#>  treatment   n events
#>          A 408    250
#>          B 592    416
#> 
#> Log hazard ratio estimates:
#>  treatment  logHR SE_robust SE_bootstrap
#>        Z:B 0.5262    0.0842       0.0844
#> 
#> Note: Use summary() for confidence intervals and exponentiated hazard ratios.

The reference_level argument is mandatory and specifies the baseline group for hazard ratio comparisons.

Multiple Treatment Groups

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)
#> 
#> === Marginal Cox Proportional Hazards Model ===
#> 
#> Treatment variable: Z 
#> Treatment levels: A, B, C, D 
#> Reference level: A 
#> Number of groups: 4 
#> Sample size: 1000 complete cases
#> Estimand: ATE 
#> Variance method: bootstrap-full 
#> Bootstrap iterations: 50 
#> Trimming: None
#> 
#> Sample sizes used in Cox model:
#>  treatment   n events
#>          A 244    128
#>          B 188    133
#>          C 208    132
#>          D 360    275
#> 
#> Log hazard ratio estimates:
#>  treatment  logHR SE_robust SE_bootstrap
#>        Z:B 0.5371    0.1320       0.1294
#>        Z:C 0.1766    0.1287       0.1303
#>        Z:D 0.8116    0.1130       0.1106
#> 
#> Note: Use summary() for confidence intervals and exponentiated hazard ratios.

Estimand Options

Same as surveff(): ATE, ATT, and overlap weights are supported. Trimming options (symmetric/asymmetric) are also available.

# ATE (default)
marCoxph(..., estimand = "ATE")

# ATT targeting specific group
marCoxph(..., estimand = "ATT", att_group = "treated")

# Overlap weights
marCoxph(..., estimand = "overlap")

Trimming Options

Propensity score trimming works identically to surveff():

# Symmetric trimming
marCoxph(..., estimand = "ATE", trim = "symmetric", delta = 0.1)

# Asymmetric trimming
marCoxph(..., estimand = "ATE", trim = "asymmetric", alpha = 0.05)

Variance Options

Two variance estimation methods:

# Bootstrap variance (default)
marCoxph(..., variance_method = "bootstrap", B = 200)

# Robust sandwich variance (no bootstrap)
marCoxph(..., variance_method = "robust")

The boot_level argument ("full" or "strata") also applies here when using bootstrap variance, with the same interpretation as in surveff().

Summary Output

The summary() method provides detailed results with confidence intervals on both log and original scales:

summary(hr_multi, conf_level = 0.95)
#> 
#> === Marginal Cox Model Summary ===
#> 
#> Treatment variable: Z 
#> Treatment levels: A, B, C, D 
#> Reference group: A 
#> Number of groups: 4 
#> Sample size: 1000 complete cases
#> Estimand: ATE 
#> Variance method: bootstrap-full 
#> Bootstrap iterations: 50 
#> Confidence level: 0.95 
#> 
#> --- Sample Sizes in Cox Model ---
#> 
#>  treatment   n events
#>          A 244    128
#>          B 188    133
#>          C 208    132
#>          D 360    275
#> 
#> --- Log Hazard Ratios (log scale) ---
#> 
#>  treatment  logHR     SE 95%CI.lower 95%CI.upper
#>        Z:B 0.5371 0.1294      0.2835      0.7908
#>        Z:C 0.1766 0.1303     -0.0788      0.4320
#>        Z:D 0.8116 0.1106      0.5948      1.0285
#> 
#> --- Hazard Ratios (original scale) ---
#> 
#>  treatment     HR 95%CI.lower 95%CI.upper
#>        Z:B 1.7111      1.3277      2.2051
#>        Z:C 1.1932      0.9242      1.5403
#>        Z:D 2.2516      1.8126      2.7968

Output styles:

summ_hr <- summary(hr_multi, style = "returns")
names(summ_hr)
#>  [1] "logHR"            "logHR_CI_lower"   "logHR_CI_upper"   "SE"              
#>  [5] "HR"               "HR_CI_lower"      "HR_CI_upper"      "variance_method" 
#>  [9] "conf_level"       "n_per_group"      "events_per_group"
summ_hr$HR  # Exponentiated hazard ratios
#>      Z:B      Z:C      Z:D 
#> 1.711070 1.193173 2.251570

Summary

Feature surveff() marCoxph()
Output Survival curves, survival differences Hazard ratios
Estimands ATE, ATT, overlap ATE, ATT, overlap
Censoring adjustment Weibull or Cox None (standard Cox)
Variance Analytical or bootstrap Robust or bootstrap
Trimming Symmetric, asymmetric Symmetric, asymmetric
Multiple groups Yes (with contrast_matrix) Yes

References

Cheng, C., Li, F., Thomas, L. E., & Li, F. (2022). Addressing extreme propensity scores in estimating counterfactual survival functions via the overlap weights. American Journal of Epidemiology, 191(6), 1140-1151.

Li, F., & Li, F. (2019). Propensity score weighting for causal inference with multiple treatments. The Annals of Applied Statistics, 13(4), 2389-2415.