Proof that the Sample Mean of the Predicted Response from a Generalized Linear Model (GLM) Fit Via MLE Equals the Sample Mean of the Response

Question

Is there a simple proof that the sample mean of the predicted response from a Generalized Linear Model (GLM) fit via maximum likelihood estimation equals the sample mean of the response? That is, a proof that

$$ \frac{1}{n} \sum_{i = 1}^{n} g^{-1}(x_{i}^{T}\hat{\beta}) = \frac{1}{n} \sum_{i = 1}^{n} Y_{i} $$

where $g$ is the link function and $\hat{\beta}$ is the MLE for $\beta$. This is trivially true for the Gaussian GLM, so I mean for other GLMs like logistic and Poisson regression.

I first noticed this for logistic regression, e.g.

set.seed(42)

n <- 100

x <- rnorm(n)

y <- rbinom(n, 1, plogis(x))

mod <- glm(y ~ x, family = binomial(link = "logit"))

phat <- predict(mod, type = 'response')

sprintf("%.20f", mean(y))
sprintf("%.20f", mean(phat))

[1] "0.47999999999999998224"
[1] "0.48000000000192888372"

where the mean of the predicted response and the mean of the outcome are equal to within floating-point fuzz.

The same is not true before convergence to the MLE:

mod <- glm(
  y ~ x, family = binomial(link = "logit"), 
  control = list(maxit = 2, trace = FALSE)
)

phat <- predict(mod, type = 'response')

sprintf("%.20f", mean(y))
sprintf("%.20f", mean(phat))

Warning: glm.fit: algorithm did not converge
[1] "0.47999999999999998224"
[1] "0.48023472280308265869"

Answer 1

I came back to this, and realized the result is a direct consequence of a first order optimality condition of maximum likelihood estimation for a distribution in the exponential family.

Consider a distribution in the exponential family with log-probability density / mass function $$ \log f(y; \theta) = \frac{y \theta - b(\theta)}{\phi} + c(y; \phi). $$ Then the log-likelihood of the data is $$ \ell(\boldsymbol{\theta}) = \sum_{i = 1}^{n} \left \{ \frac{Y_{i} \theta_{i} - b(\theta_{i})}{\phi}+ c(Y_{i}; \phi) \right \}. $$ Assuming that under the link function $\theta_{i} = g(\mu(x_{i})) = \beta_{0} + \beta_{1}x_{i}$, the log-likelihood becomes $$ \ell(\beta_{0}, \beta_{1}) = \sum_{i = 1}^{n} \left\{ \frac{Y_{i} (\beta_{0} + \beta_{1} x_{i}) - b(\beta_{0} + \beta_{1}x_{i}) }{\phi} + c(Y_{i}; \phi) \right\}. $$ Using the first-order optimality condition for $\beta_{0}$ under maximum likelihood estimation, we have $$ \frac{\partial}{\partial \beta_{0}} \ell(\beta_{0}, \beta_{1}) = \sum_{i = 1}^{n} \frac{Y_{i} - b'(\beta_{0} + \beta_{1}x_{i})}{\phi} = 0. $$ Rearranging the left hand side, $$ \frac{1}{n} \sum_{i = 1}^{n} Y_{i} = \frac{1}{n} \sum_{i = 1}^{n} b'(\beta_{0} + \beta_{1}x_{i}). $$ Recalling that for a distribution in the exponential family, the derivative of the log-partition function gives the expected value of the random variable, $$ b'(\theta(x)) = E[Y \mid X = x], $$ and thus also gives the inverse link function, i.e. $b'(\theta(x)) = g^{-1}(\theta(x))$. Substituting into the result from the first order optimality condition, we therefore find $$ \frac{1}{n} \sum_{i = 1}^{n} Y_{i} = \frac{1}{n} \sum_{i = 1}^{n} g^{-1}(\beta_{0} + \beta_{1}x_{i}) \implies \bar{Y} = \overline{g^{-1}(\beta_{0} + \beta_{1}x)}, $$ establishing the result.