In Lecture 1, we compared CVD risk for smokers and non-smokers using a simple difference in observed risks.
In Lecture 2, we ask a more careful question:
Can we compare smokers and non-smokers after adjusting for baseline covariates?
Our treatment is:
| Symbol | Meaning in this lecture |
|---|---|
T = 1 |
Smoker |
T = 0 |
Non-smoker |
Our outcome is:
| Symbol | Meaning in this lecture |
|---|---|
Y = 1 |
Had CVD before the interview |
Y = 0 |
Did not have CVD before the interview |
Our baseline covariates are:
| Symbol | Meaning in this lecture |
|---|---|
X |
Age, gender, race, education, income-to-poverty ratio, and BMI |
We will again keep the lecture simple and not use survey weights. The goal is to understand inverse probability weighting, not the full NHANES survey design.
Smokers and non-smokers might be different before we even look at CVD. For example, they might have different ages, gender distributions, education levels, or income levels.
Those baseline covariates can be related to both smoking and CVD. When a covariate is related to both the treatment and the outcome, it can act like a confounder.
The orange arrow is the causal question we care about: smoking to
CVD. The two black arrows from X remind us that baseline
covariates can affect who smokes and who develops CVD.
IPW stands for inverse probability weighting.
The idea is to give each person a weight based on how surprising their actual smoking status was, given their baseline covariates.
If a person was unlikely to be a smoker but actually was a smoker, that person gets a larger weight in the smoker group. If a person was unlikely to be a non-smoker but actually was a non-smoker, that person gets a larger weight in the non-smoker group.
This is why IPW is necessary: if the smoker and non-smoker groups have different baseline covariate mixes, then a simple CVD comparison can mix together the effect of smoking and the effect of those covariates.
data_path <- file.path("..", "Data", "NHANES_data.csv")
if (!file.exists(data_path)) {
data_path <- file.path("Data", "NHANES_data.csv")
}
nhanes <- read.csv(data_path, stringsAsFactors = FALSE)
baseline_vars <- c(
"age_yr",
"gender",
"race",
"educ_lvl",
"inc_to_pov_ratio",
"bmi"
)
analysis_vars <- c("smoker_indicator", "cvd_indicator", baseline_vars)
ipw_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]
ipw_data$smoker_indicator <- as.integer(ipw_data$smoker_indicator)
ipw_data$cvd_indicator <- as.integer(ipw_data$cvd_indicator)
ipw_data$smoking_group <- ifelse(
ipw_data$smoker_indicator == 1,
"Smoker",
"Non-smoker"
)
ipw_data$smoking_group <- factor(
ipw_data$smoking_group,
levels = c("Non-smoker", "Smoker")
)
ipw_data$gender_label <- factor(
ipw_data$gender,
levels = c("F", "M"),
labels = c("Female", "Male")
)
ipw_data$race_label <- factor(
ipw_data$race,
levels = c(
"mex_american",
"other_hispanic",
"nh_white",
"nh_black",
"nh_asian"
),
labels = c(
"Mexican American",
"Other Hispanic",
"Non-Hispanic White",
"Non-Hispanic Black",
"Non-Hispanic Asian"
)
)
ipw_data$educ_label <- factor(
ipw_data$educ_lvl,
levels = c(
"lt_9th_grade",
"9_to_12th_grade_no_diploma",
"hs_grad_or_ged",
"start_college_to_aa",
"college_grad_or_above"
),
labels = c(
"Less than 9th grade",
"9th-12th, no diploma",
"High school or GED",
"Some college or AA",
"College graduate or above"
)
)
nrow(ipw_data)
## [1] 6299
We use complete cases for this first IPW lesson. That means each person in this analysis has smoking, CVD, and all baseline covariates available.
covariate_dictionary <- data.frame(
Covariate = c(
"Age",
"Gender",
"Race/ethnicity",
"Education",
"Income-to-poverty ratio",
"BMI"
),
Why_it_might_matter = c(
"CVD risk changes strongly with age.",
"Smoking and CVD patterns can differ by gender.",
"Race/ethnicity can reflect social and health differences.",
"Education can be related to health behavior and health care access.",
"Income can be related to living conditions and health care access.",
"BMI can be related to CVD risk."
)
)
knitr::kable(covariate_dictionary)
| Covariate | Why_it_might_matter |
|---|---|
| Age | CVD risk changes strongly with age. |
| Gender | Smoking and CVD patterns can differ by gender. |
| Race/ethnicity | Race/ethnicity can reflect social and health differences. |
| Education | Education can be related to health behavior and health care access. |
| Income-to-poverty ratio | Income can be related to living conditions and health care access. |
| BMI | BMI can be related to CVD risk. |
This is a teaching example. In a real study, choosing covariates requires careful scientific thinking about timing and biology.
Before using IPW, we first look at whether smokers and non-smokers look different on baseline covariates.
mean_sd <- function(x) {
paste0(round(mean(x), 1), " (", round(sd(x), 1), ")")
}
numeric_covariate_table <- data.frame(
Covariate = c("Age", "BMI", "Income-to-poverty ratio"),
Non_smokers = c(
mean_sd(ipw_data$age_yr[ipw_data$smoking_group == "Non-smoker"]),
mean_sd(ipw_data$bmi[ipw_data$smoking_group == "Non-smoker"]),
mean_sd(ipw_data$inc_to_pov_ratio[ipw_data$smoking_group == "Non-smoker"])
),
Smokers = c(
mean_sd(ipw_data$age_yr[ipw_data$smoking_group == "Smoker"]),
mean_sd(ipw_data$bmi[ipw_data$smoking_group == "Smoker"]),
mean_sd(ipw_data$inc_to_pov_ratio[ipw_data$smoking_group == "Smoker"])
)
)
knitr::kable(
numeric_covariate_table,
col.names = c("Covariate", "Non-smokers: mean (SD)", "Smokers: mean (SD)")
)
| Covariate | Non-smokers: mean (SD) | Smokers: mean (SD) |
|---|---|---|
| Age | 46.8 (16.1) | 51.2 (15.6) |
| BMI | 30.1 (7.7) | 30.2 (7.7) |
| Income-to-poverty ratio | 2.8 (1.7) | 2.3 (1.6) |
The boxplots compare the numerical covariates. If the boxes are in different places, the groups differ on that covariate.
format_count_percent <- function(x, group) {
tab <- table(x[group])
pct <- 100 * prop.table(tab)
paste0(as.integer(tab), " (", round(as.numeric(pct), 1), "%)")
}
make_categorical_table <- function(x, label) {
non_smoker <- ipw_data$smoking_group == "Non-smoker"
smoker <- ipw_data$smoking_group == "Smoker"
levels_x <- levels(x)
data.frame(
Covariate = label,
Category = levels_x,
Non_smokers = format_count_percent(x, non_smoker)[levels_x],
Smokers = format_count_percent(x, smoker)[levels_x],
row.names = NULL
)
}
categorical_covariate_table <- rbind(
make_categorical_table(ipw_data$gender_label, "Gender"),
make_categorical_table(ipw_data$race_label, "Race/ethnicity"),
make_categorical_table(ipw_data$educ_label, "Education")
)
knitr::kable(
categorical_covariate_table,
col.names = c("Covariate", "Category", "Non-smokers: n (%)", "Smokers: n (%)")
)
| Covariate | Category | Non-smokers: n (%) | Smokers: n (%) |
|---|---|---|---|
| Gender | Female | NA | NA |
| Gender | Male | NA | NA |
| Race/ethnicity | Mexican American | NA | NA |
| Race/ethnicity | Other Hispanic | NA | NA |
| Race/ethnicity | Non-Hispanic White | NA | NA |
| Race/ethnicity | Non-Hispanic Black | NA | NA |
| Race/ethnicity | Non-Hispanic Asian | NA | NA |
| Education | Less than 9th grade | NA | NA |
| Education | 9th-12th, no diploma | NA | NA |
| Education | High school or GED | NA | NA |
| Education | Some college or AA | NA | NA |
| Education | College graduate or above | NA | NA |
These plots show why we might want adjustment. The two treatment groups are not identical at baseline.
The propensity score is the probability of being treated, based on baseline covariates.
In this lecture:
Propensity score = P(smoker | age, gender, race, education, income, BMI)
We estimate this probability using a logistic regression model.
propensity_model <- glm(
smoker_indicator ~ age_yr + gender_label + race_label +
educ_label + inc_to_pov_ratio + bmi,
data = ipw_data,
family = binomial()
)
ipw_data$propensity_score <- predict(propensity_model, type = "response")
summary(ipw_data$propensity_score)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05416 0.26826 0.39915 0.40625 0.53934 0.85971
The two groups have different propensity score distributions. That means the baseline covariates help predict who was a smoker.
For each person, we use the inverse probability of receiving the treatment status they actually received.
| Person’s actual group | IPW weight |
|---|---|
| Smoker | 1 / propensity score |
| Non-smoker | 1 / (1 - propensity score) |
ipw_data$ipw_weight <- ifelse(
ipw_data$smoker_indicator == 1,
1 / ipw_data$propensity_score,
1 / (1 - ipw_data$propensity_score)
)
weight_summary <- data.frame(
Group = c("Non-smoker", "Smoker", "Overall"),
Minimum = c(
min(ipw_data$ipw_weight[ipw_data$smoking_group == "Non-smoker"]),
min(ipw_data$ipw_weight[ipw_data$smoking_group == "Smoker"]),
min(ipw_data$ipw_weight)
),
Median = c(
median(ipw_data$ipw_weight[ipw_data$smoking_group == "Non-smoker"]),
median(ipw_data$ipw_weight[ipw_data$smoking_group == "Smoker"]),
median(ipw_data$ipw_weight)
),
Mean = c(
mean(ipw_data$ipw_weight[ipw_data$smoking_group == "Non-smoker"]),
mean(ipw_data$ipw_weight[ipw_data$smoking_group == "Smoker"]),
mean(ipw_data$ipw_weight)
),
Maximum = c(
max(ipw_data$ipw_weight[ipw_data$smoking_group == "Non-smoker"]),
max(ipw_data$ipw_weight[ipw_data$smoking_group == "Smoker"]),
max(ipw_data$ipw_weight)
)
)
weight_summary[, -1] <- round(weight_summary[, -1], 2)
knitr::kable(weight_summary)
| Group | Minimum | Median | Mean | Maximum |
|---|---|---|---|---|
| Non-smoker | 1.06 | 1.51 | 1.68 | 5.91 |
| Smoker | 1.16 | 2.05 | 2.46 | 14.65 |
| Overall | 1.06 | 1.70 | 2.00 | 14.65 |
A larger weight means a person represents more people like them in the weighted pseudo-population.
IPW should make the smoker and non-smoker groups more similar in their baseline covariates.
We will use the absolute standardized mean difference, usually called absolute SMD.
Think of SMD as a ruler for covariate imbalance:
SMD close to 0 means the smoker and non-smoker groups look similar on that covariate. Smaller is better.
For a numerical covariate, SMD compares the two group means and divides by a pooled standard deviation. For a categorical covariate, we turn each category into a yes/no variable first. For example, “Male” becomes 1 for male and 0 for not male.
absolute_smd_numeric <- function(x, treatment, weights) {
treated <- treatment == 1
mean_treated <- weighted_mean(x[treated], weights[treated])
mean_control <- weighted_mean(x[!treated], weights[!treated])
var_treated <- weighted_var(x[treated], weights[treated])
var_control <- weighted_var(x[!treated], weights[!treated])
pooled_sd <- sqrt((var_treated + var_control) / 2)
if (pooled_sd == 0) {
0
} else {
abs(mean_treated - mean_control) / pooled_sd
}
}
absolute_smd_binary <- function(x, treatment, weights) {
treated <- treatment == 1
prop_treated <- weighted_mean(as.numeric(x[treated]), weights[treated])
prop_control <- weighted_mean(as.numeric(x[!treated]), weights[!treated])
pooled_sd <- sqrt(
(
prop_treated * (1 - prop_treated) +
prop_control * (1 - prop_control)
) / 2
)
if (pooled_sd == 0) {
0
} else {
abs(prop_treated - prop_control) / pooled_sd
}
}
unweighted_weights <- rep(1, nrow(ipw_data))
make_numeric_smd_row <- function(label, x) {
data.frame(
Covariate = label,
Before_weighting = absolute_smd_numeric(
x,
ipw_data$smoker_indicator,
unweighted_weights
),
After_IPW = absolute_smd_numeric(
x,
ipw_data$smoker_indicator,
ipw_data$ipw_weight
)
)
}
make_factor_smd_rows <- function(label, x, keep_levels = levels(x)) {
rows <- lapply(keep_levels, function(current_level) {
indicator <- x == current_level
data.frame(
Covariate = paste0(label, ": ", current_level),
Before_weighting = absolute_smd_binary(
indicator,
ipw_data$smoker_indicator,
unweighted_weights
),
After_IPW = absolute_smd_binary(
indicator,
ipw_data$smoker_indicator,
ipw_data$ipw_weight
)
)
})
do.call(rbind, rows)
}
smd_table <- rbind(
make_numeric_smd_row("Age", ipw_data$age_yr),
make_numeric_smd_row("BMI", ipw_data$bmi),
make_numeric_smd_row("Income-to-poverty ratio", ipw_data$inc_to_pov_ratio),
make_factor_smd_rows("Gender", ipw_data$gender_label, "Male"),
make_factor_smd_rows("Race/ethnicity", ipw_data$race_label),
make_factor_smd_rows("Education", ipw_data$educ_label)
)
smd_table <- smd_table[order(smd_table$Before_weighting, decreasing = TRUE), ]
smd_table$Before_weighting <- round(smd_table$Before_weighting, 3)
smd_table$After_IPW <- round(smd_table$After_IPW, 3)
summary_smd_table <- data.frame(
Summary = c("Largest absolute SMD", "Average absolute SMD"),
Before_weighting = c(
max(smd_table$Before_weighting),
mean(smd_table$Before_weighting)
),
After_IPW = c(
max(smd_table$After_IPW),
mean(smd_table$After_IPW)
)
)
summary_smd_table$Before_weighting <- round(
summary_smd_table$Before_weighting,
3
)
summary_smd_table$After_IPW <- round(summary_smd_table$After_IPW, 3)
knitr::kable(
smd_table,
col.names = c("Covariate", "Before weighting", "After IPW"),
caption = "Absolute standardized mean differences before and after IPW"
)
| Covariate | Before weighting | After IPW | |
|---|---|---|---|
| 14 | Education: College graduate or above | 0.401 | 0.016 |
| 4 | Gender: Male | 0.373 | 0.015 |
| 9 | Race/ethnicity: Non-Hispanic Asian | 0.337 | 0.012 |
| 7 | Race/ethnicity: Non-Hispanic White | 0.326 | 0.001 |
| 3 | Income-to-poverty ratio | 0.293 | 0.012 |
| 1 | Age | 0.278 | 0.008 |
| 11 | Education: 9th-12th, no diploma | 0.193 | 0.001 |
| 12 | Education: High school or GED | 0.174 | 0.001 |
| 13 | Education: Some college or AA | 0.090 | 0.008 |
| 6 | Race/ethnicity: Other Hispanic | 0.083 | 0.014 |
| 5 | Race/ethnicity: Mexican American | 0.066 | 0.009 |
| 10 | Education: Less than 9th grade | 0.036 | 0.013 |
| 2 | BMI | 0.014 | 0.009 |
| 8 | Race/ethnicity: Non-Hispanic Black | 0.000 | 0.009 |
knitr::kable(
summary_smd_table,
col.names = c("Summary", "Before weighting", "After IPW")
)
| Summary | Before weighting | After IPW |
|---|---|---|
| Largest absolute SMD | 0.401 | 0.016 |
| Average absolute SMD | 0.190 | 0.009 |
After weighting, the absolute SMDs are much smaller. This means the weighted smoker group and weighted non-smoker group are more comparable on the measured baseline covariates.
Now we estimate the CVD risk for each treatment group using the IPW weights.
weighted_risk <- function(outcome, weights) {
sum(weights * outcome) / sum(weights)
}
smoker_rows <- ipw_data$smoking_group == "Smoker"
non_smoker_rows <- ipw_data$smoking_group == "Non-smoker"
observed_risk_smoker <- mean(ipw_data$cvd_indicator[smoker_rows])
observed_risk_non_smoker <- mean(ipw_data$cvd_indicator[non_smoker_rows])
ipw_risk_smoker <- weighted_risk(
ipw_data$cvd_indicator[smoker_rows],
ipw_data$ipw_weight[smoker_rows]
)
ipw_risk_non_smoker <- weighted_risk(
ipw_data$cvd_indicator[non_smoker_rows],
ipw_data$ipw_weight[non_smoker_rows]
)
observed_difference <- observed_risk_smoker - observed_risk_non_smoker
ipw_difference <- ipw_risk_smoker - ipw_risk_non_smoker
observed_ratio <- observed_risk_smoker / observed_risk_non_smoker
ipw_ratio <- ipw_risk_smoker / ipw_risk_non_smoker
result_table <- data.frame(
Method = c("Simple observed comparison", "IPW-adjusted comparison"),
Risk_non_smokers_percent = round(
100 * c(observed_risk_non_smoker, ipw_risk_non_smoker),
1
),
Risk_smokers_percent = round(
100 * c(observed_risk_smoker, ipw_risk_smoker),
1
),
Difference_percentage_points = round(
100 * c(observed_difference, ipw_difference),
1
),
Risk_ratio = round(c(observed_ratio, ipw_ratio), 2)
)
knitr::kable(
result_table,
col.names = c(
"Method",
"Non-smoker risk (%)",
"Smoker risk (%)",
"Difference (percentage points)",
"Risk ratio"
)
)
| Method | Non-smoker risk (%) | Smoker risk (%) | Difference (percentage points) | Risk ratio |
|---|---|---|---|---|
| Simple observed comparison | 6.9 | 15.1 | 8.2 | 2.18 |
| IPW-adjusted comparison | 8.7 | 11.7 | 3.0 | 1.34 |
The simple observed comparison uses the original dataset. The IPW-adjusted comparison uses the weighted pseudo-population.
This plot has two different comparisons.
The blue bars are the simple observed comparison from the original complete-case dataset. In that comparison, CVD risk is 6.9% for non-smokers and 15.1% for smokers. The observed difference is 8.2 percentage points.
The orange bars are the IPW-adjusted comparison. These bars use the weighted pseudo-population, where measured baseline covariates are more balanced between smokers and non-smokers. In the weighted comparison, CVD risk is 8.7% for non-smokers and 11.7% for smokers. The IPW-adjusted difference is 3 percentage points.
The important teaching point is that the gap becomes smaller after IPW. This suggests that part of the large raw smoker versus non-smoker difference was connected to baseline covariate differences between the groups. After balancing the measured covariates, the remaining difference is smaller, but it is still not automatic proof of causation. The IPW interpretation depends on important assumptions, especially that we measured and adjusted for the important confounders.
IPW did not change the actual people in the dataset. It changed how much each person counted.
The goal is to make the comparison closer to the causal question:
What would CVD risk be if everyone were compared under smoking versus no smoking, after accounting for measured baseline covariates?
In this teaching analysis:
We should still be careful. IPW adjusts only for the covariates we measured and included in the model. It does not solve unmeasured confounding, measurement error, or timing problems.
This appendix gives the mathematical notation for the IPW method used in this lecture.
For each person \(i = 1, \ldots, n\), define:
\[ \begin{aligned} T_i &= \text{smoking indicator}, \\ Y_i &= \text{CVD indicator}, \\ X_i &= \text{baseline covariates}. \end{aligned} \]
Here \(T_i = 1\) means smoker, \(T_i = 0\) means non-smoker, and \(Y_i = 1\) means the person has CVD.
Our goal is to estimate two adjusted risks:
\[ \begin{aligned} \psi_1 &= \text{CVD risk if everyone smoked}, \\ \psi_0 &= \text{CVD risk if everyone did not smoke}. \end{aligned} \]
Then the adjusted risk difference and adjusted risk ratio are:
\[ \begin{aligned} RD &= \psi_1 - \psi_0, \\ RR &= \frac{\psi_1}{\psi_0}. \end{aligned} \]
The propensity score is the probability of treatment given baseline covariates:
\[ e(X_i) = P(T_i = 1 \mid X_i). \]
In this lecture, treatment means smoking, so:
\[ e(X_i) = P(\text{smoker} \mid X_i). \]
We estimate this probability using logistic regression:
\[ \text{logit}\{e(X_i)\} = \alpha_0 + \alpha^\top X_i. \]
The fitted propensity score is:
\[ \hat e_i = \hat e(X_i). \]
IPW gives each person a weight based on the inverse probability of the treatment status they actually received.
For smokers:
\[ w_i = \frac{1}{\hat e_i} \quad \text{if } T_i = 1. \]
For non-smokers:
\[ w_i = \frac{1}{1 - \hat e_i} \quad \text{if } T_i = 0. \]
Equivalently:
\[ w_i = \frac{T_i}{\hat e_i} + \frac{1 - T_i}{1 - \hat e_i}. \]
The intuition is:
If a person’s observed smoking status was unlikely for someone with their covariates, that person receives a larger weight.
The lecture uses normalized weighted averages.
The IPW estimate of the smoker risk is:
\[ \hat\psi_1^{IPW} = \frac{ \sum_{i=1}^{n} \dfrac{T_iY_i}{\hat e_i} }{ \sum_{i=1}^{n} \dfrac{T_i}{\hat e_i} }. \]
The IPW estimate of the non-smoker risk is:
\[ \hat\psi_0^{IPW} = \frac{ \sum_{i=1}^{n} \dfrac{(1 - T_i)Y_i}{1 - \hat e_i} }{ \sum_{i=1}^{n} \dfrac{1 - T_i}{1 - \hat e_i} }. \]
These formulas match the R code in the lecture:
\[ \text{weighted mean} = \frac{\sum_i w_iY_i}{\sum_i w_i}. \]
Once we estimate the two weighted risks, the IPW risk difference is:
\[ \widehat{RD}^{IPW} = \hat\psi_1^{IPW} - \hat\psi_0^{IPW}. \]
The IPW risk ratio is:
\[ \widehat{RR}^{IPW} = \frac{\hat\psi_1^{IPW}} {\hat\psi_0^{IPW}}. \]
In the lecture, the risk difference is reported in percentage points:
\[ 100 \times \widehat{RD}^{IPW}. \]
IPW is trying to create a weighted pseudo-population where measured baseline covariates look similar between smokers and non-smokers.
For a numerical covariate \(Z_i\), the weighted mean among smokers is:
\[ \bar Z_1^{w} = \frac{ \sum_{i=1}^{n} T_i w_i Z_i }{ \sum_{i=1}^{n} T_i w_i }. \]
The weighted mean among non-smokers is:
\[ \bar Z_0^{w} = \frac{ \sum_{i=1}^{n} (1 - T_i) w_i Z_i }{ \sum_{i=1}^{n} (1 - T_i) w_i }. \]
The absolute standardized mean difference is:
\[ \text{SMD} = \left| \frac{ \bar Z_1^{w} - \bar Z_0^{w} }{ s_{pooled}^{w} } \right|, \]
where:
\[ s_{pooled}^{w} = \sqrt{ \frac{ (s_1^{w})^2 + (s_0^{w})^2 }{2} }. \]
Here \(s_1^{w}\) and \(s_0^{w}\) are the weighted standard deviations among smokers and non-smokers.
For a categorical covariate, the lecture turns each category into a yes/no indicator and then uses the same idea.
Smaller SMD values mean better balance:
\[ \text{SMD close to } 0 \quad \Longrightarrow \quad \text{better covariate balance}. \]
IPW needs causal assumptions.
The main assumptions are:
If the propensity score model is badly wrong, or if some people have extremely small probabilities of their observed treatment status, IPW can become unstable.
Austin, P. C., & Stuart, E. A. (2015). Moving towards best practice when using inverse probability of treatment weighting using the propensity score to estimate causal treatment effects in observational studies. Statistics in Medicine, 34(28), 3661-3679.
Cole, S. R., & Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology, 168(6), 656-664.
Hernan, M. A., & Robins, J. M. (2020). Causal Inference: What If. Chapman & Hall/CRC.
Rosenbaum, P. R., & Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41-55.