Pooling Average Marginal Effects

Why Average Marginal Effects Need Their Own Pooling Step

A regression coefficient is not always the effect a reader wants to interpret. In a nonlinear model, such as logistic regression, the change in predicted probability associated with a one-unit change in a covariate depends on the covariate values and model parameters. An average marginal effect (AME) summarizes those individual changes by averaging over a specified set of cases. With multiple imputation, the AME must be calculated in each completed data set, with its uncertainty, before the estimates are pooled.

The key is to keep the target and scale fixed. For example, if the target is the average change in probability per one-unit increase in a continuous exposure among the analysis cases, calculate that same quantity in every imputation. Pooling the regression coefficients first and then calculating a marginal effect is generally a different operation.

Definition: Let $\hat{\boldsymbol\beta}_j$ denote the estimated model parameters in imputation $j$, and let $g_i(\hat{\boldsymbol\beta}_j)$ be the model-based marginal effect for target case $i$. The imputation-specific average marginal effect is $$ \hat Q_j=\frac{1}{n}\sum_{i=1}^{n}g_i(\hat{\boldsymbol\beta}_j). $$ For a continuous covariate, $g_i$ is usually the derivative of the model’s conditional mean with respect to that covariate. For a binary exposure, it is often the discrete difference between predictions with exposure set to 1 and 0. The target cases and the definition of $g_i$ must be specified and used consistently across imputations.

The AME is an average of case-specific effects, not necessarily the effect evaluated at average covariate values. The two can differ when the model is nonlinear. For a continuous covariate, its units are the outcome-scale units per one unit of that covariate; for a binary exposure and a probability outcome, the effect is a probability difference.

Estimate the Effect and Its Complete-Data Variance

Within each imputation, calculate $\hat Q_j$ and an estimated complete-data variance $U_j$. For a smooth effect function, a common variance estimate uses the delta method. Write the AME as $Q(\boldsymbol\beta)$, let $\widehat{\mathbf V}_j$ be the fitted model’s covariance matrix for $\hat{\boldsymbol\beta}_j$, and let $\mathbf G_j$ be the gradient of the AME with respect to the parameters, evaluated at that imputation’s estimate. Then

U_j=\mathbf G_j^{\mathsf T}\widehat{\mathbf V}_j\mathbf G_j, \qquad \mathbf G_j= \left. \frac{\partial Q(\boldsymbol\beta)}{\partial\boldsymbol\beta} \right|_{\boldsymbol\beta=\hat{\boldsymbol\beta}_j}.

This calculation accounts for uncertainty in the fitted parameters while treating the target covariate values as fixed. It requires a suitable model-based covariance estimate and a differentiable effect function in the relevant parameter region. The delta-method approximation is most credible when the estimator is sufficiently regular and the function is not behaving erratically near the fitted values. If the target is explicitly a superpopulation average and uncertainty in the target covariate distribution is part of the estimand, that additional source of uncertainty must also be addressed; it is not included merely by using the formula above.

Worked Example: Logistic-Regression AME and Delta-Method Variance

Suppose, hypothetically, that a logistic model for a binary outcome has linear predictor $\beta_0+\beta_1x$, and the target cases have $x=0$ and $x=1$. In one imputation, let $\hat\beta_0=0$, $\hat\beta_1=0.8$, and let the estimated covariance matrix be $$ \widehat{\mathbf V}= \begin{pmatrix} 0.04 & -0.005\\ -0.005 & 0.01 \end{pmatrix}. $$ Calculate the average derivative of predicted probability with respect to $x$, and its delta-method variance.

For $p_i=\operatorname{logit}^{-1}(\beta_0+\beta_1x_i)$, the derivative is $g_i=\beta_1p_i(1-p_i)$. At $x=0$, $\hat p_0=0.5000$, so $\hat g_0=0.8(0.5)(0.5)=0.2000$. At $x=1$, $\hat p_1=\operatorname{logit}^{-1}(0.8)=0.6900$, giving $\hat g_1=0.8(0.6900)(0.3100)\approx0.1711$. Thus $$ \hat Q=\frac{0.2000+0.1711}{2}=0.1856. $$

