--- title: "Causal Inference with IPTW in riskdiff" author: "John D. Murphy" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{Causal Inference with IPTW in riskdiff} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, warning = FALSE, message = FALSE ) library(riskdiff) ``` ## Introduction This vignette demonstrates how to use the **riskdiff** package to perform causal inference using Inverse Probability of Treatment Weighting (IPTW). IPTW is a powerful method for estimating causal effects from observational data by creating a pseudo-population where treatment assignment is independent of measured confounders. ### When to Use IPTW IPTW is particularly useful when: - You want to estimate **causal effects** (not just associations) - You have measured all important confounders - The treatment assignment mechanism is complex - You want to estimate population-level effects (ATE) or effects on specific subgroups (ATT/ATC) ### Key Assumptions IPTW relies on three main assumptions: 1. **No unmeasured confounding**: All confounders are measured and included 2. **Positivity**: All subjects have non-zero probability of receiving either treatment 3. **Correct model specification**: The propensity score model is correctly specified ## Basic IPTW Analysis Let's start with a simple example using the Cachar cancer screening dataset: ```{r basic_example} # Load the data data(cachar_sample) # Quick look at the data head(cachar_sample) table(cachar_sample$areca_nut, cachar_sample$abnormal_screen) ``` ### Step 1: Calculate IPTW Weights First, we estimate propensity scores and calculate weights: ```{r calculate_weights} # Calculate ATE weights for areca nut use iptw_result <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking", "tobacco_chewing"), weight_type = "ATE", verbose = TRUE ) # Examine the results print(iptw_result) ``` ### Step 2: Check IPTW Assumptions Before interpreting results, we should check key assumptions: ```{r check_assumptions} # Comprehensive assumption checking assumptions <- check_iptw_assumptions(iptw_result, verbose = TRUE) ``` ### Step 3: Visualize Balance Visual diagnostics help assess whether weighting achieved balance: ```{r balance_plots, fig.width=8, fig.height=6} # Create balance plots (requires ggplot2) if (requireNamespace("ggplot2", quietly = TRUE)) { plots <- create_balance_plots(iptw_result, plot_type = "both") print(plots$love_plot) print(plots$ps_plot) } ``` ### Step 4: Estimate Causal Risk Difference Now we can estimate the causal effect: ```{r causal_effect} # Estimate ATE using IPTW rd_causal <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking", "tobacco_chewing"), weight_type = "ATE", verbose = TRUE ) print(rd_causal) summary(rd_causal) ``` ## Different Types of Causal Effects IPTW can estimate different causal estimands depending on the research question: ### Average Treatment Effect (ATE) The effect of treatment if the entire population received treatment vs. if none received treatment: ```{r ate_example} rd_ate <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking"), weight_type = "ATE" ) cat("ATE: The average causal effect of areca nut use in the population\n") cat("Risk Difference:", scales::percent(rd_ate$rd_iptw, accuracy = 0.01), "\n") ``` ### Average Treatment Effect on the Treated (ATT) The effect among those who actually received treatment: ```{r att_example} rd_att <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking"), weight_type = "ATT" ) cat("ATT: The average causal effect among areca nut users\n") cat("Risk Difference:", scales::percent(rd_att$rd_iptw, accuracy = 0.01), "\n") ``` ### Average Treatment Effect on the Controls (ATC) The effect among those who did not receive treatment: ```{r atc_example} rd_atc <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking"), weight_type = "ATC" ) cat("ATC: The average causal effect among non-users of areca nut\n") cat("Risk Difference:", scales::percent(rd_atc$rd_iptw, accuracy = 0.01), "\n") ``` ## Causal Number Needed to Treat For intervention planning, causal NNT estimates are often more actionable than risk differences: ```{r causal-nnt} # Calculate causal NNT for population (ATE) nnt_ate <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking"), weight_type = "ATE", nnt = TRUE ) print(nnt_ate) ``` This tells us that if we could eliminate areca nut use in the population, we would need to intervene on approximately r if(!is.infinite(nnt_ate$rd_iptw)) round(nnt_ate$rd_iptw) else "many" individuals to prevent one abnormal screening result. ### NNT for Different Causal Estimands Different weight types give different causal interpretations: ```{r nnt_different weights} # Among those who currently use areca nut (ATT) nnt_att <- calc_risk_diff_iptw( cachar_sample, "abnormal_screen", "areca_nut", c("age", "sex", "residence", "smoking"), weight_type = "ATT", nnt = TRUE ) cat("NNT among current users (ATT):", if(!is.infinite(nnt_att$rd_iptw)) round(nnt_att$rd_iptw) else "Undefined", "\n") # Among those who don't currently use areca nut (ATC) nnt_atc <- calc_risk_diff_iptw( cachar_sample, "abnormal_screen", "areca_nut", c("age", "sex", "residence", "smoking"), weight_type = "ATC", nnt = TRUE ) cat("NNT among current non-users (ATC):", if(!is.infinite(nnt_atc$rd_iptw)) round(nnt_atc$rd_iptw) else "Undefined", "\n") ``` ## Advanced Options ### Bootstrap Confidence Intervals For small samples or when assumptions are questionable, bootstrap confidence intervals may be more robust: ```{r bootstrap_example} rd_bootstrap <- calc_risk_diff_iptw( data = cachar_sample, outcome = "head_neck_abnormal", treatment = "tobacco_chewing", covariates = c("age", "sex", "residence", "areca_nut"), bootstrap_ci = TRUE, boot_n = 500, # Use more in practice (1000+) verbose = FALSE ) print(rd_bootstrap) ``` ### Different Propensity Score Models You can specify different link functions for the propensity score model: ```{r ps_models} # Logistic regression (default) ps_logit <- calc_iptw_weights( data = cachar_sample, treatment = "tobacco_chewing", covariates = c("age", "sex", "residence", "areca_nut"), method = "logistic" ) # Probit regression ps_probit <- calc_iptw_weights( data = cachar_sample, treatment = "tobacco_chewing", covariates = c("age", "sex", "residence", "areca_nut"), method = "probit" ) # Compare propensity score distributions cat("Logistic PS range:", round(range(ps_logit$ps), 3), "\n") cat("Probit PS range:", round(range(ps_probit$ps), 3), "\n") ``` ### Weight Stabilization and Trimming Stabilized weights often have better properties: ```{r stabilization} # Unstabilized weights ps_unstab <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex", "residence"), stabilize = FALSE ) # Stabilized weights (default) ps_stab <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex", "residence"), stabilize = TRUE ) cat("Unstabilized weight variance:", round(var(ps_unstab$weights), 2), "\n") cat("Stabilized weight variance:", round(var(ps_stab$weights), 2), "\n") ``` Extreme weights can be problematic and may need trimming: ```{r trimming} # Check for extreme weights summary(ps_stab$weights) # Trim at 1st and 99th percentiles ps_trimmed <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex", "residence"), trim_weights = TRUE, trim_quantiles = c(0.01, 0.99) ) cat("Original weight range:", round(range(ps_stab$weights), 2), "\n") cat("Trimmed weight range:", round(range(ps_trimmed$weights), 2), "\n") ``` treatment = "smoking", covariates = c("maternal_age", "race", "education"), trim_weights = TRUE, trim_quantiles = c(0.01, 0.99) ) cat("Original weight range:", round(range(ps_stab$weights), 2), "\n") cat("Trimmed weight range:", round(range(ps_trimmed$weights), 2), "\n") ``` ## Comparison with Traditional Regression Let's compare IPTW results with traditional regression adjustment: ```{r comparison} # Traditional regression-based risk difference rd_regression <- calc_risk_diff( data = cachar_sample, outcome = "abnormal_screen", exposure = "areca_nut", adjust_vars = c("age", "sex", "residence", "smoking"), link = "auto" ) # IPTW-based causal risk difference rd_iptw <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "smoking"), weight_type = "ATE" ) # Compare results comparison_table <- data.frame( Method = c("Regression Adjustment", "IPTW (ATE)"), Risk_Difference = scales::percent(c(rd_regression$rd, rd_iptw$rd_iptw), accuracy = 0.01), CI_Lower = scales::percent(c(rd_regression$ci_lower, rd_iptw$ci_lower), accuracy = 0.01), CI_Upper = scales::percent(c(rd_regression$ci_upper, rd_iptw$ci_upper), accuracy = 0.01), P_Value = sprintf("%.3f", c(rd_regression$p_value, rd_iptw$p_value)) ) print(comparison_table) ``` ## Troubleshooting Common Issues ### Poor Balance If covariates remain imbalanced after weighting: ```{r poor_balance, eval=FALSE} # Check which variables have poor balance assumptions <- check_iptw_assumptions(iptw_result) poor_balance_vars <- assumptions$balance$poor_balance_vars if (length(poor_balance_vars) > 0) { cat("Variables with poor balance:", paste(poor_balance_vars, collapse = ", "), "\n") # Try including interactions or polynomial terms iptw_improved <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "I(age^2)", "sex", "residence", "smoking", "age:sex"), # Add interactions weight_type = "ATE" ) } ``` ### Extreme Propensity Scores When subjects have very high or low propensity scores: ```{r extreme_ps, eval=FALSE} # Check propensity score distribution iptw_result <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex", "residence") ) # Identify subjects with extreme scores extreme_low <- which(iptw_result$ps < 0.05) extreme_high <- which(iptw_result$ps > 0.95) if (length(extreme_low) > 0 || length(extreme_high) > 0) { cat("Consider trimming sample to region of common support\n") # Restrict to common support common_support <- iptw_result$ps >= 0.05 & iptw_result$ps <= 0.95 data_restricted <- cachar_sample[common_support, ] # Re-analyze with restricted sample rd_restricted <- calc_risk_diff_iptw( data = data_restricted, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence") ) } ``` ### Model Specification Testing different model specifications: ```{r model_spec, eval=FALSE} # Simple model ps_simple <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "sex") ) # Complex model with interactions ps_complex <- calc_iptw_weights( data = cachar_sample, treatment = "areca_nut", covariates = c("age", "I(age^2)", "sex", "residence", "smoking", "tobacco_chewing", "age:sex") ) # Compare balance check_iptw_assumptions(ps_simple, verbose = FALSE) check_iptw_assumptions(ps_complex, verbose = FALSE) ``` ## Sensitivity Analysis IPTW assumes no unmeasured confounding. Sensitivity analysis can assess robustness: ```{r sensitivity, eval=FALSE} # Simulate an unmeasured confounder set.seed(123) cachar_sample$unmeasured_confounder <- rbinom(nrow(cachar_sample), 1, 0.3) # Compare results with and without the unmeasured confounder rd_without_u <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence") ) rd_with_u <- calc_risk_diff_iptw( data = cachar_sample, outcome = "abnormal_screen", treatment = "areca_nut", covariates = c("age", "sex", "residence", "unmeasured_confounder") ) cat("Without unmeasured confounder:", scales::percent(rd_without_u$rd_iptw), "\n") cat("With unmeasured confounder:", scales::percent(rd_with_u$rd_iptw), "\n") cat("Difference:", scales::percent(abs(rd_without_u$rd_iptw - rd_with_u$rd_iptw)), "\n") ``` ## Best Practices 1. **Always check assumptions** before interpreting results 2. **Visualize balance** using love plots and propensity score distributions 3. **Consider the estimand** (ATE vs. ATT vs. ATC) based on your research question 4. **Use stabilized weights** unless you have a specific reason not to 5. **Trim extreme weights** cautiously - they may indicate model misspecification 6. **Compare with regression adjustment** as a sensitivity check 7. **Report effective sample size** to assess precision loss 8. **Consider bootstrap CIs** for small samples or non-normal distributions ## Reporting IPTW Results When reporting IPTW analyses, include: 1. **Propensity score model specification** and diagnostics 2. **Balance assessment** before and after weighting 3. **Weight distribution** and any trimming performed 4. **Effective sample size** and precision considerations 5. **Sensitivity analyses** when possible 6. **Clear statement of estimand** (ATE/ATT/ATC) ## Conclusion IPTW provides a powerful framework for causal inference from observational data. The **riskdiff** package makes these methods accessible while providing essential diagnostics and visualizations. Remember that causal inference requires careful thought about study design, confounders, and assumptions - the methods are only as good as these foundational elements. For more advanced applications, consider methods like: - Marginal structural models for time-varying treatments - Doubly robust estimation combining IPTW with outcome modeling - Machine learning approaches for propensity score estimation ## References Austin, P. C. (2011). An introduction to propensity score methods for reducing the effects of confounding in observational studies. *Multivariate Behavioral Research*, 46(3), 399-424. Hernán, M. A., & Robins, J. M. (2020). *Causal inference: What if*. Boca Raton: Chapman & Hall/CRC. Robins, J. M., Hernán, M. A., & Brumback, B. (2000). Marginal structural models and causal inference in epidemiology. *Epidemiology*, 11(5), 550-560.