[R] Questing on fitting Baseline category Logit model

[Michael Friendly]

Not sure why VGAM::vglm doesn't work here, but most likely it is the 
small & zero counts cited on the page you quote below.  This data set
is very sparse.  You should communicate with the author of VGAM.

You can fit this model with nnet::multinom instead, something like

library(nnet)
# multinomial logit model
(mod1 <- multinom(food ~ lake+size+sex, data=Alligator, weights=count))

 > multinom(food ~ lake+size+sex, data=Alligator, weights=count)
# weights:  35 (24 variable)
initial  value 352.466903
iter  10 value 270.397070
iter  20 value 268.958046
final  value 268.932740
converged
Call:
multinom(formula = food ~ lake + size + sex, data = Alligator,
     weights = count)

Coefficients:
         (Intercept) lakeHancock lakeOklawaha lakeTrafford sizesmall 
  sexmale
fish     1.70178892  -0.5752403    0.5503569  -1.23679067 0.7303298 
0.60639521
invert   0.53452560  -2.3557451    1.4635491  -0.08096449 2.0665999 
0.14342792
other   -0.01957203   0.1913919    0.5764707   0.32102428 1.0209285 
0.35382356
reptile -1.15700455   0.5539263    3.0803954   1.82403973 0.1733300 
-0.02116283

Residual Deviance: 537.8655
AIC: 585.8655
 >


On 3/14/2012 4:15 AM, Christofer Bogaso wrote:
> 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
>


-- 
Michael Friendly     Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    Web:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA