1 Overview

The causal question in the smoking-CVD lectures is:

\[ \text{Is the risk of cardiovascular disease different if everyone smoked versus if everyone did not smoke?} \]

The observed data do not come from a randomized smoking experiment. Smokers and non-smokers may differ in age, gender, race/ethnicity, education, income-to-poverty ratio, BMI, and other baseline risk factors. The propensity score is one way to organize this problem.

Conceptual Goal

The propensity score is not an outcome model and not a causal assumption. It is a scalar summary of the measured baseline covariates that predicts treatment assignment: \[ e(X)=P(T=1\mid X). \] In this lecture, \(T=1\) means lifetime smoker and \(T=0\) means non-smoker. The propensity score helps us compare smokers and non-smokers who had similar measured chances of being smokers.

2 Definition

For a binary treatment \(T\in\{0,1\}\), the propensity score is

\[ e(X)=P(T=1\mid X). \]

In the smoking example,

\[ e(X) = P(\text{smoker}=1\mid \text{age},\text{gender},\text{race},\text{education},\text{income},\text{BMI}). \]

The propensity score is a probability. A participant with \(\hat e(X_i)=0.70\) is not “70% smoker” in any biological sense. It means that, among participants with similar measured covariates, the fitted treatment model expects about 70% to be smokers.

3 Why a Single Score Can Help

The remarkable property of the propensity score is that it is a balancing score. For the true score \(e(X)\),

\[ T \perp X \mid e(X). \]

This means that, within a narrow band of the propensity score, the distribution of measured covariates should be similar between treated and untreated participants. The statement is about measured \(X\), not about unmeasured confounders.

Balancing Score Result

Because \(e(X)=P(T=1\mid X)\), conditioning on \(X\) and on \(e(X)\) gives \[ P(T=1\mid X,e(X)) = e(X). \] Also, \[ P(T=1\mid e(X)) = e(X). \] Therefore, once \(e(X)\) is fixed, knowing the full measured covariate vector \(X\) does not further change the treatment probability. That is the sense in which \(T\) is independent of measured \(X\) within levels of the propensity score.

Read the picture from left to right. A high-dimensional covariate vector \(X\) is compressed into one fitted score \(\hat e(X)\). A propensity-score stratum is a narrow interval on that score axis. The balancing claim says that, if the score is correctly specified and the interval is sufficiently narrow, observed smokers and observed non-smokers inside the same band should have similar measured covariate distributions.

3.1 A More Formal Balancing Argument

The balancing score statement can be written more rigorously. Let \(A\) be any set of measured covariate values and suppose we condition on the event \(e(X)=r\), where \(0<r<1\). Then

\[ \begin{aligned} P(X\in A\mid T=1,e(X)=r) &= \frac{P(T=1\mid X\in A,e(X)=r)P(X\in A\mid e(X)=r)} {P(T=1\mid e(X)=r)} \\ &= \frac{r\,P(X\in A\mid e(X)=r)}{r} \\ &= P(X\in A\mid e(X)=r). \end{aligned} \]

The same calculation for \(T=0\) gives

\[ P(X\in A\mid T=0,e(X)=r) = P(X\in A\mid e(X)=r). \]

Thus, within a true propensity-score stratum, the distribution of measured \(X\) is the same in the treated and untreated groups:

\[ X\perp T\mid e(X). \]

This result explains why a one-dimensional score can be useful for design. It does not say that unmeasured variables are balanced; only variables included in \(X\), or variables whose imbalance is captured by \(X\), are addressed.

4 Propensity Scores and Causal Identification

The usual observational-study identification condition is conditional exchangeability:

\[ \{Y(1),Y(0)\}\perp T\mid X. \]

Under positivity, the true propensity score can be used as a lower-dimensional conditioning variable:

\[ \{Y(1),Y(0)\}\perp T\mid e(X). \]

The intuition is that if \(X\) blocks the backdoor paths between smoking and CVD, and \(e(X)\) balances the measured \(X\)’s, then conditioning on the true score can also remove measured confounding. In practice we estimate \(e(X)\), so diagnostics focus on whether the estimated score creates adequate balance and overlap.

