Standard errors

The package also computes “oracle” standard errors, in addition to the usual standard errors reported above. These can be accessed in the data frame est_f (or est_o for the overlap sample results)

print(m1$est_f, digits = 3)
#>                PL     OWN     ATE      EW      CW
#> Black     -0.2574 -0.2482 -0.2655 -0.2550 -0.2604
#> pop_se     0.0281  0.0291  0.0298  0.0289  0.0292
#> oracle_se      NA      NA  0.0298  0.0288  0.0290
#> Hispanic  -0.2931 -0.2829 -0.2992 -0.2862 -0.2944
#> pop_se     0.0260  0.0267  0.0299  0.0268  0.0279
#> oracle_se      NA      NA  0.0299  0.0268  0.0278
#> Asian     -0.2621 -0.2609 -0.2599 -0.2611 -0.2694
#> pop_se     0.0343  0.0343  0.0418  0.0343  0.0475
#> oracle_se      NA      NA  0.0418  0.0344  0.0465
#> Other     -0.1563 -0.1448 -0.1503 -0.1447 -0.1522
#> pop_se     0.0369  0.0370  0.0359  0.0368  0.0370
#> oracle_se      NA      NA  0.0359  0.0360  0.0366

These oracle standard errors (oracle_se) don’t account for estimation error in the propensity score, in contrast to the default standard errors (pop_se), see Methods section below.

Specifying the cluster argument allows for computation of clustered standard errors:

## cluster in interviewer ID
m1alt <- multe(r1, "race", cluster = factor(factor(fl$interviewer_ID_24)))
print(m1alt, digits = 3)
#> Estimates on full sample:
#>               PL     OWN     ATE      EW      CW
#> Black    -0.2574 -0.2482 -0.2655 -0.2550 -0.2604
#> SE        0.0412  0.0425  0.0410  0.0422  0.0420
#> Hispanic -0.2931 -0.2829 -0.2992 -0.2862 -0.2944
#> SE        0.0441  0.0454  0.0495  0.0457  0.0474
#> Asian    -0.2621 -0.2609 -0.2599 -0.2611 -0.2694
#> SE        0.0521  0.0523  0.0619  0.0523  0.0675
#> Other    -0.1563 -0.1448 -0.1503 -0.1447 -0.1522
#> SE        0.0404  0.0416  0.0410  0.0414  0.0428
#> 
#> P-values for null hypothesis of no propensity score variation:
#> Wald test: 0, LM test: 0 
#> 
#> SD(estimated propensity score), maximum over treatment arms:
#> Full sample: 0.196

Methods

This section describes the implementation of the bias decomposition formula and the implementation of alternative estimators. Relative to Goldsmith-Pinkham et al. (2024), we generalize the setup to allow for sampling weights \(\omega_{i}^{2}\) (setting the sampling weights to one recovers unweighted formulas). We also explicitly deal with overlap issues.

