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?

Thank you
Torvon

	[[alternative HTML version deleted]]

