Individual marginal effects for generalized linear models

Data (James et al., 2021, sec. 4.1):

data("Default", package = "ISLR2")
Def0 <- Default |> transform(balance = balance / 10000, income = income / 10000)
set.seed(123)
(Def0 <- Def0[sample(nrow(Def0), 2500), ]) |> str()

'data.frame':   2500 obs. of  4 variables:
 $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ student: Factor w/ 2 levels "No","Yes": 1 1 2 1 2 2 1 1 1 1 ...
 $ balance: num  0.0559 0.0906 0.1807 0.1615 0.1229 ...
 $ income : num  6.12 4.73 1.76 5.52 1.22 ...

default (\(y\)), student (\(x_{1}\)), balance (\(x_{2}\)), income (\(x_{3}\))

Model

A probit-link GLM for the probability of default with

default ~ student * (balance + income)

Individual marginal effects for generalized linear models

Focus

The individual marginal effects of balance on the probability of default \[ \partial P(y = \text{\rm "Yes"} \mid x_1, x_2, x_3) / \partial x_2 \] at the covariate settings with average income and average balance by students status

(avgs <- aggregate(cbind(income, balance) ~ student, data = Def0, mean))

  student   income    balance
1      No 3.991203 0.07722011
2     Yes 1.797029 0.10000943

Individual marginal effects for generalized linear models

library("marginaleffects")
Def0_fit <- glm(default ~ student * (balance + income),
                    family = binomial(probit),
                    data = Def0)
slopes(Def0_fit, variables = "balance", newdata = avgs)

 Estimate Std. Error    z Pr(>|z|)   S   2.5 % 97.5 %
   0.0110     0.0110 1.00    0.316 1.7 -0.0105 0.0325
   0.0898     0.0764 1.18    0.240 2.1 -0.0599 0.2396

Term: balance
Type: response
Comparison: dY/dX

Estimator performance

10000 samples under the ML fit (Def0_fit), estimating the same focus using various methods

Inference

Wald inference is anti-conservative \(\implies\) spuriously strong conclusions

Inference

Median bias correction \(+\) Hull-based confidence \(\implies\) near nominal finite-sample coverage

Mean and median bias reduction

Mean bias reduction

Reduced mean-bias estimators are well-established in parametric modelling¹, being particularly attractive in the estimation and inference from models for categorical data.²

Median bias reduction

Reduced median-bias estimators³ are also attractive in a wide range of problems because they

• are invariant under componentwise monotone transformations

• correspond to a zero-level confidence interval

• tend to improve the behaviour of standard plug-in asymptotic inference (e.g. Wald-type)

Mean and median bias reduction

Five independent observations from an exponential distribution with rate \(1\)

We estimate the negative rate, denoted by \(\psi\), whose true value is \(-1\).

Focused median bias reduction

Focus parameter

A smooth scalar transformation of the model parameters that is of interest.

Examples: marginal effects in GLMs, Mahalanobis distances, ordinal superiority measures in ordinal regression, joint indirect effects in multiple-mediator models, distributional quantiles.

Median bias correction for focus parameters

Current general methods for median bias reduction¹ are implicit, requiring the solution of nonlinear adjusted likelihood equations.

For arbitrary focus parameters, we need to set the focus as a component of the parameter vector, which typically involves tedious algebra and bespoke implementations.

Goal

An explicit median bias corrected estimator for general focus parameters and models.

Setup

Model, log-likelihood, and ML estimator

A regular parametric model with log-likelihood \(l(\theta)\), with \(\theta\in\Theta\subseteq\mathbb R^p\)

Focus parameter

\(\psi = h(\theta)\) is the focus parameter, where \(h:\Theta\to\mathbb R\) is a smooth scalar transformation.

ML estimator of \(\theta\)

The maximum likelihood estimator is \(\hat\theta = \arg\max_{\theta \in \Theta} \ell(\theta)\).

ML estimator of \(\psi\)

The plug-in estimator of \(\psi\) is \(\hat\psi=h(\hat\theta)\).

Explicit median bias correction: Oracle estimator

Let \(W = (\hat\psi - E(\hat\psi)) / \{\mathop{\mathrm{var}}(\hat\psi)\}^{1/2}\). The roots of the leading term of a normalizing transformation of \(W\) based on a Cornish-Fisher expansion give \[ \begin{aligned} \tilde\psi^{(a)} = \underbrace{\overbrace{\hat\psi - B_{(\psi)}}^{\text{mean bias correction}} + \overbrace{\frac{1}{6}\frac{K_3}{K_2}}^{\text{skewness adjustment}}}_{\displaystyle\tilde\psi^{(o)}} + \overbrace{O_p(n^{-3/2})}^{\text{lower order terms}} \end{aligned} \] where \[ B_{(\psi)}\equiv B_{(\psi)}(\theta)=E(\hat\psi)-\psi, \quad K_2 \equiv K_2(\theta)=\mathop{\mathrm{var}}(\hat\psi), \quad K_3 \equiv K_3(\theta) = {\rm cum}_3(\hat\psi) \, . \]

Then, \(\tilde\psi^{(o)}\) is third-order median unbiased with \(P(\tilde\psi^{(o)} \le \psi) = 1/2 + O(n^{-3/2})\).

Non-oracle estimators

\(B_{(\psi)}\), \(K_2\), and \(K_3\) typically:

• depend on the unknown parameter \(\theta\)

• are not available in closed form

Non-oracle estimators: Closed forms

