# PSsurvival: Propensity Score Weighting for Survival Analysis

## 1. Introduction

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

- `surveff()`: Estimates counterfactual survival functions and survival differences
- `marCoxph()`: Estimates marginal hazard ratios via weighted Cox regression
- `weightedKM()`: Estimates weighted Kaplan-Meier and cumulative risk curves

All functions support binary and multiple treatment groups, multiple weighting schemes (IPW, overlap weights, ATT), and propensity score trimming.

    library(PSsurvival)

## 2. Data for Vignette

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:

- `X1`, `X2`, `X3`: Continuous covariates
- `B1`, `B2`: Binary covariates
- `Z`: Treatment group
- `time`: Follow-up time (range 0-20)
- `event`: Event indicator (1 = event, 0 = censored)

## 3. `surveff()`: Counterfactual Survival Probabilities and Differences

### Key Arguments

Data and model specification:

- `data`: Data frame containing treatment, outcome, and covariates
- `ps_formula`: Propensity score model (`treatment ~ covariates`)
- `censoring_formula`: Censoring model (`Surv(time, event) ~ covariates`)
- `eval_times`: Time points for evaluation (default: all unique event times)
- `contrast_matrix`: Matrix for pairwise comparisons (multiple groups only; column names must match treatment levels)

Weighting and trimming:

- `weight_method`: Weighting method (“IPW”, “OW”, or “ATT”)
- `att_group`: Target group for ATT (required if `weight_method = "ATT"`)
- `trim`: Logical, whether to apply symmetric trimming (default: FALSE)
- `delta`: Threshold for symmetric trimming (default: `NULL` for automatic selection)

Variance estimation:

- `censoring_method`: Method for censoring adjustment (“weibull” or “cox”)
- `variance_method`: Variance estimation (“analytical” or “bootstrap”)
- `boot_level`: Bootstrap sampling level (“full” or “strata”)
- `B`: Number of bootstrap iterations
- `parallel`, `mc.cores`, `seed`: Parallel computing and reproducibility

### 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,
      weight_method = "OW",
      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` specifies the censoring models (Weibull or Cox) within each treatment group for estimating the censoring scores.

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

    result_multi <- surveff(
      data = simdata_multi,
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      censoring_formula = survival::Surv(time, event) ~ X1 + B1,
      weight_method = "IPW",
      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).

    # with symmetric trimming
    result_multi_trimmed <- surveff(
      data = simdata_multi,
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      censoring_formula = survival::Surv(time, event) ~ X1 + B1,
      weight_method = "IPW",
      trim = TRUE,
      delta = 0.15,
      contrast_matrix = contrast_mat,
      censoring_method = "cox",
      variance_method = "bootstrap",
      B = 50,
      seed = 123
    )
    #> Excluding 251 trimmed observations (weights == 0) from calculations.

    print(result_multi_trimmed, 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: 491 time points
    #>
    #> Survival function estimates:
    #>           A      B C D
    #> [1,] 1.0000 0.9928 1 1
    #> [2,] 0.9949 0.9928 1 1
    #> [3,] 0.9896 0.9928 1 1
    #>   ... (488 more rows not shown)
    #>
    #> Treatment effect estimates:
    #>       A vs B  A vs C  A vs D
    #> [1,]  0.0072  0.0000  0.0000
    #> [2,]  0.0020 -0.0051 -0.0051
    #> [3,] -0.0032 -0.0104 -0.0104
    #>   ... (488 more rows not shown)

### Weighting Methods

Three weighting methods are supported:

- IPW: Inverse probability (of treatment) weighting, targets average treatment effect (ATE) for the entire population
- ATT: Propensity score weighting tilted to a specific treatment group (requires `att_group`)
- OW: Overlap weighting, targets the overlap population at clinical equipoise where treatment groups have most similar covariates

    # IPW (targets ATE)
    surveff(..., weight_method = "IPW")

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

    # Overlap weights
    surveff(..., weight_method = "OW")

### Propensity Score Trimming

Symmetric trimming excludes observations with extreme propensity scores. Not available with overlap weights (which already downweight extremes). The Crump extension is used for multiple treatments.

    # Symmetric trimming with default delta (automatic selection)
    surveff(..., weight_method = "IPW", trim = TRUE)

    # Symmetric trimming with custom delta
    surveff(..., weight_method = "IPW", trim = TRUE, delta = 0.1)

### Censoring Methods

Two methods are available for estimating censoring distributions:

- `"weibull"` (default): Parametric Weibull AFT model. Efficient when censoring follows Weibull distribution. Allows analytical variance for binary treatment.
- `"cox"`: Semi-parametric Cox proportional hazards model. More flexible but requires bootstrap variance.

    # 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"`: M-estimation theory-based variance. Only available for binary treatment with Weibull censoring. Fast and efficient when applicable.