5 Load the Classroom NHANES Data

The tutorial uses the same complete-case classroom analytic dataset as the IPW lecture. The lectures deliberately do not use NHANES survey weights here, so the numerical values are for the classroom analytic sample.

data_candidates <- c(
  file.path("..", "Data", "NHANES_data.csv"),
  file.path("Data", "NHANES_data.csv"),
  "C:/Users/seyoo/OneDrive/GitHub_repository/Causal_Inference_Survival_Analysis/Causal_Inference_Lecture/Data/NHANES_data.csv"
)

data_path <- data_candidates[file.exists(data_candidates)][1]
if (is.na(data_path)) {
  stop("Could not find NHANES_data.csv. Check the data path.")
}

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)

ps_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]
ps_data$smoker_indicator <- as.integer(ps_data$smoker_indicator)
ps_data$cvd_indicator <- as.integer(ps_data$cvd_indicator)

ps_data$smoking_group <- factor(
  ifelse(ps_data$smoker_indicator == 1, "Smoker", "Non-smoker"),
  levels = c("Non-smoker", "Smoker")
)

ps_data$gender_label <- factor(
  ps_data$gender,
  levels = c("F", "M"),
  labels = c("Female", "Male")
)

ps_data$race_label <- factor(
  ps_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"
  )
)

ps_data$educ_label <- factor(
  ps_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",
    "9-12th grade, no diploma",
    "High school/GED",
    "Some college/AA",
    "College graduate or above"
  )
)

nrow(ps_data)
## [1] 6299
sample_table <- data.frame(
  Group = c("Non-smoker", "Smoker", "Total"),
  N = c(
    sum(ps_data$smoker_indicator == 0),
    sum(ps_data$smoker_indicator == 1),
    nrow(ps_data)
  ),
  CVD_percent = round(100 * c(
    mean(ps_data$cvd_indicator[ps_data$smoker_indicator == 0]),
    mean(ps_data$cvd_indicator[ps_data$smoker_indicator == 1]),
    mean(ps_data$cvd_indicator)
  ), 2)
)

knitr::kable(sample_table)
Group N CVD_percent
Non-smoker 3740 6.93
Smoker 2559 15.12
Total 6299 10.26

6 Baseline Imbalance Before Using the Propensity Score

Before estimating a propensity score, we check whether baseline covariates differ between smokers and non-smokers. If they do, then the crude CVD comparison mixes treatment differences with covariate differences.

baseline_summary <- data.frame(
  Covariate = c("Age", "BMI", "Income-to-poverty ratio", "Male (%)"),
  Non_smoker = c(
    round(mean(ps_data$age_yr[ps_data$smoker_indicator == 0]), 2),
    round(mean(ps_data$bmi[ps_data$smoker_indicator == 0]), 2),
    round(mean(ps_data$inc_to_pov_ratio[ps_data$smoker_indicator == 0]), 2),
    round(100 * mean(ps_data$gender_label[ps_data$smoker_indicator == 0] == "Male"), 2)
  ),
  Smoker = c(
    round(mean(ps_data$age_yr[ps_data$smoker_indicator == 1]), 2),
    round(mean(ps_data$bmi[ps_data$smoker_indicator == 1]), 2),
    round(mean(ps_data$inc_to_pov_ratio[ps_data$smoker_indicator == 1]), 2),
    round(100 * mean(ps_data$gender_label[ps_data$smoker_indicator == 1] == "Male"), 2)
  )
)

knitr::kable(baseline_summary)
Covariate Non_smoker Smoker
Age 46.84 51.24
BMI 30.14 30.25
Income-to-poverty ratio 2.81 2.34
Male (%) 40.40 58.73

For categorical covariates, the analogous object is not a density curve but a vector of category probabilities. The following panels compare the observed category proportions among non-smokers and smokers before weighting. Within each category, longer blue or orange bars mean that category is more common in that observed group.

These pictures are not about the outcome. They show that the treatment groups can have different baseline covariate distributions before adjustment. Age, BMI, gender, race/ethnicity, and education are only marginal views; imbalance can also occur in combinations of these covariates. The propensity score is designed to help repair this measured imbalance.

