[R] multinomial logistic regression with equality constraints?

[Jasjeet Singh Sekhon]

As we noted earlier and as is clearly stated in the docs, multinomRob
is estimating an OVERDISPERSED multinomial model.  And in your models
here the overdispersion parameter is not identified; you need more
observations.  Walter pointed out using the print.level trick to get
the coefs for the standard MNL model, but when the model the function
is actually trying to estimate is not identified, that trick will not
work.

As I also previously noted, it is a simple matter to change the
multinomMLE function to estimate the standard MNL model.  Since you
don't want to change that file and since nnet's multinom function
doesn't have some functionality people need, we'll add a "MLEonly"
function to multinomRob which will allow you to do what you want.
We'll post a new version on my webpage later tonight:
http://sekhon.berkeley.edu/robust.  And after some testing, we'll
forward the new version to CRAN.

Jas.

=======================================
Jasjeet S. Sekhon                     
                                      
Associate Professor             
Travers Department of Political Science
Survey Research Center          
UC Berkeley                     

http://sekhon.berkeley.edu/
V: 510-642-9974  F: 617-507-5524
=======================================




Roger Levy writes:
 > Walter Mebane wrote:
 > > Roger,
 > > 
 > >  > Error in if (logliklambda < loglik) bvec <- blambda :
 > >  > 	missing value where TRUE/FALSE needed
 > >  > In addition: Warning message:
 > >  > NaNs produced in: sqrt(sigma2GN)
 > > 
 > > That message comes from the Newton algorithm (defined in source file
 > > multinomMLE.R).  It would be better if we bullet-proofed it a bit
 > > more.  The first thing is to check the data.  I don't have the
 > > multinomLogis() function, so I can't run your code.  
 > 
 > Whoops, sorry about that -- I'm putting revised code at the end of the 
 > message.
 > 
 > > But do you really
 > > mean
 > > 
 > >  > for(i in 1:length(choice)) {
 > > and
 > >  > dim(counts) <- c(length(choice),length(choice))
 > > 
 > > Should that be
 > > 
 > >   for(i in 1:n) {
 > > and
 > >   dim(counts) <- c(n, length(choice))
 > > 
 > > or instead of n, some number m > length(choice).  As it is it seems to
 > > me you have three observations for three categories, which isn't going
 > > to work (there are five coefficient parameters, plus sigma for the
 > > dispersion).
 > 
 > I really did mean for(i in 1:length(choice)) -- once again, the proper 
 > code is at the end of this message.
 > 
 > Also, I notice that I get the same error with another kind of data, 
 > which works for multinom from nnet:
 > 
 > 
 > library(nnet)
 > library(multinomRob)
 > dtf <- data.frame(y1=c(1,1),y2=c(2,1),y3=c(1,2),x=c(0,1))
 > summary(multinom(as.matrix(dtf[,1:3]) ~ x, data=dtf))
 > summary(multinomRob(list(y1 ~ 0, y2 ~ x, y3 ~ x), data=dtf,print.level=128))
 > 
 > 
 > The call to multinom fits the following coefficients:
 > 
 > Coefficients:
 >      (Intercept)          x
 > y2 0.6933809622 -0.6936052
 > y3 0.0001928603  0.6928327
 > 
 > but the call to multinomRob gives me the following error:
 > 
 > multinomRob(): Grouped MNL Estimation
 > [1] "multinomMLE: -loglik initial: 9.48247391895106"
 > Error in if (logliklambda < loglik) bvec <- blambda :
 > 	missing value where TRUE/FALSE needed
 > In addition: Warning message:
 > NaNs produced in: sqrt(sigma2GN)
 > 
 > 
 > Does this shed any light on things?
 > 
 > 
 > Thanks again,
 > 
 > Roger
 > 
 > 
 > 
 > 
 > 
 > ***
 > 
 > set.seed(10)
 > library(multinomRob)
 > multinomLogis <- function(vector) {
 >    x <- exp(vector)
 >    z <- sum(x)
 >    x/z
 > }
 > 
 > n <- 20
 > choice <- c("A","B","C")
 > intercepts <- c(0.5,0.3,0.2)
 > prime.strength <- rep(0.4,length(intercepts))
 > counts <- c()
 > for(i in 1:length(choice)) {
 >    u <- intercepts[1:length(choice)]
 >    u[i] <- u[i] + prime.strength[i]
 >    counts <- c(counts,rmultinomial(n = n, pr = multinomLogis(u)))
 > }
 > dim(counts) <- c(length(choice),length(choice))
 > counts <- t(counts)
 > row.names(counts) <- choice
 > colnames(counts) <- choice
 > data <- data.frame(Prev.Choice=choice,counts)
 > 
 > for(i in 1:length(choice)) {
 >    data[[paste("last",choice[i],sep=".")]] <- 
 > ifelse(data$Prev.Choice==choice[i],1,0)
 > }
 > 
 > multinomRob(list(A ~ last.A ,
 >                   B ~ last.B ,
 >                   C ~ last.C - 1 ,
 >                   ),
 >              data=data,
 >              print.level=128)
 > 
 > 
 > 
 > I obtained this output:
 > 
 > 
 > Your Model (xvec):
 >                                 A B C
 > (Intercept)/(Intercept)/last.C 1 1 1
 > last.A/last.B/NA               1 1 0
 > 
 > [1] "multinomRob:  WARNING.  Limited regressor variation..."
 > [1] "WARNING.  ... A regressor has a distinct value for only one 
 > observation."
 > [1] "WARNING.  ... I'm using a modified estimation algorithm (i.e., 
 > preventing LQD"
 > [1] "WARNING.  ... from modifying starting values for the affected 
 > parameters)."
 > [1] "WARNING.  ... Affected parameters are TRUE in the following table."
 > 
 >                                     A     B     C
 > (Intercept)/(Intercept)/last.C FALSE FALSE  TRUE
 > last.A/last.B/NA                TRUE  TRUE FALSE
 > 
 > 
 > 
 > multinomRob(): Grouped MNL Estimation
 > [1] "multinomMLE: -loglik initial: 70.2764843511374"
 > Error in if (logliklambda < loglik) bvec <- blambda :
 > 	missing value where TRUE/FALSE needed
 > In addition: Warning message:
 > NaNs produced in: sqrt(sigma2GN)