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:

*

*Estimate $E[Y|A,W,V]$

*Estimate $\Pr(A|W, V)$

*Calculate the pseudo-outcomes using the augmented inverse probability weighting formula

*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.
 A: 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.
