I am trying to fit a particular binomial glm to the following data.
testData <- structure(list(lb = c(0.102, 0.128, 0.161, 0.203, 0.256, 0.323,
0.406, 0.512, 0.645, 0.813, 1.02, 1.29, 1.63, 2.05, 2.58, 3.25,
4.1, 5.16, 6.5, 8.19, 10.3, 13, 16.4, 20.6, 26), TotalParticles = c(50612,
20541, 18851, 8058, 5606, 4123, 1995, 1234, 677, 381, 202, 111,
75, 42, 64, 118, 127, 83, 9, 2, 4, 7, 67, 8, 143), binsize = c(0.026,
0.033, 0.042, 0.053, 0.067, 0.083, 0.106, 0.133, 0.168, 0.207,
0.27, 0.34, 0.42, 0.53, 0.67, 0.85, 1.06, 1.34, 1.69, 2.11, 2.7,
3.4, 4.2, 5.4, 6), vol = c(309.76, 309.76, 309.76, 309.76, 309.76,
309.76, 309.76, 309.76, 309.76, 309.76, 309.76, 309.76, 309.76,
309.76, 309.76, 309.76, 309.76, 309.76, 309.76, 309.76, 309.76,
309.76, 309.76, 309.76, 309.76)), class = c("tbl_df", "tbl",
"data.frame"), row.names = c(NA, -25L))
I have no trouble fitting a poisson glm to the data.
fit_model <- function(df) glm(TotalParticles ~ log(lb), offset = log(binsize * vol), data = df, family = "poisson")
summary(fit_model(testData))
Call: glm(formula = TotalParticles ~ log(lb), family = "poisson", data = df, offset = log(binsize * vol)) Deviance Residuals: Min 1Q Median 3Q Max -42.592 -3.216 1.539 14.514 40.162 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.276493 0.013243 96.39 <2e-16 *** log(lb) -3.216526 0.006667 -482.44 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1152482.7 on 24 degrees of freedom Residual deviance: 6489.9 on 23 degrees of freedom AIC: 6677.5 Number of Fisher Scoring iterations: 5
I can also fit a quasipoisson model. Which gives me exactly the same coefficients, though slightly /higher/ p-values.
fit_quasi <- function(df) glm(TotalParticles ~ log(lb), offset = log(binsize * vol), data = df, family = "quasipoisson")
summary(fit_quasi(testData))
Call: glm(formula = TotalParticles ~ log(lb), family = "quasipoisson", data = df, offset = log(binsize * vol)) Deviance Residuals: Min 1Q Median 3Q Max -42.592 -3.216 1.539 14.514 40.162 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.2765 0.9662 1.321 0.199 log(lb) -3.2165 0.4864 -6.613 9.54e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 5322.668) Null deviance: 1152482.7 on 24 degrees of freedom Residual deviance: 6489.9 on 23 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5
However, when I try to fit a negative binomial model, I get this cryptic error message
fit_nb = function(df) MASS::glm.nb(TotalParticles ~ log(lb) + offset(log(vol * binsize)), data = df)
fit_nb(testData)
Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset, : NA/NaN/Inf in 'x'
traceback()
3: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset, family = fam, control = list(maxit = control$maxit, epsilon = control$epsilon, trace = control$trace > 1), intercept = attr(Terms, "intercept") > 0) 2: MASS::glm.nb(TotalParticles ~ log(lb) + offset(log(vol * binsize)), data = df) at #1 1: fit_nb(testData)
The problem seems related to the non-linear nature of the data. For instance, if I zero out the last five values, I can do the binomial regression, though I get a warning.
testDataWorks <- testData
testDataWorks[(nrow(testDataWorks)-4):nrow(testDataWorks),"TotalParticles"] <- 0
fit_nb(testDataWorks)
glm.fit: algorithm did not converge Call: MASS::glm.nb(formula = TotalParticles ~ log(lb) + offset(log(vol * binsize)), data = df, init.theta = 1.66759784, link = log) Coefficients: (Intercept) log(lb) 1.765 -2.892 Degrees of Freedom: 24 Total (i.e. Null); 23 Residual Null Deviance: 4954 Residual Deviance: 29.95 AIC: 315.8
Why exactly is my binomial regression failing? Is there a work-around that I can do to get it to run anyway? I get that these data aren't really particularly linear, but I would like to be able to run the model anyway for reasons (I have a bunch of similar data andthis mostly works for them).
I see that my problem is similar to this one R: NA/NaN/Inf in X error however, the consensus seemed to be that that user's model had way too many parameters, while in my case I have many fewer paramters to predict total particles.
This post https://stats.stackexchange.com/questions/52527/unable-to-fit-negative-binomial-regression-in-r-attempting-to-replicate-publish may also be relevant, though I'm not exactly sure how I would apply it in practice.
Thanks for any advice.