- `"bootstrap"`: Resamples the entire estimation pipeline. Required for Cox censoring or multiple treatment groups. More computationally intensive but broadly applicable.

    # 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.

- `"full"` (default): Resamples from entire dataset. Standard choice for observational studies.
- `"strata"`: Resamples within treatment groups, preserving group sizes. Useful when treatment assignment follows a stratified or fixed-ratio design.

    # 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:

- `conf_level`: Confidence level (default 0.95)
- `max.len`: Maximum rows to display
- `round.digits`: Number of decimal places

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

- `survival_summary`: List of survival estimates by group
- `difference_summary`: List of treatment effect estimates

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

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

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",
         strata_to_plot = c("A", "C"))

    # Only specific contrasts (names must match contrast_matrix rownames exactly)
    plot(result_multi, type = "survdiff",
         strata_to_plot = c("A vs B", "A vs D"))

Customization options:

- `strata_to_plot`: Vector of strata names to display
- `strata_colors`: Custom colors for each stratum (length must match `strata_to_plot`)
- `max_time`: Maximum time on x-axis
- `include_CI`: Show confidence interval ribbons (TRUE/FALSE)
- `conf_level`: Confidence level for intervals
- `curve_width`: Line width for curves
- `CI_alpha`: Transparency of CI ribbons (0-1)
- `legend_position`: “right” or “bottom”
- `legend_title`, `plot_title`: Custom titles

    plot(result_multi,
         type = "surv",
         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")

## 4. `marCoxph()`: Marginal Hazard Ratios by Cox-PH model

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

### Key Arguments

Data and model specification:

- `data`: Data frame containing treatment, outcome, and covariates
- `ps_formula`: Propensity score model (`treatment ~ covariates`)
- `time_var`: Name of time-to-event variable (character string)
- `event_var`: Name of event indicator variable (character string, 1=event, 0=censored)
- `reference_level`: Reference treatment level in Cox model (MANDATORY)

Weighting and trimming:

- `weight_method`: Weighting method (“IPW”, “OW”, or “ATT”)
- `att_group`: Target group for ATT (required if `weight_method = "ATT"`)
- `trim`: Logical, whether to apply symmetric trimming (default: FALSE)
- `delta`: Threshold for symmetric trimming (default: `NULL` for automatic selection)

Variance estimation:

- `variance_method`: “bootstrap” (default) or “robust”
- `boot_level`: Bootstrap sampling level (“full” or “strata”)
- `B`: Number of bootstrap iterations
- `parallel`, `mc.cores`, `seed`: Parallel computing and reproducibility
- `robust`: Use robust variance in Cox model (default TRUE)

### 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",
      weight_method = "OW"
    )

    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",
      weight_method = "IPW",
      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.

    hr_multi_trimmed <- marCoxph(
      data = simdata_multi,
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      time_var = "time",
      event_var = "event",
      reference_level = "A",
      weight_method = "IPW",
      trim = TRUE,
      delta = 0.15,
      variance_method = "bootstrap",
      B = 50,
      seed = 456
    )

    print(hr_multi_trimmed)
    #>
    #> === 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: symmetric (delta = 0.15)
    #>
    #> Sample sizes used in Cox model:
    #>  treatment   n events
    #>          A 190     93
    #>          B 142     99
    #>          C 156     97
    #>          D 261    203
    #>
    #> Log hazard ratio estimates:
    #>  treatment  logHR SE_robust SE_bootstrap
    #>        Z:B 0.5637    0.1489       0.1879
    #>        Z:C 0.2315    0.1523       0.2169
    #>        Z:D 0.8309    0.1340       0.1788
    #>
    #> Note: Use summary() for confidence intervals and exponentiated hazard ratios.

