The general formula for Pearson's residual is given by: $$ e_i = \frac{y_i - \hat \mu_i}{\sqrt {V(\hat \mu_i)}} $$

But in the multinomial case, the sum of the squared residual, which is the Pearson statistic for goodness of fit, it seems to be written like this:

$$X^2= \sum_{i=1}^N \sum_{c=1}^C \frac{(y_{i,c} - \hat \mu_{i,c})^2}{\hat \mu_{i,c}} \stackrel{?}= \sum e_i^2 $$

I was wondering if someone can show me why this is equal to the sum square of "general" residuals? Basically why $n_i\pi_c(1-\pi_c) = V(\hat\mu_i) \stackrel{?} = \hat \mu_i = n_i\pi_c $?

  • $\begingroup$ @GordonSmyth my GLM professor's lecture+slides. In his slides $\hat \mu$ is replaced by $n_i \hat \pi_{i,c}$ though. $\endgroup$ Commented Jul 15, 2020 at 11:46
  • $\begingroup$ (The equal sign with question mark is what I'm trying to understand... It might not be correct but if so I would want to understand why.) $\endgroup$ Commented Jul 15, 2020 at 11:46
  • $\begingroup$ Also check the answer here and the comment stats.stackexchange.com/questions/47183/… though not very elaborative, seems to suggest the same. $\endgroup$ Commented Jul 15, 2020 at 11:48
  • $\begingroup$ I suggest you do a search for "Poisson trick" to find out more about the relationship between Poisson and multinomial glms. And try simply asking your professor. $\endgroup$ Commented Jul 16, 2020 at 2:10
  • 2
    $\begingroup$ Elementary mathematics tells you that $n_i\pi_c(1-\pi_c) \ne n_i\pi_c $ so that can't be the answer to anything. What happens with multinomial glms is that the residuals for different $y_{nc}$ for the same multinomial observation are negatively correlated, which makes the sum of squared residuals smaller than what you would expect if they were independent. The maths of this is very lengthy to explain. $\endgroup$ Commented Jul 16, 2020 at 2:10

2 Answers 2


The simplest example of a multinomial GLM is binomial regression. Suppose we have fitted a binomial GLM to counts $y_i$, $i=1,\dots, N$. Write $n_i$ for the number of binomial trials and $\pi_i$ for the estimated success probability. The Pearson goodness of fit statistic is $$X^2=\sum_{i=1}^N \frac{(y_i-n_i\pi_i)^2}{n_i\pi_i(1-\pi_i)}$$ where the denominator is the estimated variance var$(y_i)=n_i\pi_i(1-\pi_i)$.

Now let us view this as multinomial data with $c=2$. We define $y_{i1}=y_i$, $y_{i2}=n_i-y_i$, $\pi_{i1}=\pi_i$, $\pi_{i2}=1-\pi_i$. The Pearson goodness of fit statistic for the multinomial data can be written as $$X^2=\sum_{i=1}^N \sum_{c=1}^2 \frac{(y_{ic}-n_i\pi_{ic})^2}{n_i\pi_{ic}}$$

You will find that the two formulae for $X^2$, although they look different, give identical results, despite the fact that the second definition has twice as many terms in the sum as does the first. The reason for this is that $y_{i1}$ and $y_{i2}$ are (perfectly) negatively correlated conditional on $n_i$, which makes the sum of squared residuals smaller than one would otherwise expect. The mathematics of this is very lengthy to write out.

You will notice that the denominator in the second formula is var$(y_{ic})$ if $y_{ic}$ is viewed as Poisson instead of binomial. This arises because treating the $y_{ic}$ as Poisson is equivalent to treating them as binomial when conditioned on $n_i$. You could do a search for "Poisson trick" to find more about the trick of estimating multinomial GLMs via a Poisson GLM.

Multinomial GLMs are multivariate so you can't just plug in the familiar formulae for univariate GLMs, such as the formula you give for a univariate Pearson residual. For multinomial GLMs, the variance function $V(\mu)$ is a matrix rather than a single value. To see the multinomial variance function worked out explicitly, see my paper (Smyth, 1991) on multivariate GLMs. For a multinomial GLM, the closest you can get to residuals is to break up the vector observation up into the individual counts $y_{ic}$ and to form residuals using the marginal binomial distribution for each count. The resulting residuals from the same multinomial vector will be negatively correlated and will need some nuanced interpretation if used in model diagnostics.


Smyth, G.K. (1991). Exponential dispersion models and the Gauss-Newton algorithm. Australian Journal of Statistics. 33, 57–64. https://gksmyth.github.io/pubs/edm-gna.pdf

  • $\begingroup$ I will check this out, though I think I also read something similar in your book (Which I just realized you are one of the authors to, cool). $\endgroup$ Commented Jul 16, 2020 at 7:18
  • $\begingroup$ I just didn’t have time to go over it. Will take me a bit because I have some mess to deal with now. But don’t worry, I’ll come back to this soon. $\endgroup$ Commented Jul 19, 2020 at 8:59
  • $\begingroup$ do you know where I can find this lengthy mathematical derivation? $\endgroup$ Commented Sep 1, 2020 at 7:33

Adding on Gordon Smyth answer, here's the simple derivation for "multinomial" of 2 classes:

$\sum_i \sum_{c=1}^2 \frac{(y_{ic}-n_i\pi_{ic})^2}{n_i\pi_{ic}} =\sum_i \frac{(y_{i1}-n_i\pi_{i1})^2}{n_i\pi_{i1}} + \frac{(y_{i2}-n_i\pi_{i2})^2}{n_i\pi_{i2}}= \sum_i \frac{(y_{i}-n_i\pi_{i})^2}{n_i\pi_{i}} +\frac{(n_i-y_i - n_i(1-\pi_{i}))^2}{n_i(1-\pi_{i})} = \\ \sum_i \frac{(1-\pi_i)(y_{i}-n_i\pi_{i})^2 + \pi_i(y_{i}-n_i\pi_{i})^2}{n_i\pi_{i}(1-\pi_i)} = \sum_i \frac{(y_{i}-n_i\pi_{i})^2}{n_i\pi_{i}(1-\pi_i)} $


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.