2
$\begingroup$

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
In addition: Warning messages:
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))
Improve this question
$\endgroup$

2 Answers 2

Reset to default
0
$\begingroup$

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)
Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Michelle
    Commented Mar 10, 2012 at 18:39
  • $\begingroup$ @Michelle - unless I'm missing something, count is an integer (as it is fed into vglm()). See > class(alligator$count) [1] "integer" $\endgroup$
    – Chase
    Commented Mar 10, 2012 at 18:43
  • $\begingroup$ Hmm... did you try a traceback() when you got your error? $\endgroup$
    – Michelle
    Commented Mar 10, 2012 at 19:14
  • $\begingroup$ @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. $\endgroup$
    – Chase
    Commented Mar 10, 2012 at 22:09
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.