# Fitting polytomous logistic regression with R

I'm running into troubles fitting a polytomous logistic regression model using grouped data. The data are of the form (dput at bottom):

> head(alligator)
lake  sex  size    food count
1 Hancock male small    fish     7
2 Hancock male small  invert     1
3 Hancock male small reptile     0
4 Hancock male small    bird     0
5 Hancock male small   other     5
6 Hancock male large    fish     4


And I've tried to fit the model with vglm() from package VGAM:

> result <- vglm(food~lake+size+sex, data=alligator, fam=multinomial, weights=count)
Error in if (max(abs(ycounts - round(ycounts))) > smallno) warning("converting 'ycounts' to integer in @loglikelihood") :
missing value where TRUE/FALSE needed
1: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) : 96 elements replaced by 1.819e-12  It was also suggested to look at mlogit() from package globaltest (on Bioconductor), but it does not appear to support grouped data. It obviously doesn't support the weights parameter, but I can't find where the equivalent parameter is documented: source("http://bioconductor.org/biocLite.R") biocLite("globaltest") result <- mlogit(food~lake+size+sex, weights=count, data=alligator) Error in mlogit(food ~ lake + size + sex, weights = count, data = alligator) : unused argument(s) (weights = count)  If anyone could put me down the right path, I'd appreciate it! > dput(alligator) structure(list(lake = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("George", "Hancock", "Oklawaha", "Trafford"), class = "factor"), sex = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("female", "male"), class = "factor"), size = structure(c(2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("large", "small"), class = "factor"), food = structure(c(2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L), .Label = c("bird", "fish", "invert", "other", "reptile"), class = "factor"), count = c(7L, 1L, 0L, 0L, 5L, 4L, 0L, 0L, 1L, 2L, 16L, 3L, 2L, 2L, 3L, 3L, 0L, 1L, 2L, 3L, 2L, 2L, 0L, 0L, 1L, 13L, 7L, 6L, 0L, 0L, 3L, 9L, 1L, 0L, 2L, 0L, 1L, 0L, 1L, 0L, 3L, 7L, 1L, 0L, 1L, 8L, 6L, 6L, 3L, 5L, 2L, 4L, 1L, 1L, 4L, 0L, 1L, 0L, 0L, 0L, 13L, 10L, 0L, 2L, 2L, 9L, 0L, 0L, 1L, 2L, 3L, 9L, 1L, 0L, 1L, 8L, 1L, 0L, 0L, 1L)), .Names = c("lake", "sex", "size", "food", "count"), class = "data.frame", row.names = c(NA, -80L))  ## 2 Answers Figured it out, I think the class notes need to be updated. This seems to work and matches the output from the course notes. First, reshape the date into wide format: require(reshape2) a2 <- dcast(lake + sex + size ~ food, data = alligator, value.var="count")  Then fit the model with a matrix of responses. The last level passed in becomes the reference level for the models that are fit (fish in this case): result2 <- vglm(cbind(bird, invert, other, reptile, fish) ~ lake + sex + size, data = a2, family = multinomial)  • Did you try bringing your Count variable in as an integer type rather than a long type? That was what the error was telling you. R looks at the variable type, not the actual variable data. Mar 10, 2012 at 18:39 • @Michelle - unless I'm missing something, count is an integer (as it is fed into vglm()). See > class(alligator$count) [1] "integer" Mar 10, 2012 at 18:43
• Hmm... did you try a traceback() when you got your error? Mar 10, 2012 at 19:14
• @Michelle. Good thought. I looked at traceback and it looks like it drops down into an internal function vglm.fitter(). I'll do some more digging and report back. Mar 10, 2012 at 22:09

According to the vglm documentation,

The values of weights must be positive; try setting a very small value such as 1.0e-8 to effectively delete an observation.

Meaning, they can't be zero as in your dataset. So if you'd do something like

alligator$$count <- pmax(alligator$$count, 1e-8)


before calling the vglm you should be good to go.