### Weighting Methods

Same as `surveff()`: IPW, ATT, and OW are supported.

    # IPW (targets ATE)
    marCoxph(..., weight_method = "IPW")

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

    # Overlap weights
    marCoxph(..., weight_method = "OW")

### Trimming Options

Propensity score trimming works identically to `surveff()`:

    # Symmetric trimming with default delta
    marCoxph(..., weight_method = "IPW", trim = TRUE)

    # Symmetric trimming with custom delta
    marCoxph(..., weight_method = "IPW", trim = TRUE, delta = 0.1)

### Variance Options

Two variance estimation methods:

- `"bootstrap"` (default): Resamples the entire analysis pipeline (propensity scores, weights, and Cox model)
- `"robust"`: Uses sandwich variance from `coxph()` directly (faster, no resampling)

    # 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:

- `style = "prints"` (default): Displays formatted tables with log hazard ratios and exponentiated hazard ratios
- `style = "returns"`: Returns a list containing:
  - `logHR`, `logHR_CI_lower`, `logHR_CI_upper`: Log-scale estimates and CIs
  - `HR`, `HR_CI_lower`, `HR_CI_upper`: Original-scale estimates and CIs
  - `SE`: Standard errors for log-scale hazard ratio (from specified variance method)
  - `n_per_group`, `events_per_group`: Sample sizes and event counts

    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

## 5. `weightedKM()`: Kaplan-Meier or Cumulative Risk Curves with Propensity Score Weights

`weightedKM()` estimates the survival probabilities at each observed time points, using weighted Kaplan-Meier (KM) estimator with propensity score weights. It uses the classic weighted Greenwood variance formula. `plot.weightedKM()` creates weighted KM or Cumulative Risk (CR) curves for visualization.

### Key Arguments

Data specification:

- `data`: Data frame
- `treatment_var`: Name of treatment variable (character string)
- `time_var`: Name of time variable (character string)
- `event_var`: Name of event indicator (character string, 1=event, 0=censored)

Weighting:

- `weight_method`: “IPW”, “OW”, “ATT”, or “none” (classical unweighted KM)
- `ps_formula`: Propensity score model (required unless `weight_method = "none"`)
- `att_group`: Target group for ATT (required if `weight_method = "ATT"`)
- `trim`: Logical, whether to apply symmetric trimming (default: FALSE)
- `delta`: Threshold for symmetric trimming

### Fit a `weightedKM()` Object

#### Without Propensity Score Weights

    # Classical unweighted KM
    km_result <- weightedKM(
      data = simdata_bin,
      treatment_var = "Z",
      time_var = "time",
      event_var = "event",
      weight_method = "none"
    )

    print(km_result)
    #>
    #> ===============================================================================
    #> Weighted Kaplan-Meier Survival Estimates
    #> ===============================================================================
    #>
    #> Weight Method: none
    #> Treatment Levels: A, B
    #> Number of Complete Cases: 1000
    #> Number of Unique Event Times: 666
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: A
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.049    0.998 0.002
    #>  [2,] 0.066    0.995 0.003
    #>  [3,] 0.188    0.995 0.003
    #>  [4,] 0.312    0.995 0.003
    #>  [5,] 0.320    0.995 0.003
    #>  [6,] 0.396    0.993 0.004
    #>  [7,] 0.500    0.993 0.004
    #>  [8,] 0.520    0.993 0.004
    #>  [9,] 0.534    0.993 0.004
    #> [10,] 0.562    0.990 0.005
    #> ... (656 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: B
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.049    1.000 0.000
    #>  [2,] 0.066    1.000 0.000
    #>  [3,] 0.188    0.998 0.002
    #>  [4,] 0.312    0.997 0.002
    #>  [5,] 0.320    0.995 0.003
    #>  [6,] 0.396    0.995 0.003
    #>  [7,] 0.500    0.993 0.003
    #>  [8,] 0.520    0.991 0.004
    #>  [9,] 0.534    0.990 0.004
    #> [10,] 0.562    0.990 0.004
    #> ... (656 more rows)
    #>
    #> ===============================================================================
    #> Use summary() for confidence intervals. Use plot() to visualize.
    #> ===============================================================================

