Christofer Bogaso
2012-Mar-14 08:15 UTC
[R] Questing on fitting Baseline category Logit model
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
2012-Mar-14 15:19 UTC
[R] Questing on fitting Baseline category Logit model
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