Inference With Estimated Propensity Scores

Why Estimating the Scores Matters

An inverse-probability-weighted treatment effect uses propensity scores to construct weights. If those scores are estimated from the same sample, their estimation is part of the procedure: changing the sample changes both the weights and the weighted outcome means. A variance calculation that treats the fitted scores as fixed can therefore miss part of the sampling variation, and can misstate uncertainty in either direction.

This tutorial derives a joint sandwich variance for a weighted average treatment effect (ATE). The setup uses a binary treatment $A_i$, outcome $Y_i$, and covariate vector $X_i$ for individual $i$. Let $Z_i$ denote the propensity-model design vector, including an intercept if the model uses one. For concreteness, use a logistic propensity model. The same estimating-equation logic applies to other smooth fitted propensity models when their score equations are included.

Definition: The propensity score model is $e_i(\beta)=P(A_i=1\mid X_i)=\operatorname{expit}(Z_i^\mathsf{T}\beta)$, where $\operatorname{expit}(u)=1/(1+\exp(-u))$. The inverse-probability-weighted Hájek means are $$ \hat\mu_1= \frac{\sum_{i=1}^n A_iY_i/e_i(\hat\beta)} {\sum_{i=1}^n A_i/e_i(\hat\beta)},\qquad \hat\mu_0= \frac{\sum_{i=1}^n (1-A_i)Y_i/\{1-e_i(\hat\beta)\}} {\sum_{i=1}^n (1-A_i)/\{1-e_i(\hat\beta)\}}, $$ and the estimated treatment effect is $\hat\tau=\hat\mu_1-\hat\mu_0$. These are normalized weighted means; the denominators are estimated from the observed sample.

The target here is the ATE, under consistency, conditional exchangeability given $X_i$, and positivity. For causal interpretation, the propensity model must correctly represent the treatment probability given the included covariates. For the variance approximation, assume independent and identically distributed observations, a smooth model and estimating equations, finite relevant second moments, and a nonsingular population derivative matrix. Severe near-zero or near-one fitted scores can make the approximation unstable even when the equations are formally defined.

Stack the Propensity and Outcome Equations

The logistic maximum-likelihood estimate solves the score equation. The two Hájek means solve weighted residual equations: setting each weighted residual sum to zero yields exactly its corresponding weighted mean. Stack these equations so the propensity estimate and outcome means are treated as one joint estimator.

U_i(\theta)= \begin{pmatrix} Z_i\{A_i-e_i(\beta)\}\\[2pt] \dfrac{A_i}{e_i(\beta)}(Y_i-\mu_1)\\[6pt] \dfrac{1-A_i}{1-e_i(\beta)}(Y_i-\mu_0) \end{pmatrix}, \qquad \theta= \begin{pmatrix} \beta\\ \mu_1\\ \mu_0 \end{pmatrix}, \qquad \sum_{i=1}^n U_i(\hat\theta)=0.

The first block is a vector with as many components as $\beta$; the other blocks are scalars. In particular, solving the second equation gives $\hat\mu_1$, and solving the third gives $\hat\mu_0$, using the fitted propensity scores. Define the empirical derivative, often called the bread, and the empirical variability of the equations, called the meat, by

\hat A=\frac{1}{n}\sum_{i=1}^n \left.\frac{\partial U_i(\theta)}{\partial\theta^\mathsf{T}}\right|_{\hat\theta}, \qquad \hat B=\frac{1}{n}\sum_{i=1}^n U_i(\hat\theta)U_i(\hat\theta)^\mathsf{T}.

The joint sandwich estimate of the covariance matrix is

\widehat{\operatorname{Var}}(\hat\theta) =\frac{1}{n}\hat A^{-1}\hat B(\hat A^{-1})^\mathsf{T}.

Let $c$ select the treatment-effect contrast: it has zeros for the propensity coefficients, a 1 for $\mu_1$, and a $-1$ for $\mu_0$. Then $\hat\tau=c^\mathsf{T}\hat\theta$, and its estimated variance is $c^\mathsf{T}\widehat{\operatorname{Var}}(\hat\theta)c$. The standard error is the square root of that quantity. For a large-sample Wald interval, use $\hat\tau\pm z_{1-\alpha/2}SE(\hat\tau)$; its coverage is approximate and relies on the stated regularity conditions.

Why This Sandwich Includes Score Estimation

The key point is that the outcome equations depend on $\beta$ through the weights. Their derivatives with respect to $\beta$ appear in the bread, while the meat retains the observed covariance between the propensity score contributions and the weighted outcome contributions. Thus the calculation does not simply add a separate propensity-score variance to an outcome variance; it also accounts for their covariance.