#### With Propensity Score Weights

    # Overlap-weighted KM
    km_ow <- weightedKM(
      data = simdata_bin,
      treatment_var = "Z",
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      time_var = "time",
      event_var = "event",
      weight_method = "OW"
    )

    print(km_ow)
    #>
    #> ===============================================================================
    #> Weighted Kaplan-Meier Survival Estimates
    #> ===============================================================================
    #>
    #> Weight Method: OW
    #> Treatment Levels: A, B
    #> Number of Complete Cases: 1000
    #> Number of Unique Event Times: 666
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: A
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.049    0.998 0.003
    #>  [2,] 0.066    0.996 0.004
    #>  [3,] 0.188    0.996 0.004
    #>  [4,] 0.312    0.996 0.004
    #>  [5,] 0.320    0.996 0.004
    #>  [6,] 0.396    0.993 0.005
    #>  [7,] 0.500    0.993 0.005
    #>  [8,] 0.520    0.993 0.005
    #>  [9,] 0.534    0.993 0.005
    #> [10,] 0.562    0.992 0.006
    #> ... (656 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: B
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.049    1.000 0.000
    #>  [2,] 0.066    1.000 0.000
    #>  [3,] 0.188    0.997 0.003
    #>  [4,] 0.312    0.995 0.005
    #>  [5,] 0.320    0.994 0.005
    #>  [6,] 0.396    0.994 0.005
    #>  [7,] 0.500    0.992 0.006
    #>  [8,] 0.520    0.989 0.007
    #>  [9,] 0.534    0.988 0.007
    #> [10,] 0.562    0.988 0.007
    #> ... (656 more rows)
    #>
    #> ===============================================================================
    #> Use summary() for confidence intervals. Use plot() to visualize.
    #> ===============================================================================

    # IPW-weighted KM with trimming
    km_ipw <- weightedKM(
      data = simdata_multi,
      treatment_var = "Z",
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      time_var = "time",
      event_var = "event",
      weight_method = "IPW",
      trim = TRUE,
      delta = 0.1
    )

    print(km_ipw)
    #>
    #> ===============================================================================
    #> Weighted Kaplan-Meier Survival Estimates
    #> ===============================================================================
    #>
    #> Weight Method: IPW
    #> Trimming: Symmetric (delta = 0.1 )
    #> Treatment Levels: A, B, C, D
    #> Number of Complete Cases: 1000
    #> Number of Unique Event Times: 667
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: A
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.036    1.000 0.000
    #>  [2,] 0.051    0.996 0.002
    #>  [3,] 0.070    0.992 0.003
    #>  [4,] 0.125    0.992 0.003
    #>  [5,] 0.243    0.992 0.003
    #>  [6,] 0.284    0.992 0.003
    #>  [7,] 0.287    0.992 0.003
    #>  [8,] 0.290    0.992 0.003
    #>  [9,] 0.301    0.992 0.003
    #> [10,] 0.324    0.988 0.003
    #> ... (657 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: B
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.036    0.994 0.002
    #>  [2,] 0.051    0.994 0.002
    #>  [3,] 0.070    0.994 0.002
    #>  [4,] 0.125    0.994 0.002
    #>  [5,] 0.243    0.994 0.002
    #>  [6,] 0.284    0.994 0.002
    #>  [7,] 0.287    0.988 0.003
    #>  [8,] 0.290    0.983 0.004
    #>  [9,] 0.301    0.983 0.004
    #> [10,] 0.324    0.983 0.004
    #> ... (657 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: C
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.036    1.000 0.000
    #>  [2,] 0.051    1.000 0.000
    #>  [3,] 0.070    1.000 0.000
    #>  [4,] 0.125    1.000 0.000
    #>  [5,] 0.243    1.000 0.000
    #>  [6,] 0.284    0.995 0.002
    #>  [7,] 0.287    0.995 0.002
    #>  [8,] 0.290    0.995 0.002
    #>  [9,] 0.301    0.992 0.003
    #> [10,] 0.324    0.992 0.003
    #> ... (657 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: D
    #> -------------------------------------------------------------------------------
    #>        time survival    se
    #>  [1,] 0.036    1.000 0.000
    #>  [2,] 0.051    1.000 0.000
    #>  [3,] 0.070    1.000 0.000
    #>  [4,] 0.125    0.997 0.002
    #>  [5,] 0.243    0.994 0.002
    #>  [6,] 0.284    0.994 0.002
    #>  [7,] 0.287    0.994 0.002
    #>  [8,] 0.290    0.994 0.002
    #>  [9,] 0.301    0.994 0.002
    #> [10,] 0.324    0.994 0.002
    #> ... (657 more rows)
    #>
    #> ===============================================================================
    #> Use summary() for confidence intervals. Use plot() to visualize.
    #> ===============================================================================

