Main user interface for estimating marginal hazard ratios using propensity score weighting. Supports binary and multiple treatment groups with various weighting schemes (IPW, OW, or ATT) and optional trimming. Variance can be estimated via bootstrap or robust sandwich estimator.
marCoxph(
data,
ps_formula,
time_var,
event_var,
reference_level,
weight_method = "IPW",
att_group = NULL,
trim = FALSE,
delta = NULL,
variance_method = "bootstrap",
boot_level = "full",
B = 100,
parallel = FALSE,
mc.cores = 2,
seed = NULL,
ps_control = list(),
robust = TRUE
)Object of class "marCoxph" containing:
Fitted coxph model object.
Named vector of estimated log hazard ratios. Names are formatted as "treatment_var:level" (e.g., "Z:B" for treatment Z, level B vs reference).
Named vector of robust standard errors from coxph.
Named vector of bootstrap standard errors. NULL if
variance_method = "robust".
Named vector of sample sizes per treatment group used in Cox model fitting (after trimming).
Named vector of event counts per treatment group used in Cox model fitting (after trimming).
Variance method used: "bootstrap-full", "bootstrap-strata", or "robust".
Target estimand used.
Target group for ATT (NULL if not applicable).
Trimming method (NULL if no trimming).
Symmetric trimming threshold (NULL if not applicable).
Asymmetric trimming threshold (NULL if not applicable).
Name of treatment variable.
Sorted unique treatment values.
Reference level used in Cox model.
Number of treatment groups.
Number of complete cases used in analysis.
Propensity score estimation results.
Weight estimation results.
Bootstrap results (NULL if variance_method = "robust").
Contains: boot_samples, boot_allocation, n_success_by_group, B.
Data frame containing treatment, survival outcome, and covariates.
Formula for propensity score model: treatment ~ covariates.
Character string specifying the time-to-event variable name.
Character string specifying the event indicator variable name. Should be coded as 1=event, 0=censored.
Treatment level to use as reference in Cox model. MANDATORY. Must be one of the treatment levels.
Weighting method: "IPW" (inverse probability weighting), "OW" (overlap weighting), or "ATT" (average treatment effect on the treated). Default "IPW".
Target group for ATT. Required if weight_method = "ATT".
Logical. Perform symmetric propensity score trimming? Default FALSE.
If TRUE, symmetric trimming is applied (Crump extension for multiple treatments).
See estimate_weights for trimming details. Ignored if
weight_method = "OW". Asymmetric trimming is no longer supported due to
poor statistical performance.
Threshold for symmetric trimming in \((0, 1/J]\), where \(J\) is the number
of treatment levels. Default NULL uses recommended values: 0.1 for binary
treatment, 0.067 for 3 groups, \(1/(2J)\) for \(J \ge 4\) (Yoshida et al., 2019).
Used only if trim = TRUE.
Variance estimation method: "bootstrap" (default) or "robust".
"bootstrap" resamples the entire analysis pipeline. "robust" uses the sandwich
variance estimator from coxph() without bootstrap.
Bootstrap sampling level: "full" (default) or "strata".
"full" resamples from entire dataset (standard for observational studies). "strata"
resamples within each treatment group preserving group sizes (useful when treatment assignment
follows a stratified or fixed-ratio design). Only used if variance_method = "bootstrap".
Number of bootstrap iterations. Default 100. Used only if variance_method = "bootstrap".
Logical. Use parallel bootstrap computation? Default FALSE.
Number of cores for parallel bootstrap. Default 2.
Random seed for bootstrap reproducibility. Default NULL.
Control parameters for propensity score model. Default list().
Logical. Use robust (sandwich) variance in Cox model fitting?
Default TRUE. When TRUE, coxph() is called with robust = TRUE.
**Weighting Methods:**
The weight_method parameter specifies the target population for causal
inference:
IPW (Inverse Probability Weighting): Observations are weighted by the inverse probability of their observed treatment, \(w_i = 1/e_j(X_i)\) where j is the observed treatment group. Inference targets the combined population (ATE type).
OW (Overlap Weighting): Observations are weighted by overlap weights, which extends to multiple treatment groups (Li et al., 2018; Li and Li, 2019). Inference targets the population at clinical equipoise (overlap population).
ATT (Average Treatment Effect on the Treated): IPW weights tilted toward a specified target group. Observations in the target group receive weight 1, others receive \(w_i = e_{\text{target}}(X_i) / e_j(X_i)\). Inference targets the specified treatment group population (ATT type).
**Analysis Workflow:**
1. Extract treatment variable from ps_formula.
2. Estimate propensity scores using multinomial logistic regression (or logistic
for binary treatment).
3. Calculate propensity score weights based on weight_method and optional trim.
4. Fit marginal Cox model Surv(time, event) ~ treatment with weights.
5. Estimate variance via bootstrap (resampling full pipeline) or robust sandwich
estimator.
**Variance Estimation:**
- bootstrap: Resamples data (full or stratified), re-estimates PS and weights,
re-fits Cox model. Provides bootstrap SE for log hazard ratios.
- robust: Uses robust sandwich variance from coxph() directly. No
bootstrap performed (faster but may be less accurate with extreme weights).
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.
Yoshida, K., et al. (2019). Multinomial extension of propensity score trimming methods: A simulation study. American Journal of Epidemiology, 188(3), 609-616.
# \donttest{
# Example 1: Binary treatment with overlap weighting
data(simdata_bin)
result1 <- marCoxph(
data = simdata_bin,
ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
time_var = "time",
event_var = "event",
reference_level = "A",
weight_method = "OW"
)
summary(result1)
# Example 2: Multiple treatments with ATT and robust variance
data(simdata_multi)
result2 <- marCoxph(
data = simdata_multi,
ps_formula = Z ~ X1 + X2 + X3 + B1 + B2,
time_var = "time",
event_var = "event",
reference_level = "C",
weight_method = "ATT",
att_group = "C",
variance_method = "robust"
)
summary(result2)
# }
Run the code above in your browser using DataLab