To obtain the gradient, the derivative of $g_i$ with respect to $\beta_0$ is $\beta_1p_i(1-p_i)(1-2p_i)$; with respect to $\beta_1$, it is $p_i(1-p_i)+\beta_1p_i(1-p_i)(1-2p_i)x_i$. Averaging these gradients over the two cases gives $\mathbf G\approx(-0.03250,\ 0.19945)^{\mathsf T}$. Therefore, $$ U=\mathbf G^{\mathsf T}\widehat{\mathbf V}\mathbf G \approx (0.04)(-0.03250)^2 +2(-0.005)(-0.03250)(0.19945) +(0.01)(0.19945)^2 \approx0.0005049. $$

The estimated standard error within this imputation is $\sqrt{U}\approx0.02247$. The AME is about 0.1856 probability units per one-unit increase in $x$, or about 18.56 percentage points per unit. This example assumes the two listed cases define the target and that the fitted covariance matrix is an appropriate estimate of parameter uncertainty.

Pool the Imputation-Specific Effects

After calculating one estimate and one complete-data variance in each of $m$ imputations, pool the AMEs on their effect scale. Let $\bar Q$ be the average of the estimates, $W$ the average complete-data variance, $B$ the between-imputation variance, and $T$ the total variance.

\bar Q=\frac{1}{m}\sum_{j=1}^{m}\hat Q_j,\qquad W=\frac{1}{m}\sum_{j=1}^{m}U_j,\qquad B=\frac{1}{m-1}\sum_{j=1}^{m}(\hat Q_j-\bar Q)^2,\qquad T=W+\left(1+\frac{1}{m}\right)B.

The pooled standard error is $\sqrt{T}$. The term $W$ represents average within-imputation uncertainty; $B$ represents variation among the imputation-specific AMEs. The factor $1+1/m$ accounts for the finite number of imputations when estimating the missing-data contribution to uncertainty.

For a large complete-data degrees-of-freedom setting, a common Rubin degrees-of-freedom approximation is $$ df=(m-1)\left(1+\frac{W}{(1+1/m)B}\right)^2. $$ A confidence interval can then be formed as $\bar Q\pm t_{0.975,df}\sqrt{T}$. When complete-data degrees of freedom are limited, an appropriate finite-sample MI adjustment is needed rather than relying automatically on this large-sample form.

Worked Example: Pooling Four Average Marginal Effects

Consider four hypothetical imputations. Their AMEs are $0.12,\ 0.14,\ 0.13,\ 0.15$, each in probability units per one-unit increase in a continuous exposure. Their complete-data variance estimates are $0.0004,\ 0.0005,\ 0.0004,\ 0.0005$. Assume the same model, target cases, and effect definition were used in all four imputations.

The pooled estimate is $\bar Q=(0.12+0.14+0.13+0.15)/4=0.135$. The average within-imputation variance is $W=(0.0004+0.0005+0.0004+0.0005)/4=0.00045$. The deviations from the pooled estimate are $-0.015,\ 0.005,\ -0.005,\ 0.015$, so $$ B=\frac{(-0.015)^2+(0.005)^2+(-0.005)^2+(0.015)^2}{3} =\frac{0.000500}{3}=0.0001667. $$

Consequently, $$ T=0.00045+\left(1+\frac14\right)(0.0001667) =0.0006583, \qquad SE(\bar Q)=\sqrt{0.0006583}=0.02566. $$

Using the large-sample degrees-of-freedom approximation gives $df=3[1+0.00045/(1.25\times0.0001667)]^2\approx29.96$. With $t_{0.975,29.96}\approx2.043$, the 95% interval is $0.135\pm2.043(0.02566)$, approximately $(0.0826,\ 0.1874)$. In context, the pooled estimated average increase in outcome probability is 13.5 percentage points per exposure unit, with a 95% interval from about 8.26 to 18.74 percentage points per unit. The interval reflects both within- and between-imputation uncertainty under the stated MI assumptions; it is not a 95% probability statement about the fixed effect.

