Causal Analysis Workflow & Estimating ATE Using skmiscpy¶

Here we will show how we can do a very basic causal analysis using the python package skmiscpy. The following contents and the data used in this example are taken from (Heiss, 2024). It is highly recommended to read the post,

Heiss, Andrew. 2024. “Demystifying Causal Inference Estimands: ATE, ATT, and ATU.” March 21, 2024. https://doi.org/10.59350/c9z3a-rcq16.

Does anti-mosquito bed nets reduce the risk of malaria?¶

Here we will use the mosquito_nets_v2.csv data. The data is observational and non-experimental since people choose to use the nets or not based on personal preferences. In this case, income and health confound the net → malaria risk relationship: income influences health, and both income and health influence the choice to use a net and overall malaria risk.

We have four main variables in this data:

• Net use (Treatment): Binary 0/1, TRUE/FALSE variable indicating if the person uses a bed net.

• Malaria risk (Outcome): Scale from 0–100, with higher values representing greater risk.

• Income (Confounder): Weekly income, measured in dollars. Higher income causes a higher probability of using a net; people above the 80th percentile get an extra boost in probability.

• Health (Confounder): Health status, scale from 0–100, with higher values representing better health. Better health causes a higher probability of using a net; people above the 80th percentile get an extra boost in probability.

Setup¶

import pandas as pd
import statsmodels.formula.api as smf
from skmiscpy import compute_smd, plot_smd
from skmiscpy import plot_mirror_histogram

nets = pd.read_csv("../example_data/mosquito_nets_v2.csv")

nets.head()

Since the data used here is non-experimental, we will calculate ATE using IP weighting. For that at first we need to predict the propensity score (i.e. probability of using the bed nets, which is the treatment in this example) and weight the data using this propensity score.

Propensity Score Estimation¶

ps_model = smf.logit(
    "net ~ income + health",
    data = nets
).fit(disp = 0)

data_ps = nets.assign(ps = ps_model.predict(nets))

data_ps[['net', 'malaria_risk', 'ps']].head()

We can draw a mirror histogram to assess the propensity scores. To draw the mirror histogram we will use a function plot_mirror_histogram() from the skmiscpy package.

Visualizing the propensity score¶

plot_mirror_histogram(
    data = data_ps,
    var = 'ps',
    group = 'net',
    title = "Distribution of Propensity Score in the Observed Population",
    xlabel = "Propensity Score"
)

The above graph shows a clear imbalance between the treated group (Group 1:, who chose to use the mosquito net) and the untreated group (Group 0: who didn’t use the mosquito net). Such imbalance in the probability of getting treatment between the treated and untreated group does make sense. Because treatment assignment (which is mosquito nets in this case) was not randomized and assignment is biased by the users’ personal perference towards the mosquito nets. And we are assuming that users’ perference is affected by their health and income.

Now let’s calculate the inverse probability weights (IPW) using which we will create a weighted pseudo population. The formula of IPW for the causal estimand ATE is,

where \(p_i = P(T_i = 1 \vert X)\) is the predicted probability of receiving treatment from a treatment ~ confounder regression model (aka, propensity scores) for each person or unit.

data_wt = (
    data_ps
      .assign(
        wts_ate = data_ps['net'] / data_ps['ps'] + (1 - data_ps['net']) / (1 - data_ps['ps'])
      )
)

data_wt.head()

The plot_mirror_histogram function has an optional argument weights. We will refer the wts_ate column to this argument so that we can see the weighted distribution of propensity score.

plot_mirror_histogram(
    data = data_wt,
    var = 'ps',
    group = 'net',
    weights = 'wts_ate',
    title = "Distribution of Propensity Score in the Weighted Population",
    xlabel = "Propensity Score"
)

These weighted treated and untreated groups look like free of the confounding coming from health and income and are arguably more comparable, acting more like randomized treated and untreated groups.

Checking Covariate Balance¶

Now we will check and compare the standardized mean difference (SMD) for the confounders health & income in the both unweighted and weighted (i.e. adjusted by wts_ate weight) data. We will use compute_smd function from skmiscpy to do this.

smd_df = compute_smd(
    data = data_wt,
    vars = ['health', 'income'],
    group = 'net',
    wt_var = 'wts_ate',
    estimand = 'ATE'
)

smd_df

The adjusted (i.e. weighted) SMDs for both counfounders are less than the which denotes really good covariate balance. We can see the balance using love plot. The plot_smd function from skmiscpy package helps us to draw a love plot.

plot_smd(smd_df, add_ref_line = True, ref_line_value = 0.1)

Estimating the ATE¶

