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@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. ______________________________________________ R-help@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.