# Variance for a doubly-robust CATE estimator

I am interested in how the variance for the conditional average treatment effect (CATE) is calculated for the doubly robust pseudo-outcome approach. Below are the exact details of the problem and my question

Let $$A$$ be the treatment of interest, $$Y^a$$ be the potential outcome under treatment $$a$$, and $$Y$$ be the observed treatment. There is some set of confounders $$W$$ and we are interested in estimating the CATE by $$V$$, where $$V$$ is a subset of $$W$$. Specifically, the CATE we are interested in is $$E[Y^1-Y^0|V]$$ In the context I am interested in, I am using a parametric model for the above quantity. So something like $$E[Y^1-Y^0|V] = \beta_0 + \beta_1V$$ is sufficient. The doubly-robust pseudo-outcome approach works like the following:

1. Estimate $$E[Y|A,W,V]$$
2. Estimate $$\Pr(A|W, V)$$
3. Calculate the pseudo-outcomes using the augmented inverse probability weighting formula
4. Calculate $$\hat{Y}^1_i - \hat{Y}^0_i$$ from the pseudo-outcomes in step 3, and estimate the $$\beta$$'s in the previous model by fitting some parametric model (OLS in this case)

My question is how is the variance for the CATE / $$\beta$$ calculated?. As far as I am aware, there isn't any reason to use the standard standard error estimates from the regression.

Most papers I have come across use machine learning to estimate step 4, but that is not what I am interested in at this point. Any sources on the variance (or code) would be deeply appreciated.

• If you do everything parametrically, you can fit the whole thing as a system of generalized estimating equations. The variance can then be derived using M-estimation theory.
– Noah
Jul 1, 2020 at 5:29

Returning to answer my own question, one way to estimate the variance is to use M-estimation (or estimating equations as stated by Noah). This method relies on the use of parametric models for $$E[Y|A,W,V]$$ and $$\Pr(A | W,V)$$. Below is a description of how this can be done and an example.

## M-estimation

For an intro to M-estimation, see Stefanski & Boos (2002) or Cole et al. (2022). Essentially, we are going to simultaneously solve a series of estimating equations. Then we use the sandwich variance to estimate the variance for $$\beta$$.

The stacked estimating equations are $$\sum_{i=1}^{n} \psi(O_i; \theta) = \sum_{i=1}^{n} \begin{bmatrix} \left(A_i - \text{expit}(X_i^T \alpha)\right)X_i \\ \left(Y_i - Z_i^T \gamma\right)Z_i \\ \left[(\tilde{Y}_i^1 - \tilde{Y}_i^0) - V_i^T \beta\right]Z_i \\ \end{bmatrix} = 0$$ where $$O_i = (Y_i, A_i, W_i, V_i)$$, $$\theta = (\alpha, \gamma, \beta)$$, $$X_i$$ is the design matrix for the propensity score model, $$Z_i$$ is the design matrix for the outcome model, and $$\tilde{Y}_i^a = \frac{Y_i I(A_i=a)}{A_i \text{expit}(X_i^T \alpha) + (1-A_i)(1-\text{expit}(X_i^T \alpha))} - \frac{{Z_i^a}^T (\text{expit}(X_i^T \alpha) - A_i)}{\text{expit}(X_i^T \alpha)}$$ where $$Z_i^a$$ indicating the design matrix with $$A_i=a$$. Note that $$\tilde{Y}_i^a$$ are the pseudo-potential-outcomes under $$a$$ mentioned in step 3 of the question.

The first estimating equation is the propensity score model (i.e., logistic regression), the second is the outcome model (i.e., linear regression), and the final is the CATE above. To solve the estimating equations, we use a root-finding procedure to find the values of $$\hat{\theta}$$ that result in the sum of the estimating functions being zero (or nearly so).

After, the empirical sandwich variance estimator is used to estimate the full covariance matrix (and the variance for $$\beta$$). The sandwich is $$V(O_i, \hat{\theta}) = B(O_i, \hat{\theta})^{-1} M(O_i, \hat{\theta}) \left\{ B(O_i, \hat{\theta})^{-1} \right\}^T$$ where $$B(O_i, \hat{\theta}) = n^{-1} \sum_i -\psi'(O_i, \hat{\theta})$$ with $$\psi'$$ indicating the matrix derivative and $$M(O_i, \hat{\theta}) = n^{-1} \sum_i \psi(O_i, \hat{\theta}) \psi(O_i, \hat{\theta})^T$$ This variance estimator allows the uncertainty of $$\hat{\beta}$$ to depend on the uncertainty of $$\hat{\alpha}$$ and $$\hat{\gamma}$$.

