EPLM

Robins-Newey partially linear E-estimator

1 Where it fits

Group: Causal inference

EPLM targets a scalar treatment effect in a partially linear model. It combines an outcome equation with a working model for \(E[D\mid W]\) and solves the resulting stacked moment system.

The intended estimand is the coefficient on the scalar treatment after accounting for controls \(W\) without treating the nuisance regression as the object of interest.

2 Estimating equations

Let \(\tilde W=[\mathbf1,W]\). The class first fits the linear treatment projection

\[ \hat\pi=\arg\min_\pi\|D-\tilde W\pi\|_2^2, \qquad \hat v=D-\tilde W\hat\pi, \]

and then estimates

\[ \hat\beta = \frac{\hat v'Y}{\hat v'D}. \]

Because OLS makes \(\hat v\) orthogonal to \(\tilde W\), the denominator equals \(\hat v'\hat v\) up to numerical error. This is the partialling-out coefficient in the partially linear specification

\[ Y=\beta D+g(W)+U \]

when both the treatment projection and control adjustment are restricted to the supplied linear control basis.

Inference treats \(\theta=(\pi',\beta)'\) as the exactly identified solution to the stacked per-observation moments

\[ \psi_i(\theta) = \begin{bmatrix} \tilde W_i(D_i-\tilde W_i'\pi)\\ (D_i-\tilde W_i'\pi)(Y_i-D_i\beta) \end{bmatrix}. \]

The treatment is not required to be binary; it may be continuous.

3 Implementation walkthrough

The point estimate and stacked-moment inference are both native dense linear algebra.

1. fit() validates aligned finite arrays and prepends an explicit intercept to \(W\). It solves the treatment projection by rectangular least squares, then stores \(v=D-\tilde W\hat\pi\).
2. The scalar coefficient is computed directly as \(v'Y/(v'D)\) rather than by a second regression object. A denominator below \(10^{-12}\) in absolute value is treated as lack of residual treatment variation and raises.
3. The stored parameter vector is ordered $(‘,)’`. At summary time, the code reconstructs every observation’s full moment row by scaling \(\tilde W_i\) with \(v_i\) and appending \(v_i(Y_i-D_i\beta)\).
4. For each parameter coordinate it perturbs by \(h_j=\text{fd eps}\max(|\theta_j|,1)\), reevaluates the complete moment matrix at plus and minus perturbations, averages columns, and fills one central-difference Jacobian column. No closed-form Jacobian is used even though this particular system is analytic.
5. The square Jacobian is inverted once. Moment rows are mapped to parameter-score rows as \(\psi_iA^{-T}/n\); the shared iid, HC1, Newey-West, or cluster aggregator then supplies the covariance. Only the final scalar row and column are exposed publicly.

The implementation deliberately treats nuisance estimation and target inference as one exactly identified system. This accounts for the linear treatment projection algebraically, at the cost of repeated moment construction during numerical differentiation.

4 Inference

The mean-moment Jacobian \(A=\partial\bar\psi/\partial\theta'\) is evaluated by central differences. The uncorrected iid option, named vanilla, is already an empirical sandwich:

\[ \widehat V = \frac1n\hat A^{-1} \left\{\frac1n\sum_i\psi_i\psi_i'\right\} \hat A^{-T}. \]

HC1 multiplies the corresponding parameter-score covariance by \(n/(n-k)\), where \(k\) includes nuisance parameters. Newey-West applies Bartlett lag weights and that same correction; cluster covariance sums parameter scores by cluster and applies the package’s \(G/(G-1)\times(n-1)/(n-k)\) correction. The summary returns only the \(\beta\) variance block and nuisance point estimates. There is no bootstrap or built-in Wald method.

5 Performance and numerical behavior

Point estimation is a pair of dense least-squares calculations with approximately \(O(np^2+p^3)\) work for \(p\) control columns. Summary-time numerical differentiation evaluates the full moment system twice per element of \((\pi,\beta)\) and then inverts its dense Jacobian, so inference can cost more than fitting. Near-zero residual treatment variation raises. There is no sample splitting, nonlinear nuisance learner, or regularization, and collinear controls can make both projection and covariance unstable.

6 Python API

Constructor: cm.EPLM

Call fit(y, d, w) with scalar treatment d and 2D controls w. summary(vcov=None, lags=None, clusters=None) returns the coefficient, standard error, covariance matrix, and nuisance coefficients. There is no predict() method.

print(inspect.signature(cm.EPLM))

(fd_eps=1e-06)

cls = cm.EPLM
display(HTML(html_table(["Public method"], public_methods(cls))))

Public method
EPLM(fd_eps=1e-06)
fit(self, /, y, d, w)
summary(self, /, vcov=None, lags=None, clusters=None)

7 Minimal example

rng = np.random.default_rng(11)
w = rng.normal(size=(300, 3))
d = 0.4 + w @ np.array([0.6, -0.2, 0.3]) + rng.normal(size=300)
y = 1.1 * d + w @ np.array([0.2, 0.1, -0.2]) + rng.normal(scale=0.5, size=300)
model = cm.EPLM()
model.fit(y, d, w)
print(model.summary()['coef'])
print(model.summary()['se'])

1.0900606444971932
0.027112508384756203

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

rng = np.random.default_rng(111)
w = rng.normal(size=(120, 3))
d = 0.4 + w @ np.array([0.6, -0.2, 0.3]) + rng.normal(size=120)
y = 1.1 * d + w @ np.array([0.2, 0.1, -0.2]) + rng.normal(size=120) * 0.5
model = cm.EPLM()
model.fit(y, d, w)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))

summary() key	shape
coef	()
se	()
vcov	(1, 1)
nuisance_coef	(4,)