Workshop: Inverse Probability of Treatment Weighting (IPTW) in R

Author

Armando Teixeira-Pinto

Published

July 1, 2025

1 Quick Background

  • IPTW Overview: Inverse Probability of Treatment Weighting (IPTW) is a causal inference method used to estimate the average treatment effect (ATE) in observational data by creating a pseudo-population where confounders are balanced between treatment groups. It weights each observation by the inverse of the probability of receiving the observed treatment (propensity score, PS).
  • Assumptions: No unmeasured confounding, positivity (0 < PS < 1 for all), consistency, and correct PS model specification.
  • For Binary Exposure: Exposure A (1=treated, 0=untreated), binary outcome Y (1=event, 0=no). PS = P(A=1 | confounders). Weights: 1/PS for A=1, 1/(1-PS) for A=0.
  • Stabilization: Use stabilized weights (SW = P(A)/PS for A=1, P(A)/(1-PS) for A=0) to reduce variance.

We’ll simulate a clinical dataset on 500 patients: exposure = drug (binary: 1=yes, 0=no), outcome = recovery (binary: 1=yes, 0=no), confounders = age (continuous), gender (binary), severity (continuous). Confounding: older/male/higher severity less likely to get drug but affect recovery.

Note: Run the code in RStudio or your R environment. We’ll use packages like ipw for weights, cobalt for balance checks, and survey for weighted analysis.

2 Example: Two-Category Exposure

2.1 Install and Load Packages

We need ipw for IPTW, cobalt for balance diagnostics, survey for weighted models, and ggplot2 for plots.

# Install packages if not already installed
#install.packages(c("ipw", "cobalt", "survey", "ggplot2"))

# Load libraries
library(ipw)
library(cobalt)
library(survey)
library(ggplot2)
library(gtsummary)
  • What this does: Installs/loads tools. ipw computes weights; cobalt assesses balance; survey handles weighted regressions.

2.2 Simulate the Clinical Dataset

We’ll create data with confounding: drug less likely for older, male, severe patients; recovery lower for older, severe.

# Set seed for reproducibility
set.seed(456)

# Simulate 500 patients
n <- 500
age <- rnorm(n, mean = 60, sd = 10)  # Age ~ N(60, 10)
gender <- rbinom(n, 1, 0.5)  # Gender: 0=female, 1=male
severity <- rnorm(n, mean = 5, sd = 2)  # Severity score ~ N(5, 2)

# Propensity: lower for older, male, severe
logit_ps <- -.1 - 0.02 * age - 0.5 * gender - 0.1 * severity
ps <- plogis(logit_ps)
drug <- rbinom(n, 1, ps)  # Exposure: drug (1=yes)

# Outcome: recovery higher with drug, lower with age/severity
logit_y <- 2.5 + 1.5 * drug  -0.03 * age  -0.2 * gender  -0.4 * severity
recovery <- rbinom(n, 1, plogis(logit_y))  # Outcome: recovery (1=yes)

# Combine into data frame
data <- data.frame(age, gender, severity, drug, recovery)

# View summary
tbl_summary(data)
Characteristic N = 5001
age 61 (54, 68)
gender 261 (52%)
severity 5.04 (3.81, 6.47)
drug 65 (13%)
recovery 129 (26%)
1 Median (Q1, Q3); n (%)
  • What this does: Generates confounded data. Approximately true ATE: drug increases recovery odds by exp(1.5) ≈ 4.48.
  • Expected output : Table with summaries

2.3 Estimate Propensity Scores

Fit logistic regression for PS: P(drug=1 | age, gender, severity).

# Fit PS model
ps_model <- glm(drug ~ age + gender + severity, 
                family = binomial, 
                data = data)
data$ps <- predict(ps_model, type = "response") #PS

# Summary of PS
summary(data$ps)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0194  0.0821  0.1185  0.1300  0.1646  0.4005 
ggplot(data, aes(x = ps, fill = factor(drug))) + 
    geom_density(alpha = 0.5) + 
  labs(title = "Propensity Score Distribution by Treatment")

  • What this does: Estimates PS. Check overlap (densities should overlap for positivity).
  • Expected output: PS range 0.1-0.6 (mean ~0.3). Plot shows overlap, treated have higher PS.