## Code

To show how this would be implemented programmatically, the following is a simple example using the Python M-estimation library, delicatessen (disclosure: I am the creator of delicatessen).

Initializing

import numpy as np
import pandas as pd
from delicatessen import MEstimator
from delicatessen.estimating_equations import ee_regression
from delicatessen.utilities import inverse_logit


Generating some generic data

np.random.seed(20230109)
n = 1000
d = pd.DataFrame()
d['W'] = np.random.normal(size=n)
d['V'] = np.random.binomial(n=1, p=0.5, size=n)
pr_a = logistic.cdf(-2 + 2*d['V'] + d['W'])
d['A'] = np.random.binomial(n=1, p=pr_a, size=n)
d['Y1'] = 3 + 1*1 + 1*1*d['V'] + d['V'] - 1.2*d['W'] + np.random.normal(size=n)
d['Y0'] = 3 + 0*1 + 0*1*d['V'] + d['V'] - 1.2*d['W'] + np.random.normal(size=n)
d['Y'] = np.where(d['A'] == 1, d['Y1'], d['Y0'])
d['intercept'] = 1
d['AV'] = d['A']*d['V']


Setting up design matrices as NumPy arrays for delicatessen

X = np.asarray(d[['intercept', 'V', 'W']])
a = np.asarray(d['A'])
Z = np.asarray(d[['intercept', 'A', 'AV', 'V', 'W']])
y = np.asarray(d['Y'])
V = np.asarray(d[['intercept', 'V']])

da = d.copy()
da['A'] = 1
da['AV'] = da['A']*d['V']
Z1 = np.asarray(da[['intercept', 'A', 'AV', 'V', 'W']])
da['A'] = 0
da['AV'] = da['A']*d['V']
Z0 = np.asarray(da[['intercept', 'A', 'AV', 'V', 'W']])


Defining stacked estimating equations

def psi(theta):
# Separating theta into the constituent components
alpha = theta[0:3]    # Prop score params
gamma = theta[3:8]    # Outcome params
beta = theta[8:]      # CATE params

# Estimating equations for propensity score model
ee_ps = ee_regression(theta=alpha,        # Model parameters
X=X,                # ... X design matrix
y=a,                # ... treatment
model='logistic')   # ... with logistic model

# Estimating equations for outcome model
ee_out = ee_regression(theta=gamma,       # Model parameters
X=Z,               # ... Z design matrix
y=y,               # ... observed outcome
model='linear')    # ... with linear model

# Generating pseudo potential outcomes
pi_a = inverse_logit(np.dot(X, alpha))             # Propensity score
y_a1 = np.dot(Z1, gamma)                           # Predicted Y^1
y_a0 = np.dot(Z0, gamma)                           # Predicted Y^0
yt_a1 = y*a/pi_a - y_a1*(a-pi_a)/pi_a              # Pseudo Y^1
yt_a0 = y*(1-a)/(1-pi_a) + y_a0*(a-pi_a)/(1-pi_a)  # Pseudo Y^0

# Estimating CATE model
ee_cate = ee_regression(theta=beta,        # Parameters
X=V,               # ... with CATE design matrix
y=yt_a1 - yt_a0,   # ... difference is pseudo outcomes
model='linear')    # ... with linear model

# Returning stack of estimating equations
return np.vstack([ee_ps, ee_out, ee_cate])


Applying the M-estimator

estr = MEstimator(psi, init=[0, ]*10)
estr.estimate(solver='lm')

estr.theta[8:]         # Point estimates
estr.variance[8:, 8:]  # Covariance matrix


which results in $$\hat{\beta}_0 = 1.00$$ (variance: 0.020) and $$\hat{\beta}_1 = 0.85$$ (variance: 0.029).

More examples of M-estimators and an introduction to delicatessen is provided in Zivich et al. (2022).

## Caveat

This approach may not work when machine learning is used to estimate the nuisance models. Instead, the influence function could be used to derive a variance estimator.

## References

Cole, S. R., Edwards, J. K., Breskin, A., Rosin, S., Zivich, P. N., Shook-Sa, B. E., & Hudgens, M. G. (2022). Illustration of Two Fusion Designs and Estimators. American Journal of Epidemiology.

Stefanski, L. A., & Boos, D. D. (2002). The calculus of M-estimation. The American Statistician, 56(1), 29-38.

Zivich, P. N., Klose, M., Cole, S. R., Edwards, J. K., & Shook-Sa, B. E. (2022). Delicatessen: M-Estimation in Python. arXiv preprint arXiv:2203.11300.