Projection for discrete exponential dispersion (DED) families

The families currently supported by projpred (gaussian(), binomial(), brms::bernoulli(), poisson()) belong to the class of "exponential dispersion (ED) families" (Jørgensen, 1987; see also section "Remarks" below) used in GLMs (see McCullagh and Nelder, 1989, chapter 2, equation (2.4) at the beginning of section 2.2.2, p. 28) (and also used in GLMMs, GAMs, and GAMMs). As shown by Piironen et al. (2020, section 3.4), the projection onto an ED-family submodel can be performed quite easily by fitting to the fit of the reference model, at least when regarding only the projection onto the location-specific parameter vector and disregarding the projection onto the dispersion parameter.

Jørgensen (1987) also presents the class of "discrete exponential dispersion families" (hereafter abbreviated by "DED families"). For example, the negative binomial distribution is a DED family. DED families are closely related to ED families (Jørgensen, 1987), but strictly speaking, the class of DED families is not a subset of the class of ED families, as can be seen from comparing equations (1.1) and (2.12) of Jørgensen (1987). However, the projection of a reference model onto the location-specific parameters of a DED-family submodel can be performed just as easily as for an ED-family submodel because for DED families, equation (12) of Piironen et al. (2020) simplifies to

$\theta_{c} = \mathop{\mathrm{argmax}}\limits_{\theta \in \Theta} \sum\limits_{i = 1}^{N} \left((\mu^{*}_{c, i} \, \mu_i(\dot{\theta}) - K(\phi) B(\mu_i(\dot{\theta}))) + \mathbb{E}(H(\tilde{y}_i, \phi)|I_c)\right)$

where I have adapted the notation a bit:

• $\theta_{c}$ corresponds to $\theta_{\perp}$ from Piironen et al. (2020).
• I have split up $\theta$ into the location-specific parameter vector $\dot{\theta} \in \dot{\Theta}$ and the dispersion parameter $\phi \in (0, \infty)$, i.e., $\theta = (\dot{\theta}^{\scriptscriptstyle \mathsf{T}}, \phi)^{\scriptscriptstyle \mathsf{T}}$ and accordingly $\Theta = \{\theta = (\dot{\theta}^{\scriptscriptstyle \mathsf{T}}, \phi)^{\scriptscriptstyle \mathsf{T}}: \dot{\theta} \in \dot{\Theta}, \phi \in (0, \infty)\}$.
• $N$ corresponds to $n$ from Piironen et al. (2020). (I picked $N$ to agree with the notation from projpred's documentation.)
• $\mu^{*}_{c, i}$ corresponds to $\mu^{*}_i$ from Piironen et al. (2020).
• $\mu_i(\cdot)$ corresponds to $\eta_i(\cdot)$ from Piironen et al. (2020). (I picked $\mu$ instead of $\eta$ to avoid confusion with the linear or latent predictor.)
• $\mathbb{E}(\cdot|I_c)$ corresponds to $\mathbb{E}_{\tilde{y}_i | I_c}(\cdot)$ from Piironen et al. (2020). (The GitHub Markdown math mode has some issues with the latter notation when used above.)
• $K(\phi)$ corresponds to $\lambda$ from Jørgensen (1987). (I chose this functional notation in analogy to $A(\phi)$ from Piironen et al. (2020).)

When calculating the gradient of the right-hand side expression above (which is to be maximized), one can see that the term $\mathbb{E}(H(\tilde{y}_i, \phi)|I_c)$ is unimportant for the maximization with respect to the location-specific parameter vector $\dot{\theta}$ and so the projection solution with respect to $\dot{\theta}$ is obtained by a maximum-likelihood fit to the fit of the reference model, analogously to what Piironen et al. (2020, section 3.4) derived for ED families. For example, in case of a negative binomial submodel without multilevel and without additive terms, MASS::glm.nb() could be used for calculating the projection solution with respect to $\dot{\theta}$.

For the dispersion parameter, things get more complicated: The projection solution for $\phi$ is given by

$\phi_{c} = \mathop{\mathrm{argmax}}\limits_{\phi \in (0, \infty)} \sum\limits_{i = 1}^{N} \left((\mu^{*}_{c, i} \, \mu_i(\dot{\theta}_c) - K(\phi) B(\mu_i(\dot{\theta}_c))) + \mathbb{E}(H(\tilde{y}_i, \phi)|I_c)\right)$

with $\dot{\theta}_c$ denoting the projection solution with respect to $\dot{\theta}$ calculated as described above (fitting to the fit of the reference model). For the negative binomial distribution, I haven't been able so far to simplify this $\phi_c$ expression to a more tractable one because $\mathbb{E}(H(\tilde{y}_i, \phi)|I_c)$ essentially consists of expectations of log-Gamma function values. But perhaps future work for the negative binomial family could continue here, possibly using an approximation for this more complicated $\mathbb{E}(H(\tilde{y}_i, \phi)|I_c)$ part.

Remarks

1. Jørgensen (1987) himself only uses the term "exponential dispersion models", but the discussion for that article mentions the term "exponential dispersion (ED) families", so I think it should be fine to use the term "ED families" also in the context of projpred.
2. With regard to the GLM definition from McCullagh and Nelder (1989, chapter 2) that is also cited by Piironen et al. (2020), it is not really correct to say that "GLMs have their observation model in the exponential family" (Piironen et al., 2020, section 3.5) because the class of exponential families is different from the class of ED families: For a known dispersion parameter, ED families are indeed (special, i.e., linear, see Jørgensen, 1987, section 2) exponential families, but in general, this is not the case if the dispersion parameter is unknown (McCullagh and Nelder, 1989, chapter 2). I admit, however, that the wording of McCullagh and Nelder (1989, chapter 2, beginning of section 2.2.2, p. 28) is mistakable. Others seem to have been confused by this as well, see, e.g., here and here.

References

Jørgensen, B. (1987). Exponential Dispersion Models. Journal of the Royal Statistical Society. Series B (Methodological), 49(2), 127–162.

McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (D. R. Cox, D. V. Hinkley, N. Reid, D. B. Rubin, & B. W. Silverman, Eds.; 2nd ed., Vol. 37). Chapman & Hall.

Piironen, J., Paasiniemi, M., & Vehtari, A. (2020). Projective inference in high-dimensional problems: Prediction and feature selection. Electronic Journal of Statistics, 14(1), 2155–2197. https://doi.org/10.1214/20-EJS1711

EDIT: In remark 2, I emphasized that I'm referring to the GLM definition from McCullagh and Nelder (1989, chapter 2). In general, GLMs can have different definitions.

The issue names no source files, tests, or implementation entry points. Start by locating projpred's existing projection support for Gaussian, binomial, Bernoulli, and Poisson families, then compare it with the proposed MASS::glm.nb() approach. Done would require a defined implementation and validation path for discrete exponential dispersion projections, including the unresolved dispersion parameter.