Note
-
Download Jupyter notebook:
https://docs.doubleml.org/stable/examples/did/py_panel_simple.ipynb.
Python: Panel Data Introduction#
In this example, we replicate the results from the guide Getting Started with the did Package of the did-R-package.
As the did-R-package the implementation of DoubleML is based on Callaway and Sant’Anna(2021).
The notebook requires the following packages:
[1]:
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression, LogisticRegression
from doubleml.data import DoubleMLPanelData
from doubleml.did import DoubleMLDIDMulti
Data#
The data we will use is simulated and part of the CSDID-Python-Package.
A description of the data generating process can be found at the CSDID-documentation.
[2]:
dta = pd.read_csv("https://raw.githubusercontent.com/d2cml-ai/csdid/main/data/sim_data.csv")
dta.head()
[2]:
| G | X | id | cluster | period | Y | treat | |
|---|---|---|---|---|---|---|---|
| 0 | 3 | -0.876233 | 1 | 5 | 1 | 5.562556 | 1 |
| 1 | 3 | -0.876233 | 1 | 5 | 2 | 4.349213 | 1 |
| 2 | 3 | -0.876233 | 1 | 5 | 3 | 7.134037 | 1 |
| 3 | 3 | -0.876233 | 1 | 5 | 4 | 6.243056 | 1 |
| 4 | 2 | -0.873848 | 2 | 36 | 1 | -3.659387 | 1 |
To work with the DoubleML-package, we initialize a DoubleMLPanelData object.
Therefore, we set the never-treated units in group column G to np.inf (we have to change the datatype to float).
[3]:
# set dtype for G to float
dta["G"] = dta["G"].astype(float)
dta.loc[dta["G"] == 0, "G"] = np.inf
dta.head()
[3]:
| G | X | id | cluster | period | Y | treat | |
|---|---|---|---|---|---|---|---|
| 0 | 3.0 | -0.876233 | 1 | 5 | 1 | 5.562556 | 1 |
| 1 | 3.0 | -0.876233 | 1 | 5 | 2 | 4.349213 | 1 |
| 2 | 3.0 | -0.876233 | 1 | 5 | 3 | 7.134037 | 1 |
| 3 | 3.0 | -0.876233 | 1 | 5 | 4 | 6.243056 | 1 |
| 4 | 2.0 | -0.873848 | 2 | 36 | 1 | -3.659387 | 1 |
Now, we can initialize the DoubleMLPanelData object, specifying
y_col: the outcomed_cols: the group variable indicating the first treated period for each unitid_col: the unique identification column for each unitt_col: the time columnx_cols: the additional pre-treatment controls
[4]:
dml_data = DoubleMLPanelData(
data=dta,
y_col="Y",
d_cols="G",
id_col="id",
t_col="period",
x_cols=["X"]
)
print(dml_data)
================== DoubleMLPanelData Object ==================
------------------ Data summary ------------------
Outcome variable: Y
Treatment variable(s): ['G']
Covariates: ['X']
Instrument variable(s): None
Time variable: period
Id variable: id
No. Unique Ids: 3979
No. Observations: 15916
------------------ DataFrame info ------------------
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 15916 entries, 0 to 15915
Columns: 7 entries, G to treat
dtypes: float64(3), int64(4)
memory usage: 870.5 KB
ATT Estimation#
The DoubleML-package implements estimation of group-time average treatment effect via the DoubleMLDIDMulti class (see model documentation).
The class basically behaves like other DoubleML classes and requires the specification of two learners (for more details on the regression elements, see score documentation). The model will be estimated using the fit() method.
[5]:
dml_obj = DoubleMLDIDMulti(
obj_dml_data=dml_data,
ml_g=LinearRegression(),
ml_m=LogisticRegression(),
control_group="never_treated",
)
dml_obj.fit()
print(dml_obj)
================== DoubleMLDIDMulti Object ==================
------------------ Data summary ------------------
Outcome variable: Y
Treatment variable(s): ['G']
Covariates: ['X']
Instrument variable(s): None
Time variable: period
Id variable: id
No. Unique Ids: 3979
No. Observations: 15916
------------------ Score & algorithm ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
------------------ Machine learner ------------------
Learner ml_g: LinearRegression()
Learner ml_m: LogisticRegression()
Out-of-sample Performance:
Regression:
Learner ml_g0 RMSE: [[1.42736826 1.411621 1.3986617 1.42307038 1.40920164 1.41990816
1.42751955 1.40426596 1.42302084]]
Learner ml_g1 RMSE: [[1.40240005 1.43530525 1.39648269 1.41632322 1.42569556 1.38749145
1.45769169 1.41532362 1.40880946]]
Classification:
Learner ml_m Log Loss: [[0.69062986 0.69155021 0.69141641 0.68011078 0.67936708 0.67914101
0.66241554 0.66239325 0.66247537]]
------------------ Resampling ------------------
No. folds: 5
No. repeated sample splits: 1
------------------ Fit summary ------------------
coef std err t P>|t| 2.5 % 97.5 %
ATT(2.0,1,2) 0.922690 0.064099 14.394704 0.000000 0.797058 1.048322
ATT(2.0,1,3) 1.990212 0.064565 30.825136 0.000000 1.863667 2.116756
ATT(2.0,1,4) 2.960996 0.063272 46.797714 0.000000 2.836985 3.085008
ATT(3.0,1,2) -0.044061 0.065863 -0.668970 0.503514 -0.173150 0.085029
ATT(3.0,2,3) 1.105389 0.065660 16.835148 0.000000 0.976698 1.234079
ATT(3.0,2,4) 2.058165 0.065357 31.491098 0.000000 1.930068 2.186262
ATT(4.0,1,2) 0.002427 0.068537 0.035410 0.971753 -0.131903 0.136757
ATT(4.0,2,3) 0.063886 0.066388 0.962303 0.335897 -0.066233 0.194004
ATT(4.0,3,4) 0.950029 0.067515 14.071342 0.000000 0.817702 1.082356
The summary displays estimates of the \(ATT(g,t_\text{eval})\) effects for different combinations of \((g,t_\text{eval})\) via \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\), where
\(\mathrm{g}\) specifies the group
\(t_\text{pre}\) specifies the corresponding pre-treatment period
\(t_\text{eval}\) specifies the evaluation period
This corresponds to the estimates given in att_gt function in the did-R-package, where the standard choice is \(t_\text{pre} = \min(\mathrm{g}, t_\text{eval}) - 1\) (without anticipation).
Remark that this includes pre-tests effects if \(\mathrm{g} > t_{eval}\), e.g. \(ATT(4,2)\).
As usual for the DoubleML-package, you can obtain joint confidence intervals via bootstrap.
[6]:
level = 0.95
ci = dml_obj.confint(level=level)
dml_obj.bootstrap(n_rep_boot=5000)
ci_joint = dml_obj.confint(level=level, joint=True)
ci_joint
[6]:
| 2.5 % | 97.5 % | |
|---|---|---|
| ATT(2.0,1,2) | 0.748107 | 1.097273 |
| ATT(2.0,1,3) | 1.814361 | 2.166062 |
| ATT(2.0,1,4) | 2.788666 | 3.133327 |
| ATT(3.0,1,2) | -0.223448 | 0.135327 |
| ATT(3.0,2,3) | 0.926556 | 1.284222 |
| ATT(3.0,2,4) | 1.880156 | 2.236174 |
| ATT(4.0,1,2) | -0.184243 | 0.189097 |
| ATT(4.0,2,3) | -0.116932 | 0.244703 |
| ATT(4.0,3,4) | 0.766142 | 1.133916 |
A visualization of the effects can be obtained via the plot_effects() method.
Remark that the plot used joint confidence intervals per default.
[7]:
fig, ax = dml_obj.plot_effects()
/opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
return math.isfinite(val)
/opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
return math.isfinite(val)
Effect Aggregation#
As the did-R-package, the \(ATT\)’s can be aggregated to summarize multiple effects. For details on different aggregations and details on their interpretations see Callaway and Sant’Anna(2021).
The aggregations are implemented via the aggregate() method.
Group Aggregation#
To obtain group-specific effects it is possible to aggregate several \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\) values based on the group \(\mathrm{g}\) by setting the aggregation="group" argument.
[8]:
aggregated = dml_obj.aggregate(aggregation="group")
print(aggregated)
_ = aggregated.plot_effects()
================== DoubleMLDIDAggregation Object ==================
Group Aggregation
------------------ Overall Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
1.487811 0.03422 43.477706 0.0 1.42074 1.554881
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
2.0 1.957966 0.052246 37.476086 0.0 1.855566 2.060366
3.0 1.581777 0.056272 28.109300 0.0 1.471485 1.692069
4.0 0.950029 0.067515 14.071342 0.0 0.817702 1.082356
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
warnings.warn(
The output is a DoubleMLDIDAggregation object which includes an overall aggregation summary based on group size.
Time Aggregation#
This aggregates \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\), based on \(t_\text{eval}\), but weighted with respect to group size. Corresponds to Calendar Time Effects from the did-R-package.
For calendar time effects set aggregation="time".
[9]:
aggregated_time = dml_obj.aggregate("time")
print(aggregated_time)
fig, ax = aggregated_time.plot_effects()
================== DoubleMLDIDAggregation Object ==================
Time Aggregation
------------------ Overall Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
1.481688 0.035099 42.214828 0.0 1.412896 1.55048
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
2 0.922690 0.064099 14.394704 0.0 0.797058 1.048322
3 1.549157 0.051398 30.140522 0.0 1.448419 1.649895
4 1.973217 0.046615 42.329927 0.0 1.881852 2.064581
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
warnings.warn(
Event Study Aggregation#
Finally, aggregation="eventstudy" aggregates \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\) based on exposure time \(e = t_\text{eval} - \mathrm{g}\) (respecting group size).
[10]:
aggregated_eventstudy = dml_obj.aggregate("eventstudy")
print(aggregated_eventstudy)
fig, ax = aggregated_eventstudy.plot_effects()
================== DoubleMLDIDAggregation Object ==================
Event Study Aggregation
------------------ Overall Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
1.992295 0.038734 51.435383 0.0 1.916378 2.068212
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
-2.0 0.002427 0.068537 0.035410 0.971753 -0.131903 0.136757
-1.0 0.011341 0.040481 0.280148 0.779364 -0.068001 0.090682
0.0 0.991804 0.030734 32.270097 0.000000 0.931565 1.052042
1.0 2.024084 0.045602 44.385436 0.000000 1.934705 2.113463
2.0 2.960996 0.063272 46.797714 0.000000 2.836985 3.085008
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
warnings.warn(
Aggregation Details#
The DoubleMLDIDAggregation objects include several DoubleMLFrameworks which support methods like bootstrap() or confint(). Further, the weights can be accessed via the properties
overall_aggregation_weights: weights for the overall aggregationaggregation_weights: weights for the aggregation
To clarify, e.g. for the eventstudy aggregation
[11]:
print(aggregated_eventstudy)
================== DoubleMLDIDAggregation Object ==================
Event Study Aggregation
------------------ Overall Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
1.992295 0.038734 51.435383 0.0 1.916378 2.068212
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
-2.0 0.002427 0.068537 0.035410 0.971753 -0.131903 0.136757
-1.0 0.011341 0.040481 0.280148 0.779364 -0.068001 0.090682
0.0 0.991804 0.030734 32.270097 0.000000 0.931565 1.052042
1.0 2.024084 0.045602 44.385436 0.000000 1.934705 2.113463
2.0 2.960996 0.063272 46.797714 0.000000 2.836985 3.085008
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
Here, the overall effect aggregation aggregates each effect with positive exposure
[12]:
print(aggregated_eventstudy.overall_aggregation_weights)
[0. 0. 0.33333333 0.33333333 0.33333333]
If one would like to consider how the aggregated effect with \(e=0\) is computed, one would have to look at the third set of weights within the aggregation_weights property
[13]:
aggregated_eventstudy.aggregation_weights[2]
[13]:
array([0.32875335, 0. , 0. , 0. , 0.32674263,
0. , 0. , 0. , 0.34450402])
Taking a look at the original dml_obj, one can see that this combines the following estimates:
\(\widehat{ATT}(2,1,2)\)
\(\widehat{ATT}(3,2,3)\)
\(\widehat{ATT}(4,3,4)\)
[14]:
print(dml_obj.summary)
coef std err t P>|t| 2.5 % 97.5 %
ATT(2.0,1,2) 0.922690 0.064099 14.394704 0.000000 0.797058 1.048322
ATT(2.0,1,3) 1.990212 0.064565 30.825136 0.000000 1.863667 2.116756
ATT(2.0,1,4) 2.960996 0.063272 46.797714 0.000000 2.836985 3.085008
ATT(3.0,1,2) -0.044061 0.065863 -0.668970 0.503514 -0.173150 0.085029
ATT(3.0,2,3) 1.105389 0.065660 16.835148 0.000000 0.976698 1.234079
ATT(3.0,2,4) 2.058165 0.065357 31.491098 0.000000 1.930068 2.186262
ATT(4.0,1,2) 0.002427 0.068537 0.035410 0.971753 -0.131903 0.136757
ATT(4.0,2,3) 0.063886 0.066388 0.962303 0.335897 -0.066233 0.194004
ATT(4.0,3,4) 0.950029 0.067515 14.071342 0.000000 0.817702 1.082356