Tutorial 5: Difference-in-Differences

Difference-in-differences (DiD) is a popular design for causal inference. In a 2×2 DiD, we compare the change in outcomes for a treated group to the change for a control group. This tutorial shows how to compute DiD on nonlinear (logit) models using smmargins.

What you will learn

  • How to set up a DiD design with a logit outcome model

  • How to compute cell predictions, simple effects, and the DiD contrast

  • How to compute profile-specific DiD estimates

  • How to interpret the results

The DiD design

In our example, we have:

  • group: treatment group ("A" = control, "B" = treated)

  • preexist_Y: pre-existing condition indicator (0 = without condition, 1 = with condition)

  • condition_X: binary outcome (presence of a health condition)

  • age: age in years

  • female: gender indicator

The pre-existing condition (preexist_Y) serves as our “time” dimension: individuals with the condition represent the “post” period, and those without represent the “pre” period. The interaction between group and preexist_Y captures the DiD effect.

Step 1: Set up the data and model

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from smmargins import Margins

rng = np.random.default_rng(42)
N = 6_000
df_did = pd.DataFrame({
    "group":      rng.choice(["A", "B"], N, p=[0.55, 0.45]),
    "preexist_Y": rng.integers(0, 2, N),
    "age":        rng.normal(55, 15, N).clip(18, 95),
    "female":     rng.integers(0, 2, N),
})
eta = (-3.5 + 0.04*df_did["age"] - 0.3*df_did["female"]
       + 0.5*(df_did["group"]=="B") + 1.1*df_did["preexist_Y"]
       + 0.8*(df_did["group"]=="B")*df_did["preexist_Y"])
df_did["condition_X"] = (rng.uniform(0,1,N) < 1/(1+np.exp(-eta))).astype(int)

fit_did = smf.logit("condition_X ~ C(group) + preexist_Y + C(group):preexist_Y + age + female",
                    data=df_did).fit(disp=False)
M_did = Margins(fit_did)

The model includes an interaction between group and pre-existing condition, which captures the DiD effect. The coefficient on C(group)[T.B]:preexist_Y is the logit-scale interaction.

Step 2: Compute the DiD

Call the did() method with the group variable, the condition (“time”) variable, and their levels:

did = M_did.did("group", "preexist_Y",
                group_levels=["A", "B"],
                condition_levels=[0, 1])

The did object contains three components:

  1. did.cells: Predictions in each of the four cells

  2. did.simple_effects: First differences (simple effects)

  3. did.did: The difference-in-differences contrast

Step 3: View the four cells

did.cells.summary()
prediction std err z P>|z| [95% Conf. Interval]
group=A, preexist_Y=0 0.198022 0.009412 21.038469 2.916422e-98 0.179574 0.216470
group=A, preexist_Y=1 0.410620 0.011475 35.784808 1.903121e-280 0.388130 0.433110
group=B, preexist_Y=0 0.299717 0.011959 25.062174 1.286254e-138 0.276278 0.323156
group=B, preexist_Y=1 0.716206 0.011964 59.865312 0.000000e+00 0.692758 0.739654

These are the predicted probabilities of condition_X in each of the four group-condition combinations. Notice that the outcome is higher for those with preexist_Y = 1 in both groups, and higher in group B overall.

Step 4: View the simple effects

The simple effects show the change from preexist_Y = 0 to preexist_Y = 1 within each group:

did.simple_effects.summary()
simple effect std err z P>|z| [95% Conf. Interval]
group: B vs A | preexist_Y=0 0.101695 0.015219 6.682253 2.352960e-11 0.071867 0.131523
group: B vs A | preexist_Y=1 0.305586 0.016580 18.431265 7.373879e-76 0.273090 0.338082

In group A (control), the predicted probability increases by 25.5 percentage points. In group B (treated), it increases by 33.9 percentage points.

Step 5: View the DiD contrast

The DiD contrast subtracts the simple effect in group A from the simple effect in group B:

did.did.summary()
DiD std err z P>|z| [95% Conf. Interval]
DiD: group(B-A) × preexist_Y(1-0) 0.203891 0.022504 9.06029 1.301042e-19 0.159784 0.247998

The DiD estimate is 0.0835. This means that the increase in the predicted probability of condition_X associated with having preexist_Y = 1 is 8.35 percentage points larger in group B than in group A. This is the causal effect of interest under the parallel trends assumption.

Step 6: Profile-specific DiD

The DiD estimate can vary by covariate profile. We can compute the DiD at specific values of age and gender using atexog:

did_profile = M_did.did("group", "preexist_Y",
                        group_levels=["A", "B"],
                        condition_levels=[0, 1],
                        atexog={"age": 60, "female": 0})

did_profile.cells.summary()
prediction std err z P>|z| [95% Conf. Interval]
group=A, preexist_Y=0 0.241826 0.012700 19.042046 7.648041e-81 0.216935 0.266717
group=A, preexist_Y=1 0.496332 0.015304 32.430587 1.017384e-230 0.466336 0.526328
group=B, preexist_Y=0 0.366280 0.015858 23.097864 4.864630e-118 0.335199 0.397360
group=B, preexist_Y=1 0.801633 0.011731 68.332726 0.000000e+00 0.778640 0.824626 
did_profile.did.summary()
DiD std err z P>|z| [95% Conf. Interval]
DiD: group(B-A) × preexist_Y(1-0) 0.180847 0.02517 7.185068 6.717382e-13 0.131515 0.230179

For a 60-year-old male, the DiD estimate is 0.0938, slightly larger than the population-averaged estimate of 0.0835. This illustrates how the treatment effect can vary across covariate profiles.

Summary: the DiD workflow

  1. Fit a model with the group-condition interaction

  2. Call M.did() with group and condition variables

  3. Inspect cells, simple_effects, and did tables

  4. Use atexog for profile-specific estimates

Recap

In this tutorial we:

  1. Set up a 2×2 DiD design with a logit outcome model

  2. Computed predictions in the four cells

  3. Computed simple effects (first differences) within each group

  4. Computed the DiD contrast (difference of differences)

  5. Showed how to get profile-specific DiD estimates with atexog

Next steps