We are interested in the effect of treatment \(D_{i}\in\{0,1,\dotsc, K\}\) on an outcome \(Y_{i}\). Let \(X_{i}=\1{D_{i}=1,\dotsc, D_{i}=K}\) denote a vector of treatment indicators, let \(X_{i0}=\1{D_{i}=0}\), and let \(Z_{i}=(1,W_{i}')'\) denote a vector of controls, including an intercept. We focus on the case where the controls enter linearly, so that control functions take the form \(\mathcal{G}=\{z'\gamma\colon \gamma\in\mathbb{R}^{1+\dim(W_i)}\}\).

We assume that \(\mu_{k}(W_i):=E[Y_{i}(k)\mid W_{i}]\in\mathcal{G}\), so that we may write \(\mu_{k}(W_{i})=W_{i}'\alpha_{k}\) for some vectors \(\alpha_{k}\), \(k=0,\dotsc, K\). The average treatment effect (ATE) conditional on \(W_{i}\) is then given by \(\tau(W_{i})=Z_{i}'(\alpha_{k}-\alpha_{0})\), and \(\alpha_{k}\) correspond to the coefficients in the interacted regression \[\begin{equation} Y_{i}=\sum_{k=0}^{K}X_{ik}Z_{i}'\alpha_{k}+\dot{U}_{i}\label{eq:linear_mu}, % =Z_{i}'\alpha_{0}+\sum_{k=1}^{K}X_{ik}Z_{i}'\gamma_{k}\label{eq:linear_gamma} % &=\alpha_{0, k}+W_{i}'\alpha_{W, 0}+\sum_{k=1}^{K}X_{ik}\tau_{k}+\sum_{k=1}^{K}X_{ik}(W_{i}-E[W_{i}])'\gamma_{W, k}. % \label{eq:linear_tau} \end{equation}\] with \(\dot{U}_{i}\) conditionally mean zero. The uninteracted partially linear (PL) estimator is given by estimating \[\begin{equation}\label{eq:uninteracted_linear} Y_i=\sum_{k=1}^{K} X_{ik}\beta + Z_i'\phi +U_{i}, \end{equation}\] by weighted least squares (WLS), yielding \(\hat{\beta}=(\sum_{i}\omega_{i}^{2}\dot{X}_{i}\dot{X}_{i}')^{-1}\sum_{i}\omega_{i}^{2}\dot{X}_{i}Y_{i}\), where \(\dot{X}\) is the sample residual from WLS regression of \(X_{i}\) onto \(Z_{i}\). By Proposition 1 in Goldsmith-Pinkham et al. (2024), the population analog of \(\hat{\beta}\), \(\beta\), has the decomposition \[\begin{equation*} \beta=E[\diag(\Lambda_{i})\tau(W_{i})]+E[\Lambda_{i}-\diag(\Lambda_{i})\tau(W_{i})], \end{equation*}\] where \(\Lambda_{i}=E[\tilde{X}_{i}\tilde{X}_{i}']^{-1}E[\tilde{X}_{i}X_{i}\mid W_{i}]\), and \(\tilde{X}_{i}\) is the population analog of \(\dot{X}_{i}\), the population residual from regressing \(X_{i}\) onto \(W_{i}\). Let \(\hat{\alpha}_{k}\) denote the WLS estimates based on \(\eqref{eq:linear_mu}\). By construction, the sample residuals from estimating \(\eqref{eq:linear_mu}\) and \(Z_{i}\) are both orthogonal to \(\dot{X}_{i}\). As a result, we obtain the exact decomposition \[\begin{equation}\label{eq:1} \begin{split} \hat{\beta} &=E_{n}[\dot{X}_{i}\dot{X}_{i}']^{-1}E_{n}[\dot{X}_{i}Y_{i}]= % E_{n}[\dot{X}_{i}\dot{X}_{i}']^{-1} % E_{n}[\dot{X}_{i}\sum_{k=1}^{K}X_{ik}Z_{i}'(\hat{\alpha}_{k}-\hat{\alpha}_{0})] \\ & = % E_{n}[\hat{\Lambda}_{i}\hat{\Gamma}'Z_{i}]= E_{n}[\diag(\hat{\Lambda}_{i})\hat{\Gamma}'Z_{i}]+E_{n}[(\hat{\Lambda}_{i}-\diag(\hat{\Lambda}_{i}))\hat{\Gamma}'Z_{i}] =:\hbo+\hbcb, \end{split} \end{equation}\] where \(\hat{\Gamma}\) is a matrix with columns \(\hat{\gamma}_{k}=\hat{\alpha}_{k}-\hat{\alpha}_{0}\), \(\hat{\Lambda}_{i}=E_{n}[\dot{X}_{i}\dot{X}_{i}']^{-1}\dot{X}_{i}X_{i}'\), and \(E_{n}[A_{i}]=\sum_{i}\omega_{i}^{2}A_{i}/\sum_{i}\omega_{i}^{2}\) denotes the weighted sample mean.

To compute this decomposition, we don’t need to explicitly compute \(\hat{\Lambda}_{i}\). Instead, we use the identity \[\begin{equation*} \hbo_{k}=e_{k}'E_{n}[\dot{X}_{i}\dot{X}_{i}']^{-1}E_{n}[\dot{X}_{i}X_{ik}Z_{i}'\hat{\gamma}_{k}]=\hat{\delta}_{k k}' \hat{\gamma}_{k}, \end{equation*}\] where \(\hat{\delta}_{k k}\) is a WLS estimator of the system of regressions \[\begin{equation}\label{eq:delta_system} Z_{i}X_{ik}=\delta_{k k}X_{ik}+\sum_{\ell\neq k}\delta_{k \ell}X_{i\ell} +\Delta_{Z, k}Z_{i}+\zeta_{ik}. \end{equation}\]

Note this decomposition and associated standard errors, in the next subsection, are purely regression-based, so they remain valid even if \(X_{i}\) is not a set of binary indicators. Likewise, misspecification of the interacted regression only affects the interpretation of the decomposition; if \(\mu_k\) is not linear, the decomposition will not consistently estimate the contamination bias.

The own treatment weights in this decomposition sum to one, and the contamination weights sum to zero, since \(E_{n}[\hat{\Lambda}_{i}]=I_{k}\), mimicking the population decomposition. If the propensity score doesn’t satisfy \(p_{k}\in\mathcal{G}\), the implied estimate of \(\Lambda(w)\), \[\begin{equation*} \hat{\Lambda}(w)=\frac{1}{n}\sum_{i=1}^{n}\1{W_{i}=w}\hat{\Lambda}_{i} \end{equation*}\] may not be positive definite; in particular, the diagonal elements may not all be positive, in line with Proposition 1 in Goldsmith-Pinkham et al. (2024).

In addition to this decomposition, the package also computes the following alternative estimators:

1. The unweighted ATE estimator, \(\hbate_{k}=E_{n}[W_{i}]'\hat{\gamma}_{k}\)
2. The one-treatment-at-a time estimator \(\hbew_{k}\) that fits \(\eqref{eq:uninteracted_linear}\) using only observations with \(D_{i}\in\{0,k\}\). In other words, it estimates the regression \[\begin{equation}\label{eq:1-at-a-time} Y_{i}=\ddot{\phi}_{k}+X_{ik}\ddot{\beta}_{k}+W_{i}^\prime\ddot{\phi}_{k}+ \ddot{U}_{ik}, \end{equation}\] among observations with \(D_{i}\in\{0,k\}\). This estimator weights the treatment effects using the variance-minimizing weighting scheme given in Corollary 1 in Goldsmith-Pinkham et al. (2024). Consequently, we refer to as the efficiently weighted ATE estimator.
3. The common weights estimator \(\hbcw\), given by the WLS regression of \(Y_{i}\) onto \(X_{i}\), weighting each observation by \[\begin{equation*} \frac{\omega_{i}^{2}\pi_{D_{i}}(1-\pi_{D_{i}})}{\hat{p}_{D_{i}}(W_{i})\sum_{k=0}^{K}\hat{p}_{k}(W_{i})^{-1}}, \end{equation*}\] where, by default, the probabilities \(\pi_{k}\) correspond to the marginal probability \(E_{n}[X_{ik}]\) in the dataset. The propensity scores \(\hat{p}_{k}(W_{i})\) are based on fitting a multinomial logit model for the treatments. This estimator estimates a weighted ATE with weights \(\lambda^{\text{CW}}(W_{i})=\left(\sum_{k=0}^{K}\frac{\pi_{k}(1-\pi_{k})}{p_{k}(W_{i})}\right)^{-1}\). By Corollary 2 in Goldsmith-Pinkham et al. (2024), this weighting scheme minimizes the average variance, under homoskedasticity, across all treatment comparisons—comparisons of outcomes under treatment \(k\) vs treatment \(\ell\), if we draw the treatments \(k\) and \(\ell\) independently from the marginal treatment distribution \((\pi_{0},\dotsc,\pi_{K})\). Option cw_uniform=TRUE in the multe function sets these probabilities to \(1/K\); setting the option to its default, FALSE, sets them to \((E_{n}[X_{i0}],\dotsc,E_{n}[X_{i K}])\).

Derivations

We first derive \(\eqref{eq:psi_cw}\). Observe first that the common weights estimator is identical to the two-step GMM estimator that in the first step, fits a multinomial logit model \[\begin{equation}\label{eq:m-n-l} P(D_{i}=k\mid W_{i})=\frac{e^{Z_{i}'\theta_{k}}}{\sum_{k'=0}^{K}e^{Z_{i}'\theta_{k'}}}=: \pi_{k}(W_{i},\theta), \end{equation}\] with the normalization \(\theta_{0}=0\). In the second step, we use the moment condition \[\begin{equation}\label{eq:moment_overlap_weights} E\left[\frac{\lcw(W_{i};\theta)X_{ik}}{\pi_{k}(W_{i};\theta)}(Y_{i}-\alpha^{\text{CW}}_{k})\right]=0. \end{equation}\] to obtain estimates \(\hat\alpha^{\text{CW}}_{k}\), and set \(\hbcw_k=\hat\alpha^{\text{CW}}_{k}-\hat\alpha^{\text{CW}}_{0}\).

Let \[\begin{equation*} \zeta_{k}(W_{i};\hat{\theta})=\frac{\lcw(W_{i};\hat{\theta})}{\pi_{k}(W_{i},\hat{\theta})} % =\frac{1}{\pi_{k}(W_{i},\hat{\theta})\sum_{k'=0}^{K}\frac{\pi_{k'}(1-\pi_{k'})}{\pi_{k}(W_{i};\hat{\theta})}} =\frac{e^{-Z_{i}'\hat{\theta}_{k}}}{\sum_{j=0}^{K}\pi_{j}(1-\pi_{j})e^{-Z_{i}'\hat{\theta}_{j}}} \end{equation*}\] By equation (6.6) in Newey and McFadden (1994), the influence function of this two-step estimator is given by \[\begin{equation*} \psi_{i}(\hat{\alpha}^{\textnormal{CW}}_{k})= \frac{1}{E[\lcw(W_{i})]}\left( \frac{\lcw(W_{i})X_{ik}}{\pi_{k}(W_{i})}(Y_{i}-\alpha^{\text{CW}}_{k}) +M_{k}(\theta)\psi_{i}(\theta)\right), \end{equation*}\] where \(\psi_{i}(\theta)\) is the influence function of the multinomial logit estimator \(\hat{\theta}\), and \(M_{k}(\theta)\) is the derivative of \(\eqref{eq:moment_overlap_weights}\) wrt \(\theta\).

Since \(\partial \zeta_{k}(W_{i};\theta)/\partial \theta_{j}= Z_{i}\zeta_{k}(W_{i};\theta)[\pi_{j}(1-\pi_{j})\zeta_{k}(W_{i};\theta)-\1{k=j>0}]\), it follows that \[\begin{equation*} M_{k}(\theta)= E[(\eta-e_{k})\otimes Z_{i}\cdot \zeta_{k}X_{ik}(Y_{i}-\alpha_{\lambda^{C},k})],\qquad\eta=( \pi_{j}(1-\pi_{j})\zeta_{j}, \dotsc,\pi_{K}(1-\pi_{K})\zeta_{K})'. \end{equation*}\] Since the multinomial logit log-likelihood is given by \(\ell_{i}=\sum_{k=0}^{K}X_{k}\log(\pi_{k})=\sum_{k=0}^{K}X_{k}Z_{i}'\theta_{k} -\log(\sum_{k=0}^{K}e^{Z_{i}'\theta_{k}})\), the score and the Hessian are \[\begin{equation*} S_{i}(\theta)=(X_{i}-\pi(W_{i};\theta))\otimes Z_{i},\qquad H(\theta)=-E_{n}[(\diag(\pi(W_{i}))-\pi(W_{i})\pi(W_{i})')\otimes Z_{i}Z_{i}'], \end{equation*}\] Since \(\psi_{i}(\theta)=-H(\theta)^{-1}S_{i}(\theta)\), this yields \[\begin{equation}\label{eq:a_i} a_{i}(\theta)=\hat{M}_{k}(\hat{\theta}) \hat{H}(\hat{\theta})^{-1}S_{i}(\theta), \qquad S_{i}(\theta)= (X_{i}-\pi(W_{i};\theta))\otimes Z_{i}, \end{equation}\] with \(\hat{M}_{k}({\theta})=\left(\sum_{i}(\zeta-e_{k})\otimes Z_{i}\cdot \zeta_{k}X_{ik}(Y_{i}-\alpha_{\lambda^{C},k})\right)\), and \(\hat{H}({\theta})= \sum_{i}(\diag(\pi(W_{i}))-\pi(W_{i})\pi(W_{i})')\otimes Z_{i}Z_{i}'\). When \(\pi_k=1/(K+1)\) this formula reduces to that in Theorem 1 in Li and Li (2019).

Next, we show \(\eqref{eq:hbew_oracle}\). Note that it follows from Proposition 2 in Goldsmith-Pinkham et al. (2024) that the efficient influence function is given by \[\begin{multline*} \psi_{i}({\alpha}_{\lambda_{k},k0})= \frac{\lambda(W_{i})}{E[\lambda(W_{i})]}\left(\frac{X_{ik}}{p_{k}(W_{i})}(Y_{i}-\mu_{k}(W_{i})) - \frac{X_{i0}}{p_{0}(W_{i})}(Y_{i}-\mu_{0}(W_{i}))\right)\\ = \frac{X_{ik}+X_{i0}}{E[\lambda(W_{i})]}\left(X_{ik}-r_{k}\right)V_{i}= \frac{(X_{ik}+X_{i0})(X_{ik}-r_{k})V_{i}}{E[(X_{ik}+X_{i0})(X_{ik}-r_{k})^{2}]}, \end{multline*}\] where \(r_{k}=r_{k}(W_{i})=E[X_{ik}\mid W_i,X_{ik}+X_{i0}=1]\). The result then follows since \(\hat{X}_{ik}\) is an estimator of \(X_{ik}-r_{k}(W_{i})\).

Finally, we show \(\eqref{eq:overlap_oracle}\). The derivative of the moment condition \(\eqref{eq:moment_overlap_weights}\) with respect to \(\pi_{k}=p_{k}\) (assuming correct specification of the propensity score) is given by \[\begin{equation*} -E[\lambda\frac{X_{ik}}{p_{k}^{2}(W_{i})} (\mu_{k}-\alpha^{\text{CW}}_{k})\dot{p}_{k}(W_{i})], \end{equation*}\] where we write \(\lambda\) for \(\lcw(W_{i})\). Since \(p_{k}\) is a projection, by Proposition 4 in Newey (1994), the influence function for \(\hat{\alpha}^{\text{CW}}_{k}\) is given by \[\begin{multline*} \frac{1}{E[\lambda]}\left( \lambda \frac{X_{ik}}{p_{k}(W_{i})}(Y_{i}-\alpha^{\text{CW}}_{k}) -\frac{\lambda}{p_{k}(W_{i})}(\mu_{k}(W_{i})-\alpha^{\text{CW}}_{k})(X_{ik}-p_{k}(W_{i}))\right)\\ = \frac{1}{E[\lambda]}\left( \lambda\frac{X_{ik}}{p_{k}(W_{i})}(Y_{i}-\mu_{k}(W_{i})) + \lambda(\mu_{k}(W_{i})-\alpha^{\text{CW}}_{k}) \right). \end{multline*}\] Next, as noted in Abadie et al. (2014), we can view \(\tilde{\alpha}^{\text{CW}}_{k}=\sum_{i}\lambda\mu_{k}(W_{i})/\sum_{i}\lambda\) as an estimator of \({\alpha}^{\text{CW}}_{k}\) based on the moment condition \(E[\lambda(\mu_{k}(W_{i})-{\alpha}^{CW}_{k})]=0\), which by standard arguments has influence function given by \(\frac{\lambda}{E[\lambda]}(\mu_{k}(W_{i})-{\alpha}^{\text{CW}}_{k})\). Since \(\hat{\alpha}^{\text{CW}}_{k}-\alpha^{CW}_{k}=(\hat{\alpha}^{\text{CW}}_{k}-{\alpha}^{\text{CW}}_{k})-({\alpha}^{\text{CW}}_{k} -\tilde{\alpha}^{\text{CW}}_{k})\), we subtract this influence function from the preceding display to obtain \(\eqref{eq:overlap_oracle}\).