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, Hull, and Kolesár (2024), we generalize the setup to allow for sampling weights ω_i² (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 ∈ {0, 1, …, 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 𝒢 = {z′γ: γ ∈ ℝ^{1 + dim (W_i)}}.

We assume that μ_k(W_i) := E[Y_i(k) ∣ W_i] ∈ 𝒢, so that we may write μ_k(W_i) = W_i′α_k for some vectors α_k, k = 0, …, K. The average treatment effect (ATE) conditional on W_i is then given by τ(W_i) = Z_i′(α_k − α₀), and α_k correspond to the coefficients in the interacted regression with U̇_i conditionally mean zero. The uninteracted partially linear (PL) estimator is given by estimating by weighted least squares (WLS), yielding β̂ = (∑_iω_i²Ẋ_iẊ_i′)⁻¹∑_iω_i²Ẋ_iY_i, where Ẋ is the sample residual from WLS regression of X_i onto Z_i. By Proposition 1 in Goldsmith-Pinkham, Hull, and Kolesár (2024), the population analog of β̂, β, has the decomposition where Λ_i = E[X̃_iX̃_i′]⁻¹E[X̃_iX_i ∣ W_i], and X̃_i is the population analog of Ẋ_i, the population residual from regressing X_i onto W_i. Let α̂_k denote the WLS estimates based on . By construction, the sample residuals from estimating and Z_i are both orthogonal to Ẋ_i. As a result, we obtain the exact decomposition where Γ̂ is a matrix with columns γ̂_k = α̂_k − α̂₀, Λ̂_i = E_n[Ẋ_iẊ_i′]⁻¹Ẋ_iX_i′, and E_n[A_i] = ∑_iω_i²A_i/∑_iω_i² denotes the weighted sample mean.

To compute this decomposition, we don’t need to explicitly compute Λ̂_i. Instead, we use the identity where δ̂_kk is a WLS estimator of the system of regressions

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 μ_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[Λ̂_i] = I_k, mimicking the population decomposition. If the propensity score doesn’t satisfy p_k ∈ 𝒢, the implied estimate of Λ(w), may not be positive definite; in particular, the diagonal elements may not all be positive, in line with Proposition 1 in Goldsmith-Pinkham, Hull, and Kolesár (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 using only observations with D_i ∈ {0, k}. In other words, it estimates the regression among observations with D_i ∈ {0, k}. This estimator weights the treatment effects using the variance-minimizing weighting scheme given in Corollary 1 in Goldsmith-Pinkham, Hull, and Kolesár (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 where, by default, the probabilities π_k correspond to the marginal probability E_n[X_ik] in the dataset. The propensity scores 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, Hull, and Kolesár (2024), this weighting scheme minimizes the average variance, under homoskedasticity, across all treatment comparisons—comparisons of outcomes under treatment k vs treatment ℓ, if we draw the treatments k and ℓ independently from the marginal treatment distribution (π₀, …, π_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], …, E_n[X_iK]).

Derivations

We first derive . 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 with the normalization θ₀ = 0. In the second step, we use the moment condition to obtain estimates α̂_k^CW, and set $\hbcw_k=\hat\alpha^{\text{CW}}_{k}-\hat\alpha^{\text{CW}}_{0}$.

Let By equation (6.6) in Newey and McFadden (1994), the influence function of this two-step estimator is given by where ψ_i(θ) is the influence function of the multinomial logit estimator θ̂, and M_k(θ) is the derivative of wrt θ.

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 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 Since ψ_i(θ) = −H(θ)⁻¹S_i(θ), this yields with M̂_k(θ) = (∑_i(ζ − e_k) ⊗ Z_i ⋅ ζ_kX_ik(Y_i − α_λ^C, k)), and $\hat{H}({\theta})= \sum_{i}(\diag(\pi(W_{i}))-\pi(W_{i})\pi(W_{i})')\otimes Z_{i}Z_{i}'$. When π_k = 1/(K + 1) this formula reduces to that in Theorem 1 in Li and Li (2019).

Next, we show . Note that it follows from Proposition 2 in Goldsmith-Pinkham, Hull, and Kolesár (2024) that the efficient influence function is given by where r_k = r_k(W_i) = E[X_ik ∣ W_i, X_ik + X_i0 = 1]. The result then follows since X̂_ik is an estimator of X_ik − r_k(W_i).

Finally, we show . The derivative of the moment condition with respect to π_k = p_k (assuming correct specification of the propensity score) is given by where we write λ for $\lcw(W_{i})$. Since p_k is a projection, by Proposition 4 in Newey (1994), the influence function for α̂_k^CW is given by Next, as noted in Abadie, Imbens, and Zheng (2014), we can view α̃_k^CW = ∑_iλμ_k(W_i)/∑_iλ as an estimator of α_k^CW based on the moment condition E[λ(μ_k(W_i) − α_k^CW)] = 0, which by standard arguments has influence function given by $\frac{\lambda}{E[\lambda]}(\mu_{k}(W_{i})-{\alpha}^{\text{CW}}_{k})$. Since α̂_k^CW − α_k^CW = (α̂_k^CW − α_k^CW) − (α_k^CW − α̃_k^CW), we subtract this influence function from the preceding display to obtain .