If \(B_{(\psi)}\), \(K_2\), and \(K_3\) are available in closed form, and \(\hat{B}_{(\psi)} = B_{(\psi)}(\hat\theta)\), \(\hat{K}_2 = K_2(\hat\theta)\), and \(\hat{K}_3 = K_3(\hat\theta)\), then \[ \tilde\psi = \hat\psi - \hat{B}_{(\psi)} + \frac{1}{6}\frac{\hat{K}_3}{\hat{K}_2} \] is third-order median unbiased.

Squared Mahalanobis distance

\(Y_1,\ldots,Y_n \stackrel{\mathrm{iid}}{\sim} N_p(\mu,\Sigma)\), with \(\Sigma\) positive definite.

Focus parameter for a fixed point \(\mu_0\)

\(\psi=(\mu-\mu_0)^\top\Sigma^{-1}(\mu - \mu_0)\)

Maximum likelihood estimator

\(\hat\psi = (\bar Y-\mu_0)^\top S^{-1}(\bar Y-\mu_0)\)

\(B_{(\psi)}\), \(K_2\), \(K_3\)

For \(\nu_k = n - p - k\) and \(\nu_6 > 0\) \[ \begin{aligned} B_{(\psi)} = \frac{m_1}{\nu_2} - \psi \,, \quad K_2 = \frac{m_2}{\nu_2\nu_4} - \frac{m_1^2}{\nu_2^2} \,, \quad K_3 = \frac{m_3}{\nu_2\nu_4\nu_6} - 3 \frac{m_1m_2}{\nu_2^2\nu_4} + 2\frac{m_1^3}{\nu_2^3} \end{aligned} \] where \(m_j\) are polynomials of \(\psi\) only.

Squared Mahalanobis distance

Non-oracle estimators: Asymptotic approximations

We can write

\[ \begin{array}{ll} B_{(\psi)} = b_{(\psi)} + O(n^{-2}) \,, \quad & b_{(\psi)}= {b}_{(\theta)}^\top\dot{\psi}+ \frac{1}{2} \mathop{\mathrm{trace}}({V}\ddot{\psi}) \\ K_2 = \kappa_2 + O(n^{-2}) \,, \quad & \kappa_2 = {\omega}^\top \dot{\psi}\\ K_3 = \kappa_3 + O(n^{-3}) \,, \quad & \kappa_3 = -6 \sum_{t = 1}^p \omega_t {\omega}^\top \left(\frac{1}{3} {P}_t + \frac{1}{2} {Q}_t\right) {\omega}+ 3 {\omega}^\top \ddot{\psi}{\omega}\\ \end{array} \]

Then, \[ \tilde\psi = \hat\psi - \hat{b}_{(\psi)}+ \frac{1}{6}\frac{\hat{\kappa}_3}{\hat{\kappa}_2} \] is third-order median unbiased.

Non-oracle estimators: Asymptotic approximations

where

\({P}_t = \mathop{\mathrm{E}}({s} {s}^\top s_t)\) and \({Q}_t = -\mathop{\mathrm{E}}({j} s_t)\) with \({s} = \nabla l({\theta}; Y)\) and \({j}= -\nabla \nabla^\top l({\theta}; {{Y}})\)

\(\dot{\psi}\) and \(\ddot{\psi}\) can be computed either analytically or numerically at \(\hat\theta\).

\({P}_t\), \({Q}_t\) are available in the mean bias reduction literature for a wealth of models, and if not they can be approximated through Monte-Carlo simulation at \(\hat\theta\).

Notes

When the focus parameter is a component of \(\theta\), the estimator is a one-step approximation to implicit median bias reduction.

For general focus parameters, it avoids the need of an explicit nuisance parameterization and solving adjusted likelihood equations.

In contrast to the implicit, the explicit estimator is not exactly equivariant, although monotone transformation preserves its median-bias order.

Instead of \(\hat\theta\) we can use \(\tilde\theta\) with \(O(n^{-2})\) mean bias by setting \(b_{(\psi)}= \mathop{\mathrm{trace}}({V}\ddot{\psi}) / 2\).¹

Hull-based Confidence (Kuchibhotla et al., 2024)

Split the sample into independent batches \(D_1,\ldots, D_B\), and compute a focus estimator \(\check\psi^{(j)}\) on each batch.

\[ \mathrm{CI}^{\mathrm{HulC}}(\psi) = \left[ \min_{1 \le j \le B} \check\psi^{(j)}, \; \max_{1 \le j \le B} \check\psi^{(j)} \right] \]

The miscoverage probability satisfies \[ P\left(\psi \notin \mathrm{CI}^{\mathrm{HulC}}\right) \le u(B;\Delta) = \left(\tfrac12-\Delta\right)^B + \left(\tfrac12+\Delta\right)^B \] where \(\Delta\in[0,1/2]\) is a uniform upper bound on the median bias of the batch estimators.

\(\implies\) Choose \(B\) such that \(u(B;\Delta)\le \alpha\).

If \(\Delta = 0\), \(B_0 = \Bigl\lceil \log_2(2/\alpha)\Bigr\rceil\) guarantees coverage of at least \(1 -\alpha\).

For third-order median-unbiased batch estimators, the coverage of \(\mathrm{CI}^{\mathrm{HulC}}(\psi)\) with \(\Delta=0\) and calibrated randomisation between \(B_0 - 1\) and \(B_0\) is sixth-order accurate.

Individual marginal effects for generalized linear models

Remarks

Explicit correction

Explicit median bias correction for arbitrary smooth scalar focus parameters, without nuisance reparameterization or further optimization, using a convenient reference parameterization

Broad applicability

Marginal effects, Mahalanobis distances, distributional quantiles, mediation effects, post-selection inference, …

Strong inferential guarantees

Combining $with HulC delivers near-nominal finite-sample coverage, without standard-error estimation.

Ongoing work

Extensions to estimating functions, including GEE and causal inference settings