4. Score functions#
We use method-of-moments estimators for the target parameter \(\theta_0\) based upon the empirical analog of the moment condition
where we call \(\psi\) the score function, \(W=(Y,D,X,Z)\), \(\theta_0\) is the parameter of interest and \(\eta\) denotes nuisance functions with population value \(\eta_0\). We use score functions \(\psi(W; \theta, \eta)\) that satisfy \(\mathbb{E}[ \psi(W; \theta_0, \eta_0)] = 0\) with \(\theta_0\) being the unique solution and that obey the Neyman orthogonality condition
An integral component for the object-oriented (OOP) implementation of
DoubleMLPLR
,
DoubleMLPLIV
,
DoubleMLIRM
,
and DoubleMLIIVM
is the linearity of the score function in the parameter \(\theta\)
Hence the estimator can be written as
The linearity of the score function in the parameter \(\theta\) allows the implementation of key components in a very
general way.
The methods and algorithms to estimate the causal parameters, to estimate their standard errors, to perform a multiplier
bootstrap, to obtain confidence intervals and many more are implemented in the abstract base class DoubleML
.
The object-oriented architecture therefore allows for easy extension to new model classes for double machine learning.
This is doable with very minor effort whenever the linearity of the score function is satisfied.
4.1. Implementation of the score function and the estimate of the causal parameter#
As an example we consider a partially linear regression model (PLR)
implemented in DoubleMLPLR
.
In [1]: import doubleml as dml
In [2]: from doubleml.datasets import make_plr_CCDDHNR2018
In [3]: from sklearn.ensemble import RandomForestRegressor
In [4]: from sklearn.base import clone
In [5]: np.random.seed(3141)
In [6]: learner = RandomForestRegressor(n_estimators=100, max_features=20, max_depth=5, min_samples_leaf=2)
In [7]: ml_l = clone(learner)
In [8]: ml_m = clone(learner)
In [9]: data = make_plr_CCDDHNR2018(alpha=0.5, return_type='DataFrame')
In [10]: obj_dml_data = dml.DoubleMLData(data, 'y', 'd')
In [11]: dml_plr_obj = dml.DoubleMLPLR(obj_dml_data, ml_l, ml_m)
In [12]: dml_plr_obj.fit();
In [13]: print(dml_plr_obj)
================== DoubleMLPLR Object ==================
------------------ Data summary ------------------
Outcome variable: y
Treatment variable(s): ['d']
Covariates: ['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10', 'X11', 'X12', 'X13', 'X14', 'X15', 'X16', 'X17', 'X18', 'X19', 'X20']
Instrument variable(s): None
No. Observations: 500
------------------ Score & algorithm ------------------
Score function: partialling out
DML algorithm: dml2
------------------ Machine learner ------------------
Learner ml_l: RandomForestRegressor(max_depth=5, max_features=20, min_samples_leaf=2)
Learner ml_m: RandomForestRegressor(max_depth=5, max_features=20, min_samples_leaf=2)
------------------ Resampling ------------------
No. folds: 5
No. repeated sample splits: 1
Apply cross-fitting: True
------------------ Fit summary ------------------
coef std err t P>|t| 2.5 % 97.5 %
d 0.462834 0.041047 11.275605 1.731842e-29 0.382383 0.543285
library(DoubleML)
library(mlr3)
library(mlr3learners)
library(data.table)
lgr::get_logger("mlr3")$set_threshold("warn")
learner = lrn("regr.ranger", num.trees = 100, mtry = 20, min.node.size = 2, max.depth = 5)
ml_l = learner$clone()
ml_m = learner$clone()
set.seed(3141)
data = make_plr_CCDDHNR2018(alpha=0.5, return_type='data.table')
obj_dml_data = DoubleMLData$new(data, y_col="y", d_cols="d")
dml_plr_obj = DoubleMLPLR$new(obj_dml_data, ml_l, ml_m)
dml_plr_obj$fit()
print(dml_plr_obj)
================= DoubleMLPLR Object ==================
------------------ Data summary ------------------
Outcome variable: y
Treatment variable(s): d
Covariates: X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20
Instrument(s):
No. Observations: 500
------------------ Score & algorithm ------------------
Score function: partialling out
DML algorithm: dml2
------------------ Machine learner ------------------
ml_l: regr.ranger
ml_m: regr.ranger
------------------ Resampling ------------------
No. folds: 5
No. repeated sample splits: 1
Apply cross-fitting: TRUE
------------------ Fit summary ------------------
Estimates and significance testing of the effect of target variables
Estimate. Std. Error t value Pr(>|t|)
d 0.54440 0.04512 12.06 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The fit()
method of DoubleMLPLR
stores the estimate \(\tilde{\theta}_0\) in its coef
attribute.
In [14]: print(dml_plr_obj.coef)
[0.46283393]
print(dml_plr_obj$coef)
d
0.5443965
The values of the score function components \(\psi_a(W_i; \hat{\eta}_0)\) and \(\psi_b(W_i; \hat{\eta}_0)\)
are stored in the attributes psi_a
and psi_b
.
In the attribute psi
the values of the score function \(\psi(W_i; \tilde{\theta}_0, \hat{\eta}_0)\) are stored.
In [15]: print(dml_plr_obj.psi[:5])
[[[-0.15351459]]
[[ 0.02031318]]
[[ 0.00576115]]
[[ 0.08091435]]
[[-0.23145105]]]
print(dml_plr_obj$psi[1:5, ,1])
[1] -0.0009695237 0.7811465543 0.0090193584 -0.4037269089 0.8646627426
4.2. Implemented Neyman orthogonal score functions#
4.2.1. Partially linear regression model (PLR)#
For the PLR model implemented in DoubleMLPLR
one can choose between
score='partialling out'
and score='IV-type'
.
score='partialling out'
implements the score function:
with \(\eta=(\ell,m)\) and where the components of the linear score are
score='IV-type'
implements the score function:
with \(\eta=(g,m)\) and where the components of the linear score are
4.2.2. Partially linear IV regression model (PLIV)#
For the PLIV model implemented in DoubleMLPLIV
one can choose between
score='IV-type'
and score='partialling out'
.
score='partialling out'
implements the score function:
with \(\eta=(\ell, m, r)\) and where the components of the linear score are
score='IV-type'
implements the score function:
with \(\eta=(g,m)\) and where the components of the linear score are
4.2.3. Interactive regression model (IRM)#
For the IRM model implemented in DoubleMLIRM
one can choose between
score='ATE'
and score='ATTE'
.
score='ATE'
implements the score function:
with \(\eta=(g,m)\) and where the components of the linear score are
score='ATTE'
implements the score function:
with \(\eta=(g, m, p)\) and where the components of the linear score are
4.2.4. Interactive IV model (IIVM)#
For the IIVM model implemented in DoubleMLIIVM
we employ for score='LATE'
the score function:
score='LATE'
implements the score function:
with \(\eta=(g, m, r)\) and where the components of the linear score are
4.3. Specifying alternative score functions via callables#
Via callables user-written score functions can be used.
This functionality is at the moment only implemented for specific model classes in Python.
For the PLR model implemented in DoubleMLPLR
an alternative score function can be
set via score
.
Choose a callable object / function with signature score(y, d, g_hat, m_hat, smpls)
which returns
the two score components \(\psi_a()\) and \(\psi_b()\).
For example, the non-orthogonal score function
can be obtained with
In [16]: import numpy as np
In [17]: def non_orth_score(y, d, l_hat, m_hat, g_hat, smpls):
....: u_hat = y - g_hat
....: psi_a = -np.multiply(d, d)
....: psi_b = np.multiply(d, u_hat)
....: return psi_a, psi_b
....:
non_orth_score = function(y, d, l_hat, m_hat, g_hat, smpls) {
u_hat = y - g_hat
psi_a = -1*d*d
psi_b = d*u_hat
psis = list(psi_a = psi_a, psi_b = psi_b)
return(psis)
}
Use DoubleMLPLR
with inf_model=non_orth_score
in order to obtain the estimator
when applying fit()
.
Note that this estimate will in general be prone to a regularization bias, see also Overcoming regularization bias by orthogonalization.