Now we can esitimate the ATE which describes the causal effect of mosquito nets on reducing the malaria risk.

outcome_model = (
    smf.wls(
        formula = 'malaria_risk ~ net', 
        weights = data_wt['wts_ate'], 
        data = data_wt
    )
    .fit()   
)

# Predict outcomes when net = 1
nets_1 = nets.assign(net = 1)
pred_1 = outcome_model.predict(nets_1)

# Predict outcomes when net = 0
nets_0 = nets.assign(net = 0)
pred_0 = outcome_model.predict(nets_0)

mean_difference = (pred_1 - pred_0).mean()
print(mean_difference)

-14.661421260549284

So we can say, using a mosquito net reduces the risk of malaria by 14.7 points, on average across all people in the country. And this is the estimate of ATE.

Now the ATE shows the effect of the net program for everyone, even people who have no need for a net. If we’re interested in making this a universal program, the ATE is useful. But what if we’re interested in what the program is currently doing for the people using it! Then we’d need to find the ATT instead.

ATT (Average Treatment Effect among the Treated)¶

The generic procedure is as same as above. At first we need to calculate the weights, then check the covariate balance from the data using these weights. If satisfying covariate balance is achieved, proceed to ATT estimation.

The formula to calculate the weights in case of ATT is,

data_wt_att = (
    data_ps
      .assign(
        wts_att = data_ps['net'] + (1 - data_ps['net']) * data_ps['ps'] / (1 - data_ps['ps'])
    )
)

data_wt_att.head()

plot_mirror_histogram(
    data = data_wt_att,
    var = 'ps',
    group = 'net',
    weights = 'wts_att',
    title = "Distribution of Propensity Score in the Weighted Population",
    xlabel = "Propensity Score"
)

Checking Covariate Balancing¶

smd_df_att = compute_smd(
    data = data_wt_att,
    vars = ['health', 'income'],
    group = 'net',
    wt_var = 'wts_att',
    estimand = 'ATT'
)

smd_df_att

plot_smd(smd_df_att, add_ref_line=True)

Estimating ATT¶

outcome_model_att = (
    smf.wls(
        "malaria_risk ~ net",
        data = data_wt_att,
        weights = data_wt_att['wts_att']
    )
    .fit()
)


treated = data_wt_att.loc[nets.net == 1]

# Predict outcomes when net = 1 for treated people
pred_1 = (
    outcome_model_att
        .predict(treated.assign(net = 1))
)

# Predict outcomes when net = 0 for treated people
pred_0 = (
    outcome_model_att
        .predict(treated.assign(net = 0))
)

mean_difference = (pred_1 - pred_0).mean()
print(mean_difference)

-17.563307778917196

So that means, People who currently use mosquito nets see a malaria risk reduction of 17.6 points, on average. This ATT estimate is helpful for understanding the effect of net usage on people who actually use them; this shows what would happen if we withheld the program or took away everyone’s nets

Now what if, we’re interested in what would happen if we expanded the program to people not using it? For this, we need to estimate the “Average Treatment Effect among the Control group” (ATC).

ATC (Average Treatment Effect among the Control group)¶

The formula to calculate the weights in case of ATC is,

data_wt_atc = (
    data_ps
      .assign(
          wts_atc = (1 - data_ps['ps']) * data_ps['net'] / data_ps['ps'] + (1 - data_ps['net'])
      )
)

data_wt_atc.head()

plot_mirror_histogram(
    data = data_wt_atc,
    var = 'ps',
    group = 'net',
    weights = 'wts_atc'
)

Checking Covariate Balance¶

smd_df_atc = compute_smd(
    data = data_wt_atc,
    vars = ['income', 'health'],
    group = 'net',
    wt_var = 'wts_atc',
    estimand = 'ATC'
)

smd_df_atc

plot_smd(smd_df_atc, add_ref_line = True)

Estimating ATC¶

outcome_model_atc = (
    smf.wls(
        "malaria_risk ~ net",
        data = data_wt_atc,
        weights = data_wt_atc['wts_atc']
    )
    .fit()
)

untreated = data_wt_atc.loc[nets.net == 0]

pred_1 = (
    outcome_model_atc
      .predict(
          untreated.assign(net = 1)
      )
)

pred_0 = (
    outcome_model_atc
      .predict(
          untreated.assign(net = 0)
      )
)

(pred_1 - pred_0).mean()

np.float64(-11.719789670723681)

So People who don’t currently use mosquito nets would see a malaria risk reduction of 11.7 points, on average. This ATC estimate is helpful for understanding the effect of net usage on people who don’t use them right now; this shows what would happen if we expanded the program or gave nets to people without them