7 Estimate the Propensity Score

The lecture uses a simple logistic treatment model:

\[ \operatorname{logit}\{e(X)\} = \alpha_0+\alpha_1\text{age}+\alpha_2\text{gender} +\alpha_3\text{race}+\alpha_4\text{education} +\alpha_5\text{income}+\alpha_6\text{BMI}. \]

propensity_model <- glm(
  smoker_indicator ~ age_yr + gender_label + race_label +
    educ_label + inc_to_pov_ratio + bmi,
  data = ps_data,
  family = binomial()
)

ps_data$e_hat <- predict(propensity_model, type = "response")

summary(ps_data$e_hat)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05416 0.26826 0.39915 0.40625 0.53934 0.85971
ps_coef <- summary(propensity_model)$coefficients
main_ps_terms <- c("age_yr", "bmi", "inc_to_pov_ratio", "gender_labelMale")
ps_coef_table <- data.frame(
  Term = main_ps_terms,
  Log_odds_coefficient = round(ps_coef[main_ps_terms, "Estimate"], 3),
  Odds_ratio_for_smoking = round(exp(ps_coef[main_ps_terms, "Estimate"]), 3),
  P_value = signif(ps_coef[main_ps_terms, "Pr(>|z|)"], 3)
)

knitr::kable(ps_coef_table)
Term Log_odds_coefficient Odds_ratio_for_smoking P_value
age_yr age_yr 0.018 1.019 0.000
bmi bmi -0.006 0.994 0.126
inc_to_pov_ratio inc_to_pov_ratio -0.169 0.845 0.000
gender_labelMale gender_labelMale 0.795 2.214 0.000

The fitted value \(\hat e(X_i)\) is each participant’s estimated probability of being a smoker given the measured covariates. It is not the probability of CVD.

8 Reading the Propensity Score Distribution

The plot should be read horizontally. Around a given score value, we ask whether both observed smokers and observed non-smokers exist. If both groups appear in a region of the score, then comparisons in that region are supported by data. If one group is absent, the analysis would have to extrapolate.

9 Overlap in the NHANES Covariate Space

The propensity score is one-dimensional, but the original overlap problem lives in the multivariate covariate space. The slide-style plot below shows one two-dimensional slice of that space for our classroom NHANES data: age and BMI. Red x’s are observed smokers and black open circles are observed non-smokers.

This plot is only a two-covariate projection, so it cannot diagnose full multivariate positivity by itself. Its purpose is conceptual: overlap means that, for a relevant covariate region, both observed treatment states are represented. If a region contains almost only smokers or almost only non-smokers, comparisons there rely heavily on modeling or weighting extrapolation.

ps_data$age_quartile <- cut(
  ps_data$age_yr,
  breaks = quantile(ps_data$age_yr, probs = seq(0, 1, by = 0.25)),
  include.lowest = TRUE
)

ps_data$bmi_quartile <- cut(
  ps_data$bmi,
  breaks = quantile(ps_data$bmi, probs = seq(0, 1, by = 0.25)),
  include.lowest = TRUE
)

age_bmi_cells <- aggregate(
  smoker_indicator ~ age_quartile + bmi_quartile,
  data = ps_data,
  FUN = function(z) c(n = length(z), smoker_percent = 100 * mean(z))
)

age_bmi_cells <- data.frame(
  age_quartile = age_bmi_cells$age_quartile,
  bmi_quartile = age_bmi_cells$bmi_quartile,
  n = age_bmi_cells$smoker_indicator[, "n"],
  smoker_percent = age_bmi_cells$smoker_indicator[, "smoker_percent"]
)

age_bmi_cells$non_smoker_percent <- 100 - age_bmi_cells$smoker_percent
age_bmi_cells$minimum_group_percent <- pmin(
  age_bmi_cells$smoker_percent,
  age_bmi_cells$non_smoker_percent
)

low_overlap_cells <- age_bmi_cells[
  order(age_bmi_cells$minimum_group_percent),
  c("age_quartile", "bmi_quartile", "n", "smoker_percent",
    "non_smoker_percent", "minimum_group_percent")
]