### Creating Weighted KM or CR Curves

`plot.weightedKM()` function visualizes a `weightedKM` object by survival (aka. Kaplan-Meier) curves or cumulative risk (aka. cumulative incidence) curves. By default, it produces curves with 0.95 confidence band and a risk table for all groups, with each group indicated by a specific color. It allows users to flexibly adjust the confidence band, risk table, colors, and other aesthetic components.

#### Kaplan-Meier vs. CR

Specified by `type`: `"Kaplan-Meier"` for Kaplan-Meier and `"CR"` for cumulative risk. Cumulative risk is defined as 1-KM.

#### Confidence intervals

Plot survival curves or cumulative risk curves with or without confidence intervals:

    # Survival curves with log CI
    plot(km_ow, type = "Kaplan-Meier", conf_type = "log")

    # Cumulative risk curves with log-log CI
    plot(km_ow, type = "CR", conf_type = "log-log")
    #> Warning: Removed 2 rows containing missing values or values outside the scale range
    #> (`geom_ribbon()`).

With or Without confidence intervals: Specified by `include_CI` (default: TRUE).

Confidence interval types: Three methods available via `conf_type`:

- `"plain"`: Normal approximation (may exceed [0,1])
- `"log"`: Log transformation (guarantees bounds in (0,1))
- `"log-log"`: Complementary log-log transformation (best for extreme probabilities, default)

The confidence intervals and corresponding transformation are calculated and performed on the scale of \(\hat{S}(t)\) for Kaplan-Meier and of \(1-\hat{S}(t)\) for cumulative risk, where \(\hat{S}(t)\) is the estimated survival probability at time \(t\).

#### Risk Tables

Display weighted number at risk and cumulative events below the curves:

    plot(km_ipw,
         risk_table = TRUE,
         risk_table_breaks = c(0, 5, 10, 15),
         risk_table_height = 0.3,
         risk_table_stats = c("n.risk", "n.acc.event"))
    #> Warning: Removed 9 rows containing missing values or values outside the scale range
    #> (`geom_ribbon()`).

Risk table customization:

- `risk_table`: Logical, add risk table (default FALSE)
- `risk_table_stats`: Statistics to display: “n.risk” (number at risk), “n.acc.event” (cumulative events)
- `risk_table_height`: Relative height of risk table (0-1, default 0.25)
- `risk_table_breaks`: Time points for columns (default: automatic)
- `risk_table_fontsize`: Font size (default 3.5)

### Summary Output

    summary(km_ow, type = "Kaplan-Meier", conf_type = "log-log", print.rows = 5)
    #>
    #> ===============================================================================
    #> Summary of Weighted Kaplan-Meier Estimates
    #> ===============================================================================
    #>
    #> Type: Survival Probability
    #> Confidence Level: 95%
    #> Confidence Type: log-log
    #> Weight Method: OW
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: A
    #> -------------------------------------------------------------------------------
    #>       time estimate    se CI_lower CI_upper
    #> [1,] 0.000    1.000 0.000    1.000    1.000
    #> [2,] 0.049    0.998 0.003    1.000    0.965
    #> [3,] 0.066    0.996 0.004    0.999    0.969
    #> [4,] 0.188    0.996 0.004    0.999    0.969
    #> [5,] 0.312    0.996 0.004    0.999    0.969
    #> ... (662 more rows)
    #>
    #> -------------------------------------------------------------------------------
    #> Treatment Group: B
    #> -------------------------------------------------------------------------------
    #>       time estimate    se CI_lower CI_upper
    #> [1,] 0.000    1.000 0.000    1.000    1.000
    #> [2,] 0.049    1.000 0.000       NA       NA
    #> [3,] 0.066    1.000 0.000       NA       NA
    #> [4,] 0.188    0.997 0.003    1.000    0.967
    #> [5,] 0.312    0.995 0.005    0.999    0.969
    #> ... (662 more rows)
    #>
    #> ===============================================================================
    #> Note: Full-precision results returned invisibly. Use print.rows to show more.
    #> ===============================================================================