Justification: Let $A=E\{\partial U_i(\theta_0)/\partial\theta^\mathsf{T}\}$ and $B=E\{U_i(\theta_0)U_i(\theta_0)^\mathsf{T}\}$, where $\theta_0$ solves $E\{U_i(\theta_0)\}=0$. A first-order Taylor expansion of the sample equations around $\theta_0$ gives $$ 0=\frac{1}{n}\sum_i U_i(\theta_0) +A(\hat\theta-\theta_0)+o_p(n^{-1/2}). $$ Therefore $\sqrt n(\hat\theta-\theta_0)=-A^{-1}n^{-1/2}\sum_i U_i(\theta_0)+o_p(1)$. The central limit theorem gives covariance $A^{-1}B(A^{-1})^\mathsf{T}$ for this scaled estimator. Replacing $A$ and $B$ by their empirical versions gives the stated covariance divided by $n$.

For the logistic model, $\partial e_i/\partial\beta=e_i(1-e_i)Z_i$. For example, the derivative of the treated-outcome equation with respect to $\beta$, holding $\mu_1$ fixed, is

\frac{\partial}{\partial\beta^\mathsf{T}} \left\{\frac{A_i}{e_i(\beta)}(Y_i-\mu_1)\right\} = -\frac{A_i\{1-e_i(\beta)\}}{e_i(\beta)} (Y_i-\mu_1)Z_i^\mathsf{T}.

The control-outcome equation also has a propensity derivative. In implementation, calculate the full derivative matrix, including these cross-derivatives, rather than using only the diagonal derivatives with respect to $\beta,\mu_1,\mu_0$. The joint inverse then propagates uncertainty in the fitted scores into both weighted means.

Fit the propensity model.
Estimate $\hat\beta$ from the treatment score equation and calculate each $\hat e_i=e_i(\hat\beta)$.

Estimate the weighted means.
Solve the two weighted residual equations for $\hat\mu_1$ and $\hat\mu_0$, then form $\hat\tau=\hat\mu_1-\hat\mu_0$.

Build the joint sandwich.
Evaluate each observation's stacked equation and derivative at $\hat\theta$, form $\hat A$ and $\hat B$, and calculate the contrast variance.

Check stability and interpretation.
Inspect fitted scores, weights, and the derivative matrix. Report the effect and standard error in outcome units, and state the causal and sampling assumptions.

Worked Calculations and Interpretation

Worked Example: Calculating a Hájek ATE From Hypothetical Data

Suppose a hypothetical sample has two treated individuals with outcomes 8 and 10 and fitted scores 0.4 and 0.5. Two control individuals have outcomes 5 and 6 and fitted scores 0.4 and 0.6. Assume these fitted scores came from a propensity model estimated on the full sample; this small table is only for illustrating the weighted point estimate.

For the treated group, the weighted outcome total is $8/0.4+10/0.5=20+20=40$, and the sum of weights is $1/0.4+1/0.5=2.5+2=4.5$. Thus $\hat\mu_1=40/4.5=8.889$.

For the controls, the weighted outcome total is $5/(1-0.4)+6/(1-0.6)=5/0.6+6/0.4=8.333+15=23.333$. The weights sum to $1/0.6+1/0.4=1.667+2.5=4.167$, so $\hat\mu_0=23.333/4.167=5.600$.

The weighted effect is $\hat\tau=8.889-5.600=3.289$ outcome units. This calculation alone does not provide its standard error: the joint estimating equations also require the propensity score contributions and the derivatives of both outcome equations with respect to the fitted propensity parameters.

Worked Example: The Nuisance-Parameter Correction in a Local Approximation

Consider a hypothetical one-dimensional local approximation to a stacked fit with a nuisance propensity parameter $\beta$ and an effect parameter $\tau$. Suppose the equation derivative matrix, in the order $(\beta,\tau)$, is

\hat A= \begin{pmatrix} -4&0\\ 1&-2 \end{pmatrix}.

Let $s_i$ be the propensity score contribution and $\psi_i$ the effect-equation contribution. Suppose their empirical second moments are $n^{-1}\sum s_i^2=4$, $n^{-1}\sum\psi_i^2=0.09$, and $n^{-1}\sum s_i\psi_i=-0.30$. From $\hat A\,IF_i=- (s_i,\psi_i)^\mathsf{T}$, the first row gives $IF_{\beta,i}=s_i/4$. The second gives $IF_{\tau,i}=(\psi_i+IF_{\beta,i})/2=\psi_i/2+s_i/8$.