knitr::kable(
  head(
    transform(
      low_overlap_cells,
      smoker_percent = round(smoker_percent, 1),
      non_smoker_percent = round(non_smoker_percent, 1),
      minimum_group_percent = round(minimum_group_percent, 1)
    ),
    8
  ),
  caption = "Age-BMI strata with the weakest empirical overlap"
)
Age-BMI strata with the weakest empirical overlap
age_quartile bmi_quartile n smoker_percent non_smoker_percent minimum_group_percent
5 [20,35] (24.9,28.9] 344 29.1 70.9 29.1
1 [20,35] [14.6,24.9] 562 32.2 67.8 32.2
9 [20,35] (28.9,34] 346 32.4 67.6 32.4
13 [20,35] (34,92.3] 392 36.5 63.5 36.5
6 (35,49] (24.9,28.9] 385 37.4 62.6 37.4
7 (49,62] (24.9,28.9] 438 38.6 61.4 38.6
2 (35,49] [14.6,24.9] 343 38.8 61.2 38.8
10 (35,49] (28.9,34] 372 39.0 61.0 39.0

The last column is the smaller of the smoker and non-smoker percentages inside an age-BMI cell. A small value means one treatment state is rare in that cell. This is a finite-sample diagnostic for practical overlap; it is not a formal proof or disproof of positivity.

10 Positivity and Extreme Scores

The positivity condition for the ATE target is

\[ 0<P(T=1\mid X=x)<1 \]

for covariate values \(x\) in the target population. In propensity-score notation,

\[ 0<e(x)<1. \]

This condition has a practical version: the estimated scores should not be too close to 0 or 1 for participants who matter to the target population.

The inverse probability weights are large exactly when a participant received a treatment that the model regarded as unlikely for their covariates. Such participants are influential because they represent many covariate-similar people who did not receive that treatment.

11 From Propensity Score to IPW

For the ATE-style marginal risks used in the IPW lecture, the population identities are

\[ E\left\{\frac{TY}{e(X)}\right\}=E\{Y(1)\}, \qquad E\left\{\frac{(1-T)Y}{1-e(X)}\right\}=E\{Y(0)\}. \]

The weights are inverse probabilities of the treatment that was actually observed:

\[ w_i^{ATE} = \frac{T_i}{e(X_i)} + \frac{1-T_i}{1-e(X_i)}. \]

ps_data$ate_weight <- ifelse(
  ps_data$smoker_indicator == 1,
  1 / ps_data$e_hat,
  1 / (1 - ps_data$e_hat)
)

weight_summary <- data.frame(
  Group = c("Non-smoker", "Smoker", "Overall"),
  Mean = round(c(
    mean(ps_data$ate_weight[ps_data$smoker_indicator == 0]),
    mean(ps_data$ate_weight[ps_data$smoker_indicator == 1]),
    mean(ps_data$ate_weight)
  ), 2),
  Median = round(c(
    median(ps_data$ate_weight[ps_data$smoker_indicator == 0]),
    median(ps_data$ate_weight[ps_data$smoker_indicator == 1]),
    median(ps_data$ate_weight)
  ), 2),
  Maximum = round(c(
    max(ps_data$ate_weight[ps_data$smoker_indicator == 0]),
    max(ps_data$ate_weight[ps_data$smoker_indicator == 1]),
    max(ps_data$ate_weight)
  ), 2)
)

knitr::kable(weight_summary)
Group Mean Median Maximum
Non-smoker 1.68 1.51 5.91
Smoker 2.46 2.05 14.65
Overall 2.00 1.70 14.65

The next displays mirror the pseudo-population picture in the slides, but use the classroom NHANES data. Because the full \(X\) is multivariate, we first look at one-dimensional slices of it: age and BMI. The same weighting idea applies jointly to all measured covariates in the propensity-score model.

In each display, the top row shows the observed treatment-specific distribution, the middle gray panel shows the analytic-sample target distribution, and the bottom row shows the weighted treatment-specific distribution after applying the ATE inverse-probability weights. If the propensity-score model is adequate, the weighted smoker and weighted non-smoker distributions should move toward the same target distribution:

