---
title: "PSsurvival: Propensity Score Weighting for Survival Analysis"
output:
  html_document:
    toc: true
    toc_float:
      collapsed: true
      smooth_scroll: true
vignette: >
  %\VignetteIndexEntry{PSsurvival: Propensity Score Weighting for Survival Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

```{r load-package}
library(PSsurvival)
```

## 2. Data for Vignette

The package includes two simulated datasets:

```{r 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)
```

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

```{r surveff-basic}
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)
```

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

```{r 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")

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

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

```{r surveff-multi-trimmed}
# 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
)

print(result_multi_trimmed, max.len = 3)
```

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

```{r weight-methods, eval=FALSE}
# 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.

```{r trimming, eval=FALSE}
# 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.

```{r 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)
```

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

```{r variance-methods, eval=FALSE}
# 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.

```{r boot-level, eval=FALSE}
# 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

```{r summary-surveff}
summary(result_bin, conf_level = 0.95, max.len = 5)
```

For programmatic access:

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

### Plotting Results

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

```{r plot-surv, fig.width=7, fig.height=5}
# Survival curves for all groups
plot(result_multi, type = "surv")
```

```{r plot-survdiff, fig.width=7, fig.height=5}
# Treatment effect curves
plot(result_multi, type = "survdiff")
```

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

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

```{r plot-subset-diff, fig.width=7, fig.height=5}
# 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

```{r plot-custom, fig.width=7, fig.height=5}
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

```{r 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",
  weight_method = "OW"
)

print(hr_result)
```

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

### Multiple Treatment Groups

```{r 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",
  weight_method = "IPW",
  variance_method = "bootstrap",
  B = 50,
  seed = 456
)

print(hr_multi)
```

```{r marCoxph-multi-trimmed}
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)
```

### Weighting Methods

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

```{r marcoxph-weight-methods, eval=FALSE}
# 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()`:

```{r marcoxph-trimming, eval=FALSE}
# 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)

```{r marCoxph-variance, eval=FALSE}
# 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:

```{r summary-marcoxph}
summary(hr_multi, conf_level = 0.95)
```

**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

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

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

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

print(km_result)
```

#### With Propensity Score Weights

```{r weightedKM-OW}
# 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)
```

```{r weightedKM-IPW}
# 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)
```

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

```{r plot-weightedKM, fig.width=7, fig.height=5}
# Survival curves with log CI
plot(km_ow, type = "Kaplan-Meier", conf_type = "log")
```

```{r plot-weightedKM-CR, fig.width=7, fig.height=5}
# Cumulative risk curves with log-log CI
plot(km_ow, type = "CR", conf_type = "log-log")
```

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

```{r plot-weightedKM-risktable, fig.width=7, fig.height=6}
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"))
```

**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

```{r summary-weightedKM}
summary(km_ow, type = "Kaplan-Meier", conf_type = "log-log", print.rows = 5)
```

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

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