Consequently, the estimated variance per observation for the effect influence contribution is

\frac{1}{n}\sum IF_{\tau,i}^2 =\frac{1}{4}(0.09)+\frac{1}{64}(4)+\frac{1}{8}(-0.30) =0.0225+0.0625-0.0375=0.0475.

For a hypothetical sample of $n=100$, the estimated variance is $0.0475/100=0.000475$, and the standard error is $\sqrt{0.000475}=0.02179$ effect units. If the propensity contribution were incorrectly ignored in this illustration, the variance would instead be $0.0225/100=0.000225$, giving standard error $0.01500$. The difference arises from both propensity-score variability and its covariance with the effect equation. In other data, the covariance can have the opposite sign, so ignoring score estimation need not always make the standard error smaller.

Worked Example: Constructing a Wald Interval

Suppose a hypothetical analysis estimates an ATE of 2.40 points with a joint-sandwich standard error of 0.80 points. For a 95% large-sample Wald interval, use $z_{0.975}=1.96$:

2.40\pm1.96(0.80)=2.40\pm1.568, \qquad [0.832,\ 3.968]\text{ points}.

The interval estimates a positive average effect of about 0.832 to 3.968 points under the model and identification assumptions. Its 95% coverage is a repeated-sampling property of the procedure under those assumptions; it is not a 95% probability statement about this fixed interval containing the true effect.

Diagnostics, Alternatives, and Limitations

The sandwich formula is not a remedy for a poor propensity model or unstable weights. Review the fitted score distribution separately for treated and control individuals, identify extreme weights, and examine whether a few observations dominate the weighted means. Also check convergence and whether the empirical derivative matrix is singular or nearly singular. A near-singular bread matrix can produce an unstable variance even if the point estimate appears numerically available.

The formula above uses the logistic likelihood score for $\beta$, together with normalized Hájek outcome equations. If the propensity model is fitted with penalization, calibration, or a different estimating procedure, its actual estimating equations must replace the logistic score equation in the stack. Treating a penalized estimate as if it solved the ordinary likelihood score equation does not generally give the right sandwich. Likewise, a variance derived for normalized Hájek means should not be silently applied to an unnormalized Horvitz–Thompson estimator; their estimating equations differ.

The main shortcut is to calculate an outcome-only weighted variance while holding $\hat e_i$ fixed. This can be useful when the scores truly are fixed by design or were obtained from an independent source under an appropriate sampling calculation, but it does not generally describe uncertainty when scores were estimated in the analysis sample. A joint sandwich is a large-sample analytic approach; resampling the entire fitting procedure is another approach, but it has its own assumptions and is not a substitute for checking the design, model, and weight stability.

Mistake: treating fitted scores as known. The scores were estimated from the observed treatment assignments, so their variability and covariance with the weighted outcome equations can matter. Include the propensity score equations in the stack.
Mistake: adding only a propensity variance term. The joint calculation also includes cross-covariance between score and outcome contributions. Use the full matrix sandwich.
Mistake: interpreting robust variance as protection against bias. The sandwich estimates sampling variability around the estimator's large-sample target. It does not correct confounding from omitted covariates, model misspecification, or positivity problems.

Key takeaway: When propensity scores and weighted outcome means are estimated from the same sample, put their estimating equations in one stack. The derivative matrix propagates score estimation into the means, and the empirical meat preserves covariance among all equation contributions.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why are the two weighted residual equations equivalent to the normalized weighted means?
For the logistic model, what is $\partial e_i/\partial\beta$, and why does it matter for the sandwich?
A point estimate is 1.5 outcome units and its joint-sandwich standard error is 0.5 units. Calculate a 95% Wald interval using 1.96.
What can happen if extreme fitted scores make the empirical derivative matrix nearly singular?
Does a sandwich standard error protect an ATE estimate from bias caused by an omitted confounder? Explain.

Solutions. (1) Setting $\sum_i A_i(Y_i-\mu_1)/e_i=0$ and solving for $\mu_1$ gives the treated weighted outcome total divided by the treated weight total; the control equation gives the analogous control ratio. (2) The derivative is $e_i(1-e_i)Z_i$. It determines how the weights and weighted means change with $\beta$, so its cross-derivatives enter the bread. (3) $1.5\pm1.96(0.5)=1.5\pm0.98$, giving $[0.52,2.48]$ outcome units. (4) Inverting the bread may amplify small numerical or sampling changes, producing an unstable, very large, or unreliable standard error. (5) No. The sandwich describes sampling variability under its regularity conditions; it does not remove confounding bias or make an incorrectly specified causal adjustment valid.