\[ \frac{\sum_{i=1}^n T_i\hat e(X_i)^{-1} I(X_i\in A)} {\sum_{i=1}^n T_i\hat e(X_i)^{-1}} \approx \frac{1}{n}\sum_{i=1}^n I(X_i\in A), \]

and

\[ \frac{\sum_{i=1}^n (1-T_i)\{1-\hat e(X_i)\}^{-1} I(X_i\in A)} {\sum_{i=1}^n (1-T_i)\{1-\hat e(X_i)\}^{-1}} \approx \frac{1}{n}\sum_{i=1}^n I(X_i\in A). \]

For categorical covariates, take \(A=\{C=c\}\), where \(C\) is a category-valued baseline variable. The plot below shows the same target-distribution logic for gender, race/ethnicity, and education. The gray bars are the full analytic-sample target proportions; the blue and orange bars are the IPW-weighted non-smoker and smoker proportions.

Thus, the weighted pseudo-population is not a new dataset with new people. It is a reweighted empirical distribution in which covariate patterns that were rare within a treatment group receive more influence.

12 Diagnosing Balance

The main practical question after estimating a propensity score is not whether the treatment model has a small p-value. The question is whether the score creates balance in measured baseline covariates.

ps_data$male_indicator <- as.integer(ps_data$gender_label == "Male")
ps_data$white_indicator <- as.integer(ps_data$race_label == "Non-Hispanic White")
ps_data$college_indicator <- as.integer(ps_data$educ_label == "College graduate or above")

balance_table <- data.frame(
  Covariate = c(
    "Age",
    "BMI",
    "Income-to-poverty ratio",
    "Male",
    "Non-Hispanic White",
    "College graduate or above"
  ),
  Before = c(
    absolute_smd_numeric(ps_data$age_yr, ps_data$smoker_indicator),
    absolute_smd_numeric(ps_data$bmi, ps_data$smoker_indicator),
    absolute_smd_numeric(ps_data$inc_to_pov_ratio, ps_data$smoker_indicator),
    absolute_smd_binary(ps_data$male_indicator, ps_data$smoker_indicator),
    absolute_smd_binary(ps_data$white_indicator, ps_data$smoker_indicator),
    absolute_smd_binary(ps_data$college_indicator, ps_data$smoker_indicator)
  ),
  After_IPW = c(
    absolute_smd_numeric(ps_data$age_yr, ps_data$smoker_indicator, ps_data$ate_weight),
    absolute_smd_numeric(ps_data$bmi, ps_data$smoker_indicator, ps_data$ate_weight),
    absolute_smd_numeric(ps_data$inc_to_pov_ratio, ps_data$smoker_indicator, ps_data$ate_weight),
    absolute_smd_binary(ps_data$male_indicator, ps_data$smoker_indicator, ps_data$ate_weight),
    absolute_smd_binary(ps_data$white_indicator, ps_data$smoker_indicator, ps_data$ate_weight),
    absolute_smd_binary(ps_data$college_indicator, ps_data$smoker_indicator, ps_data$ate_weight)
  )
)

balance_table$Before <- round(balance_table$Before, 3)
balance_table$After_IPW <- round(balance_table$After_IPW, 3)

knitr::kable(balance_table)
Covariate Before After_IPW
Age 0.278 0.008
BMI 0.014 0.009
Income-to-poverty ratio 0.293 0.012
Male 0.373 0.015
Non-Hispanic White 0.326 0.001
College graduate or above 0.401 0.016

This balance plot is one of the most important propensity score diagnostics. Good balance does not prove the causal assumptions, but poor balance shows that the estimated score has not yet made the observed treatment groups comparable on measured covariates.

13 Propensity Score Stratification

Another way to understand the propensity score is through stratification. Split participants into strata with similar \(\hat e(X)\), then compare smokers and non-smokers within each stratum.

ps_data$ps_quintile <- cut(
  ps_data$e_hat,
  breaks = quantile(ps_data$e_hat, probs = seq(0, 1, by = 0.2)),
  include.lowest = TRUE
)

strata_table <- aggregate(
  cbind(age_yr, bmi, inc_to_pov_ratio) ~ ps_quintile + smoking_group,
  data = ps_data,
  FUN = mean
)