2.4 Compute IPTW Weights

Use stabilized weights for stability.

# Marginal P(drug=1)
p_drug <- mean(data$drug)

# Stabilized weights
data$sw <- ifelse(data$drug == 1, 
                  p_drug / data$ps, 
                  (1 - p_drug) / (1 - data$ps))

# Summary of weights
summary(data$sw)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3246  0.9347  0.9783  0.9987  1.0372  2.8843 
ggplot(data, aes(x = sw)) + 
    geom_histogram() + 
    labs(title = "Distribution of Stabilized Weights")

  • What this does: Computes SW. Trim if extreme (e.g., >10 or <0.1), but here none.
  • Expected output: Weights 0.5-2.0 (mean ~1). Histogram centered around 1.

2.5 Check Balance

Assess confounder balance pre/post weighting using standardized mean differences (SMD <0.1 ideal).

# Balance check with cobalt
bal <- bal.tab(drug ~ age + gender + severity, data = data, weights = data$sw, 
               method = "weighting", 
               estimand = "ATE")
print(bal)
Balance Measures
            Type Diff.Adj
age      Contin.  -0.0157
gender    Binary  -0.0256
severity Contin.   0.0005

Effective sample sizes
           Control Treated
Unadjusted  435.     65.  
Adjusted    432.36   51.82
# Love plot
love.plot(bal, threshold = 0.1, abs = TRUE)

  • What this does: bal.tab computes SMD; love.plot visualizes.
  • Expected output: Unweighted SMD >0.2 (confounded); weighted <0.05 (balanced).

2.6 Estimate Treatment Effect with Weighted Analysis

Use weighted logistic regression via survey for ATE on odds ratio scale (or risk difference with linear).

# Create survey design
svy_design <- svydesign(ids = ~1, 
                        weights = ~sw, 
                        data = data)

# Weighted GLM for OR
fit_weighted <- svyglm(recovery ~ drug, 
                       design = svy_design, 
                       family = binomial)
tbl_regression(fit_weighted, exponentiate = T)
Characteristic OR 95% CI p-value
drug 4.16 2.29, 7.58 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
  • What this does: Estimates ATE.
  • Expected output: OR ≈ 4.2 (95% CI 2.3-7.6 ), close to true 4.48

2.7 Compare to Unweighted Analysis (For Illustration)

Naive model ignores confounding.

# Unweighted GLM
fit_unweighted <- glm(recovery ~ drug, 
                      family = binomial, 
                      data = data)
tbl_regression(fit_unweighted, 
               exponentiate = T)
Characteristic OR 95% CI p-value
drug 5.32 3.10, 9.25 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
  • Expected: Biased OR due to confounding.

3 More Complex Example: Three-Category Exposure

Now, extend to multi-level exposure: drug_level (0=none, 1=low, 2=high). Use multinomial PS.

3.1 Simulate Data with Three-Level Exposure

# Set seed
set.seed(789)

# Confounders
age <- rnorm(n, mean = 60, sd = 10)
gender <- rbinom(n, 1, 0.5)
severity <- rnorm(n, mean = 5, sd = 2)

# Multinomial logits for PS (ref=0)
logit1 <- -1 + -0.01 * age -0.3 * gender -0.05 * severity  # For level 1
logit2 <- 1.5 + -0.03 * age -0.7 * gender -0.15 * severity  # For level 2
prob0 <- 1 / (1 + exp(logit1) + exp(logit2))
prob1 <- exp(logit1) * prob0
prob2 <- exp(logit2) * prob0

# Sample drug_level
drug_level <- sapply(1:n, function(i) sample(0:2, 1, prob = c(prob0[i], prob1[i], prob2[i])))

# Outcome: dose-response
logit_y <- 1 + 2.0 * (drug_level == 1) + 2.0 * (drug_level == 2) - 0.03 * age - 0.2 * gender - 0.4 * severity
recovery <- rbinom(n, 1, plogis(logit_y))

# Data frame
data_multi <- data.frame(age, gender, severity, 
                         drug_level = factor(drug_level), 
                         recovery)

tbl_summary(data_multi)
Characteristic N = 5001
age 60 (53, 67)
gender 252 (50%)
severity 4.98 (3.54, 6.15)
drug_level
    0 351 (70%)
    1 45 (9.0%)
    2 104 (21%)
