Multinomial GLM Pearson residuals 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 $?
 A: 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.
Reference
Smyth, G.K. (1991). Exponential dispersion models and the Gauss-Newton algorithm. Australian Journal of Statistics. 33, 57–64.
http://www.statsci.org/smyth/pubs/edm-gna.pdf
A: 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)}
$
