[R] Questing on fitting Baseline category Logit model

[Christofer Bogaso]

Dear all,

I am facing some problem with how to fit a "Baseline category Logit
model" with R. Basically I am considering famous  "Alligator" data as
discussed by Agresti. This data can also be found here:

https://onlinecourses.science.psu.edu/stat504/node/174
(there is also an accompanying R file, however the underlying R code
could not load the data properly!!!)

Below are the stuffs what I have done so far:

My_Data <- structure(list(Number = c(7L, 4L, 1L, 0L, 0L, 0L, 0L, 1L, 5L,
2L, 16L, 3L, 3L, 0L, 2L, 1L, 2L, 2L, 3L, 3L, 2L, 13L, 2L, 7L,
0L, 6L, 0L, 0L, 1L, 0L, 3L, 0L, 9L, 1L, 1L, 0L, 0L, 1L, 2L, 0L,
3L, 8L, 7L, 6L, 1L, 6L, 0L, 3L, 1L, 5L, 2L, 0L, 4L, 1L, 1L, 0L,
1L, 0L, 4L, 0L, 13L, 9L, 10L, 0L, 0L, 0L, 2L, 1L, 2L, 2L, 3L,
8L, 9L, 1L, 1L, 0L, 0L, 0L, 1L, 1L), 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", "Invertebrate", "Other", "Reptile"), 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"), 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"),
    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")), .Names =
c("Number",
"Food", "Size", "Sex", "Lake"), row.names = c(NA, 80L), class = "data.frame")


library(VGAM)
vglm(Food~Size+Sex+Lake, data = My_Data, fam=multinomial, weights = Number)



However I am getting following error:

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
2: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :
  96 elements replaced by 1.819e-12
3: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :
  96 elements replaced by 1.819e-12
4: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :
  96 elements replaced by 1.819e-12
5: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :
  96 elements replaced by 1.819e-12
6: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :
  96 elements replaced by 1.819e-12

Can somebody points me why I am getting this error?

Thanks for you help