glm.nb - theta, dispersion, and errors
On Sun, 21 Oct 2012, Eiko Fried wrote:
I am running 9 negative binomial regressions with count data. The nine models use 9 different dependent variables - items of a clinical screening instrument - and use the same set of 5 predictors. Goal is to find out whether these predictors have differential effects on the items. Due to various reasons, one being that I want to avoid overfitting models, I need to employ identical types of models for all 9 regressions. The problem is that some of my dependent variables are overdispersed, others are not (parameter estimates gained by using family=quasipoisson). dispersion p theta S1 1.084 0.102 30.0 S2 0.903 0.000 4125.0 PROBLEM S3 0.997 0.926 3754.0 PROBLEM S4 0.784 0.000 10283.0 PROBLEM S5 1.108 0.003 8.6 S6 1.010 0.772 1287.0 PROBLEM S7 1.228 0.001 1.7 S8 1.222 0.005 0.4 S9 2.120 0.283 1.2 So I thought using maximum likelihood Poisson models would provide wrong results, and that using negative binomial models would be the best way to go. 3 models cause the following problem (the ones that are not overdispersed):
m2.nb <- glm.nb(t0s2 ~ Sex + HisDep + FamHis + ZEFE + ZNeuro, data=data)
Error in while ((it <- it + 1) < limit && abs(del) > eps) { :
missing value where TRUE/FALSE needed
What happens is that, looking at theta iterations, values become too large, and thereby infinite, and thereby NaN. This is an example for model 2 above, with a theta of 4125 and a significantly smaller dispersion value than 1 of .903 (p<.000):
theta.ml: iter76 theta =1.34118e+16 theta.ml: iter77 theta =3.0232e+16 theta.ml: iter78 theta =-Inf theta.ml: iter79 theta =NaN
One symptom causes this problem: 1: In sqrt(1/i) : NaNs produced How would you recommend to deal with these problems? Are results reliable although these errors occur? Again, I have to use a model that fits all 9 models overall best - maybe I should default to maximum likelihood Poisson?
From the description of your response in earlier e-mails on this list, my
feeling is that this is caused by the type of response you have. It is not really a count response and might be better modeled by an ordered response model. Below I've simulated data from an ordered probit model and then pretended they are count data. The result are the symptoms you describe: underdispersion, non-convergence for theta in NB. Fitting an ordered probit model, avoids these problems for this data and (not surprisingly) reasonably recovers the true parameters. hth, Z ## random regressor and latent Gaussian variable set.seed(1) x <- runif(200) y0 <- 5 * x + rnorm(200) ## collapse to factor with 5 levels yf <- cut(y0, c(-Inf, 3, 4, 5, 6, Inf), labels = 0:4) ## pretend levels are count variables yc <- as.numeric(as.character(yf)) ## quasi-Poisson with underdispersion summary(glm(yc ~ x, family = quasipoisson)) ## NB with problems in theta estimation summary(glm.nb(yc ~ x)) ## ordered probit recovers true parameters summary(polr(yf ~ x, method = "probit", Hess = TRUE)) ## compare with the latent linear model summary(lm(y0 ~ x))
Thank you Torvon [[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.