recovery 73 (15%)
1 Median (Q1, Q3); n (%)
  • True effects: OR=exp(1)≈2.7 for low vs none, exp(2)≈7.4 for high vs none.

3.2 Estimate Generalized Propensity Scores

Use multinomial logistic via nnet.

library(nnet)
ps_model_multi <- multinom(drug_level ~ age + gender + severity, 
                           data = data_multi)
# weights:  15 (8 variable)
initial  value 549.306144 
iter  10 value 372.308039
final  value 372.247212 
converged
# Matrix of PS for each level
ps_multi <- predict(ps_model_multi, type = "probs") 

# Observed PS
data_multi$ps_obs <- sapply(1:n, 
                            function(i) ps_multi[i, as.character(data_multi$drug_level[i])])

ggplot(data_multi, aes(x = ps, fill = factor(drug_level))) + 
    geom_density(alpha = 0.5) + 
  labs(title = "Propensity Score Distribution by Treatment")

3.3 Compute Generalized IPTW Weights

Stabilized: SW = P(A=a) / PS(A=a | confounders)

p_drug_level <- prop.table(table(data_multi$drug_level))  # Marginal prob

data_multi$sw <- ifelse(data_multi$drug_level == 0, 
                        p_drug_level[1]/data_multi$ps_obs, 
                        ifelse(data_multi$drug_level == 1, 
                               p_drug_level[2]/data_multi$ps_obs,  
                                p_drug_level[3]/data_multi$ps_obs))

summary(data_multi$sw)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3431  0.8302  0.9164  1.0011  1.1061  4.2041 
ggplot(data_multi, aes(x = sw)) + 
    geom_histogram() + 
    labs(title = "Distribution of Stabilized Weights")

3.4 Check Balance

For multi-level, compare each pair or use generalized balance.

# Balance check with cobalt
bal_multi <- bal.tab(drug_level ~ age + gender + severity, 
                     data = data_multi, 
                     weights = data_multi$sw, 
                     method = "weighting", 
                     estimand = "ATE", 
                     which.treat = .all)
print(bal_multi)
Balance by treatment pair

 - - - 0 (0) vs. 1 (1) - - - 
Balance Measures
            Type Diff.Adj
age      Contin.  -0.0102
gender    Binary  -0.0006
severity Contin.  -0.0065

Effective sample sizes
                0     1
Unadjusted 351.   45.  
Adjusted   336.92 43.69

 - - - 0 (0) vs. 2 (1) - - - 
Balance Measures
            Type Diff.Adj
age      Contin.  -0.0686
gender    Binary   0.0112
severity Contin.   0.0368

Effective sample sizes
                0      2
Unadjusted 351.   104.  
Adjusted   336.92  72.29

 - - - 1 (0) vs. 2 (1) - - - 
Balance Measures
            Type Diff.Adj
age      Contin.  -0.0584
gender    Binary   0.0117
severity Contin.   0.0433

Effective sample sizes
               1      2
Unadjusted 45.   104.  
Adjusted   43.69  72.29
 - - - - - - - - - - - - - - - -  
# Love plot
love.plot(bal_multi, threshold = 0.1, abs = TRUE)

  • Expected: Weighted SMD <0.1 across comparisons.

3.5 Weighted Analysis

Weighted GLM with factor exposure.

svy_design_multi <- svydesign(ids = ~1, 
                              weights = ~sw, 
                              data = data_multi)

fit_weighted_multi <- svyglm(recovery ~ drug_level, 
                             design = svy_design_multi, 
                             family = binomial)

tbl_regression(fit_weighted_multi, exponentiate = T)
Characteristic OR 95% CI p-value
drug_level


    0 — —
    1 5.21 2.33, 11.7 <0.001
    2 7.09 3.64, 13.8 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
  • Expected: ORs close to … and …. for level 1 and 2 vs 0.

4 Workshop Exercises

  1. Trim weights at 1st/99th percentile—how does effect change?
  2. Add interactions in PS model.
  3. Use WeightIt package for alternative weighting.
  4. Apply to real data (e.g., from twang package examples).

This covers IPTW for binary and multi-level exposures! If questions, let me know.