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

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)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## ----summary-weightedKM-------------------------------------------------------
summary(km_ow, type = "Kaplan-Meier", conf_type = "log-log", print.rows = 5)


.