10

When is it appropriate to use a Type I versus Type II negative binomial distribution in a zero-inflated negative binomial distribution?

I've found a Similar question, but without an answer I can comprehend or determine if it relates to zero-inflated negative binomial distributions

In my dataset using R code, I have determined using AIC values that the zero-inflated negative binomial (ZINB) distribution provides the best fit compared with other distribution models.

Using the R package glmmTMBi have specified ZINB models with both Type I and Type II:

library(glmmTMB);library(MuMIn)
m1 <- glmmTMB(dv ~ iv1 + iv2 + iv2,ziformula=~.,data=df,family=nbinom1)
m2 <- glmmTMB(dv ~ iv1 + iv2 + iv2,ziformula=~.,data=df,family=nbinom2)

When testing the models' AIC values, the model with Type II provides a better fit

AICc(m1,m2)
   df     AICc
m1  9 528.3359
m2  9 527.7481

When using the zeroinfl function from the pscl

zm2 <-  zeroinfl(dv ~ iv1 + iv2 + iv2,data=df,dist="negbin")
> AICc(m2,zm2)
    df     AICc
m2   9 527.7481
zm2  9 527.7481

It yields the same AIC value as Type 2 NB (and the estimtes and p-values are nearly identical), So it seems that Type 2 is assumed in the zeroinfl function.

My data set models the use of drugs (as the dv) over the past 30 days.

Is it appropriate to use a Type II negative binomial distribution and why?. Is it reasonable to justify this decision with AIC values?

Nick Cox
  • 56,404
  • 8
  • 127
  • 185
Timo Kvamme
  • 347

3 Answers3

11

The difference between these two model families is the relationship between mean and variance.

nbinom1 (also called quasi-poisson) variance = µ * phi

where µ is the mean and phi is the over-dispersion parameter

nbinom2 (the default negative binomial in most packages) variance = µ(1+µ/k) also written µ + (µ^2)/k

where µ is the mean and k is the over-dispersion parameter

When choosing between these the paper by VerHoef, J.M. & Boveng is very helpful as are pages 16 and 17 of Bolker et al 2012.

VerHoef, J.M. & Boveng say that AIC doesn't necessarily apply to quasi poisson models (nbinom1) and they are skeptical about comparing AIC and qAIC (an information criteria developed for quasi models) although you do see it done.

Instead they recommend plotting the observed values against the squared residuals. This plot can be very noisy so grouping samples with similar observed values together and making the equivalent plot for the groups is recommended. If this plot follows a linear trend it suggests quasi-poisson (nbinom1) is best whereas a quadratic trend argues for a negative binomial model (nbinom2).

If you have a decent number of samples and a finite number of possible combinations of explanatory variables you could form groups not based on response variables but on treatment combinations. This plot is demonstrated in Bolker et al 2012 (link in the references) along with code to generate the plot in R.

Ben Bolker, Mollie Brooks, Beth Gardner, Cleridy Lennert, Mihoko Minami, October 23, 2012, Owls example: a zero-inflated, generalized linear mixed model for count data. https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf

VerHoef, J.M. & Boveng, P.L., 2007. Quasi-Poisson Vs. Negative Binomial Regression: How Should We Model Overdispersed Count Data? Ecology, 88(11), pp.2766–2772.

D A Wells
  • 385
  • 2
  • 12
3

Users of binomial and Poisson models are often surprised to learn that semiparametric models work extremely well for count data while making no assumptions about the distribution of counts for any one type of observation. The question about how one can model an entire distribution empirically, which may require thousands of intercepts when there are thousands of distinct $Y$ values, comes up frequently.

Consider the special case where there are no covariates (predictors), and no ties in the data ($n$ observations). The empirical cumulative distribution function (ECDF) estimates the $Y$ distribution unbiasedly, without assumptions. It is a nonparametric estimator with $n-1$ parameters. In ordinal semiparametric models such as the proportional odds model, the ECDF is just the estimated cumulative probabilities when there are no covariates. If there are $k$ distinct $Y$ values there are $k-1$ intercepts (here $n-1$ if there are no ties in the data). The intercepts are exactly the logits of the ECDF values. The effective number of parameters in this model is low because of the order restrictions in the intercepts, and because the ECDF has effectively about 3-5 degrees of freedom. This can been seen by comparing widths of confidence bands from the ECDF for the unknown cumulative distribution with widths from fitting 2-, 3-, 4-, and 5-parameter parametric distributions.

Now consider covariates. You can easily estimate more parameters than there are observations without problem, as intercepts are order-restricted so effectively only have a few degrees of freedom. Computationally thousands of intercepts are handled (in R rms::orm and SAS JMP only) by respecting the sparse nature of the information matrix, which is tri-band diagonal). Intercepts are like covariate-adjusted ECDF values. Interestingly, Bayesian ordinal models do not need to invert information matrices so they automatically handle lots of intercepts without special programming (see R rmsb::blrm).

