[R] multinomial logistic regression with equality constraints?

[Roger Levy]

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)