[R] multinomial logistic regression with equality constraints?

[Roger Levy]

Many thanks for this, Jas.  I was successfully able to use the revised 
version of multinomRob, and it satisfies exactly the needs I was looking 
for.

Thanks once again.

Best,

Roger

Jasjeet Singh Sekhon wrote:
> 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)
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.