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.