The summary returns a list of matrices (one per treatment group) with full-precision estimates, invisibly. Customize display with `print.digits` and `print.rows`.

## 6. Multiple imputation with `mice`

Missing covariate data can be handled by multiple imputation. `surveff()` and `marCoxph()` provide `tidy()` and `df.residual()` methods, so they slot into the standard `mice` workflow — impute, fit each dataset with `with()`, combine with `pool()` — exactly like an ordinary model such as `coxph()`.

    library(mice)

    d <- simdata_bin
    set.seed(1)
    d$X1[sample(nrow(d), 100)] <- NA          # some missing covariate values
    imp <- mice(d, m = 5, printFlag = FALSE)  # impute

Marginal hazard ratio:

    fit <- with(imp, marCoxph(
      data = data.frame(Z, X1, X2, X3, B1, B2, time, event),
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      time_var = "time", event_var = "event",
      reference_level = "A", weight_method = "OW", variance_method = "robust"))
    summary(pool(fit))
    #>   term  estimate  std.error statistic       df      p.value
    #> 1  Z:B 0.5236878 0.08433858  6.209351 656.9997 9.432952e-10

Counterfactual survival is identical; pass a fixed `eval_times` so the pooled parameters line up across imputations:

    fit <- with(imp, surveff(
      data = data.frame(Z, X1, X2, X3, B1, B2, time, event),
      ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
      censoring_formula = survival::Surv(time, event) ~ X1 + B1,
      eval_times = c(0.5, 1, 1.5), weight_method = "OW", censoring_method = "weibull"))
    summary(pool(fit))
    #>                      term     estimate   std.error   statistic         df
    #> 1        survival:A@t=0.5  0.993414533 0.003816406 260.3010410  141226818
    #> 2          survival:A@t=1  0.987709181 0.005030314 196.3513874    6764865
    #> 3        survival:A@t=1.5  0.977930071 0.007098604 137.7637224   25162727
    #> 4        survival:B@t=0.5  0.991923194 0.004099385 241.9687959  780072754
    #> 5          survival:B@t=1  0.966869364 0.008018972 120.5727367 8929482505
    #> 6        survival:B@t=1.5  0.934314781 0.011010260  84.8585547    3211406
    #> 7 difference:B vs A@t=0.5 -0.001491338 0.005599100  -0.2663533   98684433
    #> 8   difference:B vs A@t=1 -0.020839817 0.009465539  -2.2016514   44685408
    #> 9 difference:B vs A@t=1.5 -0.043615290 0.013098272  -3.3298506    7003082
    #>        p.value
    #> 1 0.0000000000
    #> 2 0.0000000000
    #> 3 0.0000000000
    #> 4 0.0000000000
    #> 5 0.0000000000
    #> 6 0.0000000000
    #> 7 0.7899671567
    #> 8 0.0276899456
    #> 9 0.0008689304

The same pattern applies to multi-group data (`simdata_multi`).

## Summary

Feature | `surveff()` | `marCoxph()` | `weightedKM()`
--- | --- | --- | ---
Output | Survival curves, survival differences | Hazard ratios | Weighted KM curves
Weighting | IPW, ATT, OW | IPW, ATT, OW | IPW, ATT, OW, none
Censoring adjustment | Weibull or Cox | None (standard Cox) | None
Variance | Analytical or bootstrap | Robust or bootstrap | Weighted Greenwood
Trimming | Symmetric | Symmetric | Symmetric
Multiple groups | Yes (with contrast_matrix) | Yes | Yes
Risk tables | No | No | Yes

## References

Li, F., Morgan, K. L., & Zaslavsky, A. M. (2018). Balancing covariates via propensity score weighting. Journal of the American Statistical Association, 113(521), 390-400.

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

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.