[R] Expressing a multinomial GLM as a series of binomial GLMs

[Charles Berry]

Scherber, Christoph <cscherb1 <at> gwdg.de> writes:

> 
> Dear all,
> 
> I am trying to express a multinomial GLM (using nnet) as a series of GLM
models.
> 
> However, when I compare the multinom() predictions to those from GLM, I
see differences that I can´t
> explain. Can anyone help me out here?
> 
> Here comes a reproducible example:
> 
> ##
> # set up data: (don´t care what they are, just for playing)
> set.seed(0)
> cats=c("oligolectic","polylectic","specialist","generalist")
> explan1=c("natural","managed")
> explan2=c("meadow","meadow","pasture","pasture")
> multicats=factor(sample(cats,replace=T,100,prob=c(0.5,0.2,0.1,0.5)))
> multiplan1=factor(rep(explan1,50))
> multiplan2=factor(rep(explan2,25))
> 
> ########################
> library(nnet)
> m2=multinom(multicats~multiplan1)
> 
> # predictions from multinomial model
> predict(m2,type="probs")
> 
> ########################
> # now set up contrasts for response variable "multicats" (which has 4 levels):

[snip - Christoph's comparison]

Doing the obvious comparison:

ggen.preds <- 
    sapply( levels(multicats), 
            function(x) predict(glm(I(multicats==x)~multiplan1,
                         family=binomial),type="response"))

max(abs(ggen.preds-predict(m2,type="probs")))
## [1] 1.349607e-06

---

The predictions are the same - up to numerical issues in the algorithms.

HTH,

Chuck