[R] negative multinomial regression models

[ripley@stats.ox.ac.uk]

A multinomial regression model is not a glm.  As it name implies, glm.nb 
uses generalized linear models.  If it had been as simple as you suggest, 
I am sure the authors of glm.nb would already have done it.

One clue is that multinomial logistic regression models are not 
implemented (nor can they be) by extending glm.  The best idea is to 
follow what was done there, namely to define a log-likelihood and maximize 
it.  That's actually not a bad way to fit negative binomial glms, also.

On Wed, 22 Jan 2003, Maurice Masliah wrote:

> Hello,
> 
> I've spent a lot of time during the past month trying to get negative
> multinomial regression models for clustered event counts as described in
> (Guang Guo. 1996. "Negative Multinomial Regression Models For Clustered
> Event Counts." Sociological Methodology 26: 113-132., abstract at
> http://depts.washington.edu/socmeth2/4abst96.htm) implemented in R. A
> FORTRAN version of the method described in the Guo paper is available at
> http://www.stat.unipg.it/stat/statlib/general/negmul.
> 
> I've been approaching this problem by modifying the loglik glm.nb function
> as such:
> 
> ## original loglik function in glm.nb
> ##    loglik <- function(n, th, mu, y) {
> ##        sum(lgamma(th + y) - lgamma(th) - lgamma(y + 1) + th * 
> ##            log(th) + y * log(mu + (y == 0)) - (th + y) * log(th + 
> ##            mu))
> ##   }
> 
>       loglik <- function(n, th, mu, y, mu.grouped, y.grouped) {	
> 	sum( th*log(th))*length(y.grouped) + 
> 		sum(y*log(mu+ (y == 0)))+ 
> 		sum(lgamma(y.grouped+ th) -lgamma(th) -
> (y.grouped+th)*log(mu.grouped+ th))
>         }
> 
> where y.grouped and mu.grouped are the sums of the clusters, as in
> 
>     y.grouped <- tapply(Y,group,sum)	
>     mu.grouped <- tapply(mu,group,sum)
> 
> where group represents the level used to identify the clusters.
> 
> I've made some other small changes to glm.nb on the calls to theta.ml such
> that 
> 
> ##original call to theta.ml
> ##    th <- as.vector(theta.ml(Y, mu, n, limit = control$maxit, 
> ##        trace = control$trace > 2))
> 
>     th <- as.vector(theta.ml(y.grouped, mu.grouped, length(y.grouped), limit
> = control$maxit, 
>         trace = control$trace > 2))
> 
> 
> At first I thought all I needed was to make some small changes and modify
> the loglik function in glm.nb but implementation now appears to be more
> complex than just that.
> 
> I'm looking for any help/advice in what I should do in order to implement
> negative multinomial regression models in R. What more needs to be done
> after modifying the loglik function in glm.nb?
> 
> Of course, if it has already been done that would be nice too (though I
> don't think so).
> 
> I have not been able to use the web interface to log on to the mailing lists
> (I keep getting an error message) and can only review the archives, so it
> would be much appreciated if you would cc any response to my email address
> below.
> 
> Sincerely,
> 
> Maurice Masliah
> Senior Researcher
> iTRANS Consulting Inc.
> 100 York Boulevard
> Suite 300, Richmond Hill
> L4B 1J8, Ontario
> CANADA
> Tel:905-882-4100 ext.5295
> Fax:905-882-1557
> e-mail: mmasliah at itransconsulting.com
> 
> 	[[alternate HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595