strata_counts <- as.data.frame(table(ps_data$ps_quintile, ps_data$smoking_group))
names(strata_counts) <- c("ps_quintile", "smoking_group", "N")
strata_table <- merge(strata_table, strata_counts,
                      by = c("ps_quintile", "smoking_group"))

knitr::kable(
  data.frame(
    PS_quintile = strata_table$ps_quintile,
    Group = strata_table$smoking_group,
    N = strata_table$N,
    Mean_age = round(strata_table$age_yr, 1),
    Mean_BMI = round(strata_table$bmi, 1),
    Mean_income_to_poverty = round(strata_table$inc_to_pov_ratio, 2)
  )
)
PS_quintile Group N Mean_age Mean_BMI Mean_income_to_poverty
(0.241,0.344] Non-smoker 900 44.2 30.6 2.62
(0.241,0.344] Smoker 360 46.0 31.7 2.82
(0.344,0.45] Non-smoker 749 47.1 31.1 2.54
(0.344,0.45] Smoker 510 46.4 31.3 2.48
(0.45,0.573] Non-smoker 625 52.4 30.5 2.32
(0.45,0.573] Smoker 635 51.8 30.2 2.16
(0.573,0.86] Non-smoker 418 57.6 30.7 2.02
(0.573,0.86] Smoker 842 58.0 29.0 1.80
[0.0542,0.241] Non-smoker 1048 41.4 28.6 3.78
[0.0542,0.241] Smoker 212 43.3 30.0 3.88

Each panel reads vertically within a quintile. The blue and orange points are not meant to be identical across different quintiles; higher propensity-score quintiles can legitimately contain older or higher-risk participants. The diagnostic question is whether, inside the same quintile, smokers and non-smokers have similar covariate means. The short gray vertical segments mark the remaining within-quintile differences.

Stratification makes the balancing idea tangible: instead of comparing all smokers with all non-smokers, we compare people with similar estimated chances of smoking. The price is that balance is checked stratum by stratum and may still be imperfect, especially if a quintile contains few participants from one treatment group.

14 What the Propensity Score Is Not

Common Misinterpretations

  1. The propensity score is not an outcome risk. It predicts smoking, not CVD.
  2. A high propensity score does not mean the participant should smoke. It means their measured covariates resemble those of observed smokers.
  3. A good treatment model does not by itself prove no unmeasured confounding.
  4. A propensity score model should use pre-treatment covariates, not variables affected by treatment.
  5. Extreme scores are not harmless. They can create unstable weights and reveal poor overlap.

15 Summary

The propensity score is the treatment probability given measured baseline covariates:

\[ e(X)=P(T=1\mid X). \]

For this lecture, the main ideas are:

  1. The propensity score compresses many measured covariates into one scalar treatment probability.
  2. The true propensity score is a balancing score: within levels of \(e(X)\), measured covariates are balanced between treated and untreated groups.
  3. The propensity score is useful only under the causal assumption that the measured covariates are sufficient for exchangeability.
  4. Positivity requires overlap: both smokers and non-smokers must be possible for relevant covariate patterns.
  5. IPW uses inverse propensity-score weights to create a weighted pseudo-population in which measured covariates are more comparable.
  6. Balance and overlap diagnostics are central. Coefficients and p-values from the treatment model are secondary.

The conceptual roadmap is:

\[ \text{measured baseline }X \longrightarrow \text{propensity score }e(X) \longrightarrow \text{balance / overlap diagnostics} \longrightarrow \text{weighted or matched causal comparison.} \]

16 References

Rosenbaum, P. R., and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41-55.

Rosenbaum, P. R., and Rubin, D. B. (1984). Reducing bias in observational studies using subclassification on the propensity score. Journal of the American Statistical Association, 79(387), 516-524.

Hernan, M. A., and Robins, J. M. (2020). Causal Inference: What If. Boca Raton: Chapman & Hall/CRC.

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.

Stuart, E. A. (2010). Matching methods for causal inference: A review and a look forward. Statistical Science, 25(1), 1-21.

Cole, S. R., and Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology, 168(6), 656-664.