Intercepts encode the whole ECDF so distribution assumptions are gone. Ordinal models can handle extreme ties in the data (the limiting case is the binary logistic model which is a special case of the proportional odds model) so they handle zero inflation as well as ceiling effects and upper limits of detection.

Frank Harrell
  • 91,879
  • 6
  • 178
  • 397
  • That's quite incredible. So for even continuous variables we can use ordinal models, making less restrictive distributional assumptions, with little loss in efficiency. That would solve a lot of modelling problems that I have worked on (cognitive measurements with ceilings). Do you prefer a particular link function for ordinal models? If so, why? – Lachlan Oct 01 '23 at 03:55
  • 1
    Yes I always use ordinal models now for continuous Y because of those major advantages. Otherwise there is too much uncertainty in how to transform continuous Y, and Y with heavy ties doesn’t work with standard models. I like the logit link because of interpretability and not having mathematical constraints. Choose between that one and one much different from it - the log-log link (Cox proportional hazards model). No reason to use probit which is too similar to logit but hard to interpret. For a case study looking a multiple link functions for continuous Y see hbiostat.org/rmsc – Frank Harrell Oct 01 '23 at 12:39
  • @FrankHarrell does the same property apply for the continuation ratio model or only for the proportional odds one? Also, what is the coefficients' interpretation for the proportional odds model (with a logit link) for a continuous outcome? Is it how covariates relate to the CDF of this outcome? – Dimitris Rizopoulos Oct 01 '23 at 19:39
  • I’m 2/3 sure it doesn’t apply to the CR model unless you transform it in to log-log link form (there is a paper about that equivalence). In its usual form the intercepts can bounce around. The easy properties I mention apply to any member of the cumulative probability class of semiparametric models. – Frank Harrell Oct 02 '23 at 12:10
1

It is appropriate to use both Type I and Type II negative binominal distributions, as their AIC values are very close to each other. Internally, they differ by different likelihood function construction, where each observation's probability mass formula is given by the assumed distribution. If one formula leads to a larger log likelihood value that corresponds to a smaller AIC, then it is more likely the data generation process followed this chosen formula.

Compared with choosing a family that minimizes AIC, more important questions are what I listed in this answer https://stats.stackexchange.com/a/627172/284766. A zero-inflated negative binominal model makes sense here if your participants include both nondrug users (who caused zero inflation) and drug users (who may use no drugs during your study window) IF you don't know who is drug user and who is not. In your case, modeling the dispersion is particularly important through glmmTMB(..., disp = ~ ...) if you want to stress that one drug use induces more drug use on the same person and examine how this self-reinforcement mechanism is related to any predictors, as positively correlated event occurrence is a common reason for overdispersion.

If you do know who are drug users and nonusers among your participants, which appears more plausible a case to me, you may want to use other modeling techniques, such as Heckman sample selection model that is very popular in criminology studies, in which it is not meaningful to discuss the count of drug usage among nonusers. You can present both zero-inflated count models and Heckman sample selection models and compare whether substantive conclusions differ between them.

DrJerryTAO
  • 1,514
  • 1
    I find it hard to make the assumptions in any of the binomial models when semiparametric models work well and don't make any distributional assumptions. – Frank Harrell Sep 24 '23 at 12:55
  • Hi @FrankHarrell, are you suggesting to use instead an ordinal regression to predict the frequency of drug use? It is a good choice if the observed integer values are very limited, such as only 0, 1, 2, 3 are recorded in the response variable. By the way, I came across your 1990 paper on partial proportional odds model, that I referred to in my dissertation before expanding the technique when specifying threshold-specific coefficients. There is still some distributional assumption in the error term in the latent-response formulation (e.g. logistic distribution in partial PO). – DrJerryTAO Sep 24 '23 at 13:09
  • Cool. Yes that is what I suggest. It's better to formulate it without latent variables as then the assumptions are easier to see. More about assumptions here. – Frank Harrell Sep 24 '23 at 14:32
  • I think the assumption is much easier to comprehend from a latent-response formulation, where the distributional assumption of errors gives the probability mass function and thus the likelihood function. I agree that ordinal regression can be a very good choice to predict drug use frequency, especially if it was measured in bins, often done in frequency questions, instead of all nonnegative integers. – DrJerryTAO Sep 26 '23 at 12:23
  • I find the opposite to be true, prefering to stick with observables. But I can see how others will disagree. – Frank Harrell Sep 26 '23 at 13:12
  • @FrankHarrell, it's amazing from your web page where it claims that lrm and orm can fit thousands of distinct values in the response. Is it done through penalization of coefficients or by binning intervals of response values into larger groups? An ordinal regression model estimates (k - 1) intercepts (thresholds) if there are k different levels in the response. If data record k unique response values across k observations, I fail to see how a model of (k - 1) intercepts can be any useful, as there will be no more observation left for estimating coefficients of any predictor. – DrJerryTAO Sep 30 '23 at 08:00
  • I've expanded on this in a new answer. – Frank Harrell Sep 30 '23 at 11:34