Why the Order of Averaging Matters

For a nonlinear model, averaging case-specific effects is generally not equivalent to evaluating an effect at an average covariate value. This difference defines the target, not a numerical nuisance to be removed. Consider a hypothetical logistic model with $\beta_0=-1$, $\beta_1=1$, and two equally weighted target cases with $x=0$ and $x=2$.

Worked Example: Average Derivative Is Not the Derivative at the Mean

The logistic derivative is $p(x)\{1-p(x)\}\beta_1$. At $x=0$, $p=\operatorname{logit}^{-1}(-1)=0.2689$, so the derivative is $0.2689(0.7311)(1)\approx0.1966$. At $x=2$, $p=\operatorname{logit}^{-1}(1)=0.7311$, giving the same derivative, about $0.1966$. The AME over these two cases is therefore $(0.1966+0.1966)/2=0.1966$.

The mean covariate value is $x=1$. At that value, $p=\operatorname{logit}^{-1}(0)=0.5$, and the derivative is $0.5(0.5)(1)=0.25$. Thus the derivative at the mean is 0.25, not the AME of 0.1966. The former describes one covariate profile; the latter averages effects over the specified two-case target.

Assumptions, Diagnostics, and Alternatives

Pooling AMEs is justified when the imputations are proper for the missing-data analysis, the analysis model and its variance estimates are appropriate, and every $\hat Q_j$ estimates the same defined quantity. The missing-data assumptions and imputation-model adequacy are essential background conditions: Rubin pooling does not repair an unsuitable imputation model. Within each imputation, the delta-method variance also relies on a sufficiently smooth effect function and a credible parameter covariance estimate.

Before pooling, check that each model converged, that the AMEs and their variances are finite, and that the same target cases, weights, exposure contrast, and outcome scale were used in every completed data set. Investigate unusually large AMEs, variances, or between-imputation spread. In logistic models, separation or extreme fitted probabilities can make estimates unstable and the delta approximation questionable. A large $B$ relative to $W$ signals substantial variation among imputed-data estimates; it is a reason to inspect the imputations and model, not a reason to discard between-imputation uncertainty.

The main alternative is to pool regression parameters and then calculate an effect from those pooled parameters. In nonlinear models this need not equal the average of the imputation-specific AMEs, and it does not automatically provide the variance for the AME estimand. Conversely, computing an AME in each imputation but pooling only the point estimates ignores uncertainty from both the fitted models and missing data. The effect and its variance must travel together through the pooling calculation.

An AME depends on the population or case set over which it is averaged and may not describe every individual. It is also model-based: misspecification can distort the predicted responses and their derivatives. These limits do not invalidate the method, but they make a clear estimand definition and consistent computation across imputations indispensable.

Key takeaway: Define one effect and one target, calculate the AME and its complete-data variance within every imputation, then apply Rubin’s rules to those effect-scale quantities. Do not substitute a different estimand by changing the order of averaging and pooling.

Check Your Understanding

Answer each question, then compare with the solutions.

In a logistic model, what is averaged to obtain the AME for a continuous covariate?
Four imputations give AMEs $0.10,\ 0.12,\ 0.14,\ 0.12$. What is $\bar Q$?
If $W=0.0006$, $B=0.0002$, and $m=4$, calculate $T$ and its square-root standard error.
Why can the average derivative differ from the derivative evaluated at the average covariate value?
Name two checks to make before pooling imputation-specific AMEs.

Solutions. (1) Average the model-based derivative of predicted outcome with respect to the covariate over the specified target cases. (2) $\bar Q=(0.10+0.12+0.14+0.12)/4=0.12$. (3) $T=0.0006+(1+1/4)(0.0002)=0.00085$, so $SE=\sqrt{0.00085}\approx0.02915$. (4) In a nonlinear model, the derivative varies with covariate values, so averaging those derivatives generally differs from evaluating the derivative at one average value. (5) For example, check that every model converged and that the same target cases and effect definition were used; also inspect estimates and variances for nonfinite or unusually extreme values.