A comparison of statistical methods for modeling count data with an application to hospital length of stay

Abstract

Background

Hospital length of stay (LOS) is a key indicator of hospital care management efficiency, cost of care, and hospital planning. Hospital LOS is often used as a measure of a post-medical procedure outcome, as a guide to the benefit of a treatment of interest, or as an important risk factor for adverse events. Therefore, understanding hospital LOS variability is always an important healthcare focus. Hospital LOS data can be treated as count data, with discrete and non-negative values, typically right skewed, and often exhibiting excessive zeros. In this study, we compared the performance of the Poisson, negative binomial (NB), zero-inflated Poisson (ZIP), and zero-inflated negative binomial (ZINB) regression models using simulated and empirical data.

Methods

Data were generated under different simulation scenarios with varying sample sizes, proportions of zeros, and levels of overdispersion. Analysis of hospital LOS was conducted using empirical data from the Medical Information Mart for Intensive Care database.

Results

Results showed that Poisson and ZIP models performed poorly in overdispersed data. ZIP outperformed the rest of the regression models when the overdispersion is due to zero-inflation only. NB and ZINB regression models faced substantial convergence issues when incorrectly used to model equidispersed data. NB model provided the best fit in overdispersed data and outperformed the ZINB model in many simulation scenarios with combinations of zero-inflation and overdispersion, regardless of the sample size. In the empirical data analysis, we demonstrated that fitting incorrect models to overdispersed data leaded to incorrect regression coefficients estimates and overstated significance of some of the predictors.

Conclusions

Based on this study, we recommend to the researchers that they consider the ZIP models for count data with zero-inflation only and NB models for overdispersed data or data with combinations of zero-inflation and overdispersion. If the researcher believes there are two different data generating mechanisms producing zeros, then the ZINB regression model may provide greater flexibility when modeling the zero-inflation and overdispersion.

Background

In healthcare, length of stay (LOS) is a key indicator used to assess the hospital care management efficiency, cost of care, quality control, appropriate use of hospital services and resources, and hospital planning [1,2,3,4,5,6]. The need for efficient hospital management has been exemplified with the recent onset of the 2019 coronavirus/COVID-19 pandemic. Health crises like these show the best interest of patients, hospitals, and public health is in the efficient management of hospital stays while ensuring adequate bed capacity and that clinician time can be provided for patients with other conditions [7]. Reducing LOS improves financial, operational, and clinical outcomes by decreasing the costs of care for a patient and minimizing the risk of hospital-acquired conditions [8, 9]. In some hospitals, administrators benefit from using predictive models to assist with planning and resource allocation for deliveries [9]. Clinics optimize clinical settings by implementing analytical applications leading to timely and accurate decision making while reducing the hospital LOS [8,9,10]. Hospital LOS is often used as a measure of a post-medical procedure outcome, as a guide to the benefit of a treatment of interest, and/or as an important risk factor for adverse events, hospital readmission, and mortality [11,12,13]. Therefore, understanding hospital LOS variability across various patients’ clinical and socio-demographic characteristics and hospitals’ characteristics, such as geographic region and hospital sizes, is always an important public health focus [9, 14,15,16,17,18,19,20,21,22].

Inpatient hospital LOS is the number of nights spent in hospital, calculated from the day of admission to the day of discharge [23]. This type of data can be treated as count data, and count data values are usually nonnegative with a typically right-skewed distribution, often exhibiting excessive zeros and overdispersion [17, 24, 25]. Different analytic strategies have been used for modeling hospital LOS. However, the best way to model LOS and other right skewed data has been debated in the literature. Literature review showed that non-transformed or logarithm-transformed count outcome variable are often modeled with linear regression [26,27,28]. Linear regression is usually employed for continuous, normally, or approximately normally distributed outcomes. LOS data rarely adheres to these assumptions. Studies conducted to compare analyses of logarithm-transformed count outcome variables have reported several issues that might arise with such transformations, including zero values not considered, predicted meaningless negative values for the outcome variable, uninterpretable and biased parameter estimates, and inconsistent inferences about important policy parameters [29, 30]. Gardner et al. [31] showed that when the mean of the count outcome variable is small, linear regression produces biased standard errors and hence biased significance tests. Using simulation study, O’Hara [32] found that the log-transformations of count data often used to satisfy parametric test assumptions perform poorly, except when the dispersion was small, and the mean counts were large. When the mean count is very small and zero is the most common value in the data set, the normalization with log transformation will not work and the mode will always be at the lowest value [33]. Bryk et al. [34] stated, that there are important cases for which the assumption of linearity and normality are not realistic, and no transformation can make them so. An alternative approach that has been used as a solution to handle the non-normality of LOS outcome variable by researchers is to dichotomize LOS and use logistic regression to predict the LOS [35]. Dichotomizing count outcome variable lead to loss of information. Based on simulated and empirical data analyses, Sroka [36] concluded that more precise odds ratios estimates can be obtained using count regression models with log-odds link function. In summary, using linear regression models with or without logarithmic transformation of a count outcome variable, or logistic regression models on a dichotomized count outcome variable are subject to criticism for their inadequacy in modeling this type of data. This can lead to biased parameter estimates; prediction of meaningless negative values; and the loss of precision of inferences and important information about the underlying counts.

Common statistical methods for analysis of count data are Poisson, negative binomial (NB), zero-inflated Poisson (ZIP), and zero-inflated negative binomial (ZINB) regressions [24, 37,38,39,40]. The results from the existing research evaluating the performance of regression models for count data are conflicting regarding which model is preferred. Lambert (1992) compared ZIP to NB regression models in an experimental study concerning soldering defects on printing wiring boards where 81% of the board areas had 0 defects. He found that ZIP was better than the NB model in terms of prediction accuracy [37]. Greene (1994) compared Poisson, NB, ZIP, and ZINB models on a consumer loan behavior empirical data characterized with overdispersion and zero-inflation. In the analysis the author found that the NB model was superior to the ZIP model and the ZIP model was superior to the Poisson model in terms of model fit [41]. Slymen et al. (2006) compared Poisson, overdispersed Poisson, NB, ZIP and ZINB regression models in assessing predictors of vigorous physical activity among Latina women using data with 82% zeros in the outcome variable. They reported a little difference in ZIP and ZINB models' fits, however, overall, the ZIP model fitted best [42]. In overdispersed and zero-inflated data of the number of incidents involving human papillomavirus infection, Lee et al. (2012) found that ZIP, followed by NB, and ZINB had the smallest Akaike’s information criteria (AIC); and ZIP model showed the same results as the NB model regarding the covariates at a 0. 05 significance level. In addition, ZINB model did not always converged [43]. Tuzen et al. (2018) examined the performance in terms of fit of Poisson, NB, ZIP, ZINB, Poisson Hurdle and NB Hurdle models under various outliers and zero-inflation scenarios of simulated data and found that ZINB and NB Hurdle were superior to Poisson, NB, and ZIP models. They also reported that in some scenarios, the NB model outperformed all models in the presence of outliers and/or excess zeros [44]. Tlhaloganyang et al. [45] compared NB with ZIP and ZINB models using different real datasets characterized by overdispersion and zero-inflation. The authors found that NB provided a superior fit in all datasets [45].

Based on the reviewed literature, the question remains open to whether the different results in terms of model fit may arise from the different proportion of zeros, overdispersion, and sample size of the datasets used in these studies. In this study we had two objectives. The first objective was to compare the performance of Poisson, NB, ZIP, and ZINB regression models in simulation study. The second objective was to compare the performance of Poisson, NB, ZIP, and ZINB regression models using real life hospital data in assessing the effect of age, sex, health insurance status, and type of hospital admission on the hospital LOS. This research added to previous studies by including additional experimental scenarios, such as varying sample sizes, larger dispersion levels, various proportions of zero in the outcome variable, and data generated using Poisson and ZIP distributions, along with NB and ZINB distributions.

Methods

Overview of count data regression models

Poisson model

The most widely used and the most basic model that explicitly considers the nonnegative integer-valued aspect of the count outcome variable is the Poisson regression model [46]. Let ${Y}_{i}, i=1,\dots ,n$, be random variables for the number of occurrences of the event of interest and its realizations ${y}_{i}=0, 1, 2\dots$. Let ${{\varvec{X}}}_{{\varvec{i}}}^{\boldsymbol{^{\prime}}}=\left({X}_{1i}, \dots , {X}_{ki}\right)$ be a k-dimensional random vector of predictors and its realization ${{\varvec{x}}}_{{\varvec{i}}}^{\boldsymbol{^{\prime}}}=\left({x}_{1i}, \dots , {x}_{ki}\right), i=1,\dots ,n$. Poisson regression assumes that the dependent variable Y_i, given ${{\varvec{X}}}_{{\varvec{i}}}={{\varvec{x}}}_{{\varvec{i}}}$ i = 1, …, n, is independently Poisson-distributed with:

$$P\left({{Y}_{i}=y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}\right)=\frac{{e}^{-{\mu }_{i}}{\mu }_{i}^{{y}_{i}}}{{y}_{i}!}, {y}_{i}=0, 1, 2, \dots$$

(1)

and the mean parameter (i.e., the mean number of events per period) is given by:

$${\mu }_{i}={e}^{{{\varvec{x}}}_{{\varvec{i}}}^{\boldsymbol{^{\prime}}}\beta }$$

(2)

where $\beta$ is a column vector of parameters.

In the Poisson regression model the conditional mean and the conditional variance of Y_i are equal (equidispersion):

$${E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)=V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)=\mu }_{i}$$

(3)

Poisson regression model is also called log-linear model because the logarithm of the conditional mean is linear in the parameters:

$${\mathrm{ln}(E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right))=\mathrm{ln}(\mu }_{i})={{\varvec{x}}}_{i}^{^{\prime}}\beta$$

(4)

The marginal effect of a predictor variable ${X}_{j}$ is given by:

$$\frac{\partial E\left({Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={{\varvec{x}}}_{{\varvec{i}}}\right)}{\partial {x}_{ji}}={{\beta }_{j}e}^{{{\varvec{x}}}_{{\varvec{i}}}^{\boldsymbol{^{\prime}}}\beta }={\beta }_{j}E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)$$

(5)

and the interpretation of this effect is that a one-unit change in the j^th predictor leads to a ${\beta }_{j}$ change in the conditional mean $E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)$

Real-life count data often exhibit two (related) characteristics: overdispersion and zero-inflation. Overdispersion refers to an excess of variability in the data (i.e., the variance exceeds the mean), while zero-inflation refers to an excess of zeros [39, 47]. In the presence of overdispersion, the Poisson regression model is not adequate and can lead to biased parameter estimates and unreliable standard errors estimates [38, 39]. The most commonly used model that accounts for overdispersion is the negative binomial model.

Negative binomial model

The Poisson regression model can be generalized by introducing an unobserved heterogeneity term for observation i. The subjects are assumed to differ randomly in a manner that is not fully accounted for by the observed covariates. This is formulated as:

$${E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i},{\tau }_{i}\right)=\mu }_{i}{\tau }_{i}={e}^{{{\varvec{x}}}_{{\varvec{i}}}^{\boldsymbol{^{\prime}}}\beta +{\varepsilon }_{i}}$$

(6)

where the unobserved heterogeneity term ${\tau }_{i}={e}^{{\varepsilon }_{i}}$ is independent of the vector of predictor variables x_i. Then the conditional distribution of Y_i on ${{\varvec{X}}}_{{\varvec{i}}}={{\varvec{x}}}_{{\varvec{i}}}$ is Poisson with conditional mean and conditional variance ${\mu }_{i}{\tau }_{i}$:

$$f\left({{Y}_{i}=y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}},{\tau }_{i}\right)=\frac{{e}^{-{\mu }_{i}{\tau }_{i}}{({\mu }_{i}{\tau }_{i})}^{{y}_{i}}}{{y}_{i}!}, {y}_{i}=0, 1, 2, \dots$$

(7)

The negative binomial distribution is derived as a gamma mixture of Poisson random variables [39, 48,49,50]. By letting $g\left({\tau }_{i}\right)$ be the probability density function of ${\tau }_{i}$, the distribution $f\left({{Y}_{i}=y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={{\varvec{x}}}_{{\varvec{i}}} \right)$ is obtained by integrating $f\left({{Y}_{i}=y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}},{\tau }_{i}\right)$ with respect to ${\tau }_{i}.$ The analytical solution of the integral exists if ${\tau }_{i}$ is gamma distributed and this solution is the NB distribution. Specifically, it is necessary to assume that $E\left({\tau }_{i}\right)=1,$ and then ${\tau }_{i}$ follows gamma (θ, θ) distribution with $E\left({\tau }_{i}\right)=1$ and $V\left({\tau }_{i}\right)=\frac{1}{\theta }.$ It can be shown, that the NB distribution can be written as:

$$f\left({{Y}_{i}=y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}\right)=\frac{\Gamma \left({y}_{i}+{\alpha }^{-1}\right)}{{y}_{i}!\Gamma \left({\alpha }^{-1}\right)}{\left(\frac{{\alpha }^{-1}}{{\alpha }^{-1}+{\mu }_{i}}\right)}^{{\alpha }^{-1}}{\left(\frac{{\upmu }_{\mathrm{i}}}{{\alpha }^{-1}+{\mu }_{i}}\right)}^{{\mathrm{y}}_{\mathrm{i}}}, {y}_{i}=0, 1, 2, \dots$$

(8)

where ${\alpha }^{-1}=\theta$ and $\theta >0$ is the gamma scale parameter.

The NB conditional mean and conditional variance of the outcome variable ${y}_{i}$ are given by:

$${E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)=\mu }_{i}$$

(9)

$$V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)={\mu }_{i}(1+\alpha {\mu }_{i})>E\left({{y}_{i}|x}_{i}\right)$$

(10)

The parameter $\alpha$ is defined as the dispersion parameter. As $\alpha$ approaches zero (i.e., the gamma scale parameter $\theta$ approaches infinity), $V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)$ decreases to ${\mu }_{i}$=$E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right),$ and the NB distribution approaches the Poisson distribution. Thus, the Poisson regression model is nested within the NB regression model.

Zero-inflated count models

Zero-inflated count models provide a way to both model the excess zeros and the overdispersion (He et al. 2014) [51]. In particular, there are two possible data generation processes for the number of occurrences of the event of interest ${y}_{i}$ for each observation i = 1,…,n and the result of a Bernoulli trial is used to determine which of the two to use. For observation i, Process 1 is chosen with probability ${\varphi }_{i}$ and Process 2 with probability $1-{\varphi }_{i}$. Process 1 generates only zero counts (“structural” zeros). Process 2 generates counts from either a Poisson model [37] or a NB model [41]. $P\left({Y}_{i}={y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}\right)$ can be described as follows:

$$P\left({Y}_{i}={y}_{i}|{{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}\right)=\left\{\begin{array}{ll}{\varphi }_{i}+(1-{\varphi }_{i})g(0)& {I}_{({y}_{i}=0)}\\ (1-{\varphi }_{i})\mathrm{g}({y}_{i})& {I}_{({y}_{i}>0)}\end{array}\right.$$

(11)

where $\mathrm{g}({y}_{i})$ follows either Poisson or NB distributions, defined in (1) and (8), respectively, and therefore the zero-inflated count models are called either zero-inflated Poisson (ZIP) or zero-inflated negative binomial (ZINB) regression models, respectively.

Further, if ${\varphi }_{i}$ depends on the characteristics of observation i, then ${\varphi }_{i}={F}_{i}=F({z}_{i}^{^{\prime}}\gamma )$, where ${z}_{i}$ is a (q + 1)-dimensional vector of zero-inflated covariates and $\gamma$ is a (q + 1)-dimensional vector of zero-inflated regression coefficients to be estimated. The function F is called zero-inflated link function.

In the case of the ZIP regression model, the conditional expectation and the conditional variance of the outcome variable ${Y}_{i}$ are given by:

$${E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)=\mu }_{i}(1-{F}_{i})$$

(12)

$$V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)=E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)(1+{F}_{i}{\mu }_{i})>E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)$$

(13)

Since $V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}, {z}_{i}\right)>E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}, {z}_{i}\right)$, ZIP model exhibits overdispersion as well.

In the case of ZINB regression model, the conditional expectation and the conditional variance of the outcome variable ${Y}_{i}$ are given by:

$${E\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)=\mu }_{i}(1-{F}_{i})$$

(14)

$$V\left({{Y}_{i}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{i}, {z}_{i}\right)=E\left({{{\varvec{Y}}}_{{\varvec{i}}}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}, {z}_{i}\right)[1+{(\alpha +F}_{i}{)\mu }_{i}]$$

(15)

Since $V\left({{{\varvec{Y}}}_{{\varvec{i}}}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}, {z}_{i}\right)>E\left({{{\varvec{Y}}}_{{\varvec{i}}}|{{\varvec{X}}}_{{\varvec{i}}}={\varvec{x}}}_{{\varvec{i}}}, {z}_{i}\right)$, the ZINB model like ZIP model exhibits overdispersion as well. Just as the NB distribution converges to the Poisson distribution as $\alpha$ approaches zero, the ZINB distribution converges to the ZIP distribution as $\alpha$ approaches zero.

Generalized linear models

Poisson, NB, ZIP, and ZINB are all part of the Generalized Linear Models (GLMs). The term GLM refers to a large class of models first introduced by Nelder and Wedderburn [52] and further developed and explained by McCullagh and Nelder [53]. GLMs extend standard linear regression models to encompass non-normal response distributions and possibly nonlinear functions of the mean [40]. The ordinary linear regression model uses linearity to describe the relationship between the mean of the response variable and a set of explanatory variables, with inference assuming that the response distribution is normal [40]. GLMs have three components: 1) A random component, that specifies the response variable Y_i, for the i^th observation and its probability distribution. 2) A inear component, ${\eta }_{i}={{\varvec{X}}}_{i}^{^{\prime}}\beta$, where $\beta$ is a column vector of parameters and ${{\varvec{X}}}_{i}$ is a column vector of predictors for the ith observation. 3) A monotonic differentiable link function g(.) describing how the expected value of variable Y_i is related to the linear predictor ${\eta }_{i}$, ${g\left[E\left(Yi\right)\right]=g({\mu }_{i}) = {\varvec{X}}}_{i}^{^{\prime}}\beta$, [40]. The response variable Y_i are independent for i = 1, 2, and have a probability distribution for an exponential family. This implies that the variance of the response variable Y_i depends on the mean ${\mu }_{i}$ through a variance function V: $var\left({Y}_{i}\right)=\frac{\phi V\left({\mu }_{i}\right)}{{\omega }_{i}},$ where $\phi$ is a constant, known as dispersion parameter, and ${\omega }_{i}$ is a known weight for each observation. The link function g for Poisson, NB, ZIP and ZINB regression models is log (${\eta }_{i}=\mathrm{log}({\mu }_{i})$). The binary link function h for the model of the probability of a zero count in the case of ZIP and ZINB regression models, is one of the logit, probit, or complementary log–log.

Simulation study

Dataset generation

Several datasets with one dependent variable y and two predictor variables x₁ and x₂ were generated from the following four distributions: Poisson, NB, ZIP, and ZINB. Variable x₁ was continuous and generated from a normal distribution with mean µ = 57.3 and a variance σ² = 306.25 representing the distribution of variable age observed in the Medical Information Mart for Intensive Care (MIMIC-III) dataset for patients with an asthma diagnosis [54,55,56]. The binary variable x₂ was generated from Bernoulli distribution with probability of success p = 0.43, representing the distribution of variable sex in the MIMIC-III dataset for patients with an asthma diagnosis. The values of the population regression coefficient β₀, β₁, and β₂ were pre-specified and obtained by fitting a NB regression model for the outcome variable hospital LOS in the same MIMIC-III dataset. For each of the simulated data under Poisson, NB, ZIP, and ZINB distributions, four different sample size scenarios were considered (50, 200, 600, and 1000). In the cases of count data generated with NB distribution or ZINB distribution different levels of dispersion (0.01, 1, 5, and 10) were considered under each of the sample size simulation scenarios. In the cases of count data generated with ZIP distribution or ZINB distribution, different proportions of structural zero (0.1, 0.3, 0.5, and 0.7) were considered under each of the sample size simulation scenario and under each of the dispersion levels simulation scenarios for data generated under ZINB distribution. To minimize the impact of simulation error, each scenario was repeated 1000 times. A summary of the simulation scenarios considered in the study is shown in Table 1.

Table 1 Simulation scenarios considered in the simulation study

Models evaluation

Poisson, NB, ZIP, and ZINB regression models with dependent variable y and independent variables x₁ and x₂ were fitted on the generated data under each of the simulation scenarios and replications using the maximum likelihood estimation (MLE) method [57]. The Quasi-Newton optimization technique was used to maximize the likelihood functions to obtain the regression models' estimates. To assess the performance of the four different models under each simulation scenario, we first calculated the models’ convergence rates. When the MLE procedure converged, it means it found a unique set of values for each parameter, the combination of which returned the highest likelihood value of all parameter values examined [58]. For the converged models, we extracted the widely-used Akaike’s Information Criteria (AIC), a model selection criterion developed by Hirotsugu Akaike [59].We also used the Bayesian information criteria (BIC) (also Schwarz criterion, SBC, SBIC) first formulated by Gideon Schwarz [60]. Smaller values of these criteria indicate a better model fit. In addition to the AIC and BIC statistics, we calculated the mean absolute error (MAE) for the E(y_i|x_i), defined as $MAE=\frac{{\sum }_{i=1}^{n}|E\left({y}_{i}|{x}_{i}\right)-E({\widehat{y}}_{i}|{x}_{i})|}{n},$ where n is the sample size. The AIC, BIC, and MAE were averaged over the 1000 replications in each of the simulation scenario. All simulations and statistical analyses were conducted using Statistical Analysis System (SAS) 9.4 (SAS Institute, Inc., 2015).

Empirical study

Data description

In the empirical study we used data from the Medical Information Mart for Intensive Care (MIMIC-III) [54,55,56]. MIMIC-III is a large, single-center database comprising information relating to patients admitted at the Beth Israel Deaconess Medical Center in Boston, Massachusetts. Agreement for data use was obtained. For the purposes of our study, we extracted patients with International Classification of Diseases, ninth revision (ICD-9) Code 49,390, which is related to the diagnosis of asthma. The resultant dataset consisted of 2,195 hospital patient admission records.

In our study, the outcome variable of interest was hospital LOS, calculated as the difference in days between the date of admission and the date of discharge. The predictor variables considered in the regression analyses were age, sex, patient health insurance and type of admission. Age was measured in years. Sex was a categorical variable with two levels: male and female. Patient health insurance was a categorical variable with 5 levels: government, Medicaid, Medicare, private, and self-pay. Admission type was a categorical variable with 3 levels: elective, emergency, and urgent.

Statistical analysis

First, we conducted descriptive statistical analysis to summarize and describe the study data. Frequencies and percentages were used to describe categorical variables and means, and standard deviations were used to describe continuous variables. The distribution of the count outcome variable LOS was visually examined using histogram. In addition, we calculated the variance and the mean of the outcome variable LOS to highlight potential Poisson distribution violations and overdispersion in the data. Poisson, NB, ZIP, and ZINB regression models were fitted for LOS on the predictor variables age, sex, health insurance, and admission type. The Pearson dispersion statistic, calculated by dividing the model's Pearson Chi-square statistic by the corresponding degrees of freedom, was used as a criterion for assessing model's misspecification or an overdispersed response variable. When the resultant value is greater than one, the model is considered to be overdispersed. AIC and BIC were used to compare the models. Furthermore, the models estimated coefficients, standard errors and their significance where examined, giving special attention to the difference in findings and conclusions across the models. All regression models were fitted using PROC COUNTREG and PROC GENMOD in SAS 9.4 (SAS Institute, 2015). Statistical testing was two-sided and performed at a significance (α) level of 0.05.

Results

Simulation study

Data generated with poisson regression model

Table 2 shows the convergence rates of Poisson, NB, ZIP, and ZINB regression models for simulated data in four different sample size scenarios (n = 50, 200, 600, 1000). NB regression convergence rate was between 53.5% and 58.7% with the largest convergence rate achieved in simulated samples of size 1000 (Table 2). The convergence rate for ZINB increased from 90.4% to 94.2% as the sample size increased from 50 to 1000. Both Poisson and ZIP regression models had 100% convergence rate in all simulation scenarios.

Table 2 Convergence rates of regression models fitted on data generated with Poisson regression model

Table 3 shows the averaged AIC and BIC statistics produced for Poisson, NB, ZIP, and ZINB regression models across all replications under each of the four sample size simulation scenarios. Poisson regression model (true model) resulted with the lowest mean AIC and BIC values, followed by the NB regression model. The difference in mean AIC and mean BIC values between the fitted Poisson and the fitted NB and ZINB regression models increased as the sample size increased.

Table 3 AIC and BIC of regression models fitted on data generated with Poisson regression model

Table 4 shows the MAE values of the predicted counts by the fitted Poisson, NB, ZIP, and ZINB regression models across all replications under each of the four sample size simulation scenarios. Poisson regression model (true model) resulted with the lowest MAE followed by the ZIP, ZINB and NB regression models.

Table 4 MAE values of regression models fitted on data generated with Poisson regression model

Data generated with NB regression model

In this section we describe the analysis of data generated with NB regression model in sixteen different simulation scenarios with varying sample sizes (n = 50, 200, 600, 1000) and magnitudes of dispersion (0.01, 1, 5, 10). Table 5 shows the convergence rates of Poisson, NB, ZIP, and ZINB regression models. At very low levels of overdispersion (0.01) the NB regression model's (true model) convergence rate ranged from 64.8% in the scenario with the smallest sample size to 94% in the scenario with the largest sample size. Similarly, in the same scenario of low level of overdispersion the ZINB regression model did not achieve 100% convergence rate. Overall, the ZINB regression model’s convergence rate was slightly better than the convergence rate of the NB model, ranging between 91% and 99.6% with the increase of the sample size. Poisson and ZIP regression models converged 100% in all simulation scenarios. In the rest of the scenarios of level of dispersion (1, 5, and 10) and all the sample sizes (n = 50,200,600 and 1000), all the models achieved 100% convergence rate.

Table 5 Convergence rates of regression models fitted on data generated with NB regression model

Table 6 displays the averaged AIC and BIC model fit statistics for Poisson, NB, ZIP, and ZINB regression models fitted on data generated with NB regression model with different magnitudes of dispersion and sample sizes. In the simulation scenario with nearly nonexistent overdispersion level (0.01), the Poisson regression model had the lowest AIC and BIC values regardless of the sample size. In all other simulation scenarios, the NB regression model (true model) had the lowest AIC and BIC values. The model fit of the NB and ZINB models improved as the level of dispersion increased; conversely, the model fit of Poisson and ZIP regression modes decreased as the dispersion level increased.

Table 6 AIC and BIC of regression models fitted on data generated with NB regression model

MAE values of the predicted counts with Poisson, NB, ZIP, and ZINB regression models fitted on data generated with NB regression model with different magnitude of dispersion and sample sizes are shown in Table 7. In the simulation scenario with the smallest sample size (n = 50), the Poisson regression model resulted with the smallest MAE regardless of the dispersion level. When the sample size increased to (n = 200) and at very small levels of overdispersion (0.01), both Poisson and ZIP regression models produced the smallest MAE, followed by NB regression model. When the dispersion level was 10 and sample seizes were greater than 200, NB, followed by ZIP regression model produced the lowest MAE. In the scenarios with sample sizes greater than 200 the true model (i.e., NB regression model) produced the smallest MAE regardless of the level of dispersion.

Table 7 MAE values of regression models fitted on data generated with NB regression model

Data generated with ZIP regression model

Table 8 shows the convergence rates of Poisson, NB, ZIP, and ZINB regression models fitted on data generated with a ZIP distribution in simulation scenarios with different levels of structural zeros and sample sizes. Poisson and ZIP regression models achieved 100% convergence rate across all the simulation scenarios. NB regression model achieved 98.4% and 99.9% convergence rate in data simulated with 10% proportion of structural zeros and sample sizes n = 50 and n = 200, respectively, and achieved 100% convergence rate in all the combination of scenarios for sample sizes and proportion of structural zeros greater than 30%. ZINB regression model convergence rate varied between 94.4% to 99.4% across all the simulation scenarios.

Table 8 Convergence rates of regression models fitted on data generated with ZIP regression model

Table 9 shows the AIC and BIC fit statistics values for the fitted Poisson, NB, ZIP, and ZINB regression models on data generated with ZIP regression model. When the proportion of structural zeros was 10%, the ZINB regression model had the smallest AIC and BIC values rather than the true ZIP regression model. In simulation scenarios with a higher proportion of structural zeros in the data (30%, 50%, 70%), the true ZIP regression model had the smallest AIC and BIC values followed by the ZINB regression model. In addition, the Poisson regression model produced the largest AIC and BIC values in all the scenarios. Another finding was that Poisson models fit became worse as the proportion of structural zeros increased from 10 to 50%, then slightly improved when the proportion of structural zeros reached 70%; contrary to the rest of the models, where the fit considerably became better as the proportion of structural zeros increased.

Table 9 AIC and BIC of regression models fitted on data generated with ZIP regression model

The MAE values of the predicted counts based on Poisson, NB, ZIP, and ZINB regression models fitted on data generated with the ZIP regression model in different simulation scenarios are shown in Table 10. ZIP model had the lowest MAE in the scenarios with the smallest sample size (n = 50) and proportion of structural zeros of 10% and 30%. Both ZIP and ZINB had the lowest MAEs in the scenario with the smallest sample size (n = 50) and 50% proportion of structural zeros. ZINB had the lowest MAE in the scenarios with sample sizes 200 and 1000 and 70% proportion of structural zeros. The ZIP model produced the lowest MAEs in scenarios with sample sizes 200 and 1000 and 10% proportion of structural zeros. When the sample size was 600 and the proportion of structural zeros was 30%, the ZINB had the highest MAE while Poisson, NB and ZIP had the same MAE. In all other scenarios, NB produced the lowest MAE.

Table 10 MAE values of regression models fitted on data generated with ZIP regression model

Data generated with ZINB regression model

In this section, we present the results from the analysis of the data generated with the ZINB regression model on sixty-four different simulation scenarios with various proportions of zeros, magnitudes of dispersion, and sample sizes. Tables 11, 12, 13, 14 present the convergence rates of Poisson, NB, ZIP, and ZINB regression models fitted on the simulated data across different simulation scenarios. In the scenario where data were simulated with the smallest sample size (n = 50), the Poisson model achieved 100% convergence rate regardless of the magnitudes of dispersion or structural zero proportions (Table 11). The convergence rates for the other models were unstable across the simulation scenarios and varied between 94.6% and 100%. For instance, ZIP model convergence rate was 100% for dispersion levels 0.01, 1, and 5 regardless of the proportion of structural zeros; and varied between 98.9% and 99.9% in the scenario with the largest overdispersion (dispersion = 10). When the sample size was 50, the convergence rate of ZINB slightly reduced as the proportion of zeros increased. Only in the scenarios with dispersion levels of 1 and 5 and proportion of structural zeros of (10%, 30%, and 50%), the ZINB model achieved 100% convergence rate. However, the trend observed in the convergence rate of ZINB in scenarios with sample size 50 changed as the sample size became larger (Tables 11, 12, 13, 14). With the increase of the sample size all the models achieved 100% convergence rate in all simulation scenarios except for the NB and ZINB regression models, where in simulation scenarios with the smallest level of dispersion (dispersion = 0.01) the models’ convergence rates were slightly below 100%.

Table 11 Convergence rates of regression models fitted on data generated with ZINB regression model, n = 50

Table 12 Convergence rates of regression models fitted on data generated with ZINB regression model, n = 200

Table 13 Convergence rates of regression models fitted on data generated with ZINB regression model, n = 600

Table 14 Convergence rates of regression models fitted on data generated with ZINB regression model, n = 1000

The AIC and BIC fit statistics of Poisson, NB, ZIP, and ZINB regression models fitted on data generated with the ZINB regression model and with different magnitudes of dispersion, zero proportions, and sample sizes are displayed in Tables 15, 16, 17, 18. In the simulation scenarios with the smallest dispersion level of 0.01 and 10% structural zeros, ZINB produced the lowest AIC and BIC fit statistics values regardless of the sample size, except for the case with the smallest sample size, where NB model had the lowest AIC values. However, ZINB had the lowest BIC values across all scenarios. When the proportion of structural zeros was greater than 10% ZIP model produced the lowest AIC and BIC statistics in all simulation scenarios with a dispersion level of 0.01 regardless of the sample size. When the dispersion level reached 1, the regression model that produced the lowest AIC values was the ZINB model in nearly all the scenarios of different sample sizes and proportions of structural zeros, except for the scenario with sample size 50, where NB produced the lowest AIC values in scenarios with a proportion of structural zeros up to 50%. ZINB produced the lowest AIC in the simulation scenarios with proportion of structural zeros greater than 50% and dispersion equal to 1. Similarly, at level of dispersion equal to 1 the model that produced the lowest BIC was the NB model in all the simulation scenarios of proportion of structural zeros and a small sample size 50. When the sample size increased from 50 to 1000, the NB regression model produced the lowest BIC statistics in the simulation scenarios with proportion of structural zero below 50%. The ZINB regression model produced the lowest BIC fit statistics in the scenarios with a proportion of structural zeros exceeding 50%. In summary, based on AIC and BIC statistics NB regression model consistently produced the best fit in all simulation scenarios with a dispersion level exceeding 1 regardless of the proportion of structural zeros and sample size.

Table 15 AIC and BIC of regression models fitted on data generated with ZINB regression model, n = 50

Table 16 AIC and BIC of regression models fitted on data generated with ZINB regression model, n = 200

Table 17 AIC and BIC of regression models fitted on data generated with ZINB regression model, n = 600

Table 18 AIC and BIC of regression models fitted on data generated with ZINB regression model, n = 1000

The MAEs values of the predicted counts based on Poisson, NB, ZIP, and ZINB regression models fitted on data generated with the ZINB regression model and with different magnitudes of dispersion, zero proportions, and sample sizes are displayed in Tables 19, 20, 21, 22. In scenarios with a small sample size of 50 and dispersion levels of 0.01 and 1, the ZIP model provided the smallest MAEs in nearly all simulation scenarios of proportion of structural zero, except for the scenarios with structural zero proportion of 70% and dispersion level of 0.01. Another exception were the scenarios with structural zero proportion of 10% and dispersion level of 1, where ZINB and Poisson regression models had the smallest MAEs, respectively. In the rest of the simulation scenarios within the small sample size of 50, the Poisson model provided the lowest MAEs (Table 19). As the sample size increased from 200 to 1000, the Poisson regression model no longer resulted with the smallest MAEs. In simulation scenarios with sample sizes 200, 600, and 1000, and a small dispersion of 0.01, the NB produced the lowest MAEs in nearly every scenario of proportion of structural zeros, with the exception when the proportion of structural zeros was 10% where the ZINB in (n = 200), ZIP in (n = 600) or either NB or ZIP in (n = 1000) regression models resulted with the lowest MAEs. Another finding was that in simulation scenarios with large sample size (n = 1000), and dispersion level of 1, the NB model produced the lowest MAEs regardless of the proportion of structural zeros. In simulation scenarios with large sample size (n = 1000) and dispersion level greater than 1, ZIP model produced the lowest MAEs when the proportion of structural zeros did not exceed 30%; and ZINB produced the lowest MAEs when the proportion of structural zeros was greater than 30%.

Table 19 MAE values of regression models fitted on data generated with ZINB regression model, n = 50

Table 20 MAE values of regression models fitted on data generated with ZINB regression model, n = 200

Table 21 MAE values of regression models fitted on data generated with ZINB regression model, n = 600

Table 22 MAE values of regression models fitted on data generated with ZINB regression model, n = 1000

Empirical study

Description of the study population

The empirical study population consisted of 2,167 patients admitted in hospitals with a diagnosis of Asthma selected from MIMIC dataset using ICD-9 code 49,390. Table 23 presents the main demographic characteristics of the study population. Sixty percent of the admitted patients were females. The mean age was 62.3 (SD = 40.66). It should be noted that in the MIMIC data, patients of under89 years old were merged into the same age group. 80.66% of the study population were admitted due to emergency, 17.44% were electively admitted, and just 1.89% of the patients had an urgent type of admission. Most patients had either Medicare (44.35%) or private health insurance (36.41%) (Table 23). The distribution of the variable hospital LOS was positively skewed, with values ranging from 0 to 40 days (Fig. 1). The mean LOS, 8.0 days, was much lower than the variance of 43.10. The larger sample variance compared to the sample mean suggested a deviation from the Poisson regression model’s assumption for equal variance and mean [61].

Table 23 Demographic and clinical characteristics of the study population with asthma diagnosis, n = 2,167

Fig. 1

Histogram of hospital length of stay for patients with asthma diagnosis, n = 2,167

Comparison of fitted poisson, NB, ZIP, and ZINB regression models

Table 24 presents the results of fitted Poisson, NB, ZIP, and ZINB regression models for the outcome variable LOS on the patient level predictor variables age, sex, type of hospital admission, and health insurance status. In the zero-inflated models the same predictors were used to fit both the count model and the logistic (zero) model. Based on the results in Table 24, the NB regression model provided the best fit to the data since it resulted the smallest AIC and BIC values. The second-best model was ZINB, followed by the ZIP model. The Poisson regression model resulted with the worst fit to the data according to the AIC and BIC values. The Pearson dispersion statistic in Poisson regression model was 5.3016, greater than 1, suggesting overdispersion. The fitted NB regression model had the smallest dispersion statistic of 1.1815. The regression coefficient estimates and their respective standard errors differed across the models (Table 24). It is quite noticeable in Table 24 the tendency for the Poisson, and ZIP regression models to produced smaller standard errors of the regression coefficient estimates than NB and ZINB regression models. Overdispersion may cause standard errors of the regression coefficient estimates to be underestimated and therefore contributing to discrepancies in significant regression coefficients findings between the models [39, 43]. For instance, at a 5% significance level, only based on the fitted Poisson and ZIP regression models there were significant association between age and log LOS, controlling for the effect of sex, health insurance type, and admission type variables included in the models (Table 24). In relation to the logistic part (zero-model), none of the variables in both ZIP and ZINB regression models had significant contribution to the structural zero-generating process of LOS.

Table 24 Findings from fitted multivariable Poisson, NB, ZIP, and ZINB regression models for hospital LOS, n = 2,167

Discussion

Simulation and empirical studies were conducted to compare performance in terms of convergence rate and model fit of the Poisson, NB, ZIP, and ZINB regression models. This research added to previous studies by including additional experimental scenarios, such as varying sample sizes, larger dispersion levels, various proportions of zero in the outcome variable, and data generated using Poisson and ZIP distributions, along with NB and ZINB distributions. Our motivating real-life example was the analysis of the count outcome variable hospital length of stay.

Based on the simulation study, when the data were generated with a Poisson regression model (i.e., there was no overdispersion or zero-inflation present in the data), the results showed that regardless of the sample size of the simulated data, the Poisson and ZIP regression models did not have convergence problems. Both NB and ZINB model did not converge 100% in all the sample size simulation scenarios. Compared to NB, ZIP, and ZINB regression models, the Poisson regression model (true model) resulted with the smallest AIC, BIC, and MAE. Our findings slightly differ from the findings reported by Nekesa et al. (2019) where in simulated data with fixed sample size of 500 with no zero-inflation and very low levels of overdispersion, the fitted NB model had the lowest AIC. However, in Nekesa’s study the response variable was generated with a negative binomial distribution [62]. By allowing the variance in the data to be greater than the mean, we generated overdispersed data with a NB regression model varying the level of dispersion and the sample size. We found that when the data have very low overdispersion, the Poisson regression model provided the smallest AIC and BIC statistics values regardless the sample size and fitting a NB or ZINB model may encounter convergence problems especially in situations where the data have low overdispersion. This was expected since as the dispersion parameter approaches zero, the NB distribution approaches the Poisson distribution [63]. When the dispersion was greater than 1, the NB model (true model) provided the best fit in terms of AIC and BIC, regardless of the sample size. These results are in line with Gardner [31] and Saffari [64], who showed that NB regression should be used when there is overdispersion in the data. In the scenarios with large sample sizes (200, 600 or 1000) the NB model produced the smallest MAE regardless of the dispersion level present in the data. Interestingly, when the sample size was 50, the Poisson regression model produced the smallest MAE. The reason for this could be that the small sample size affected the accuracy of the estimation [65]. When manipulating the sample size and the proportion of structural zeros in the data generated with a ZIP regression model, we found that ZINB had less than 100% convergence rate regardless of the sample size or proportion of structural zeros in the data. However, the NB model convergence rate was better than the ZINB regression model convergence rate since NB reached 100% in most of the simulation scenarios. It is important to note that the ZINB regression model is more complex than NB, which may influence the performance of the Quasi-Newton algorithm used for MLE estimations [66]. Other researchers have reported similar convergence issues in fitting the ZINB regression model. Lee et al. (2012) reported that the model did not always converge, or a model diagnostic indicated that the estimated model was not reliable. Another finding was that in scenarios with proportion of zeros of 30% or greater, the ZIP regression model (true model) had the best fit since it had the smallest AIC and BIC statistic values regardless of the sample size. This is in line with the findings reported by Nekesa et al. (2019), where in simulated conditions of very small overdispersion and proportion of zeros equal or greater than 20% ZIP model had smaller AIC than Poisson, NB and ZINB regression models [62]. An interesting finding of our study was that in the simulation scenario with 10% proportion of structural zeros, ZINB regression model had the best fit both in terms of AIC and BIC and regardless of the sample size. With respect to MAE, the performance of the models was very similar. In addition, our findings indicated that in scenarios with the smallest sample size the ZIP regression model had smallest MAE in most of the proportion of zeros scenarios (10%, 30%, 50%). However, when the sample size increased to 1000, the NB had the smallest MAE in most scenarios of proportion of zeros (30%, 50%, 70%). When varying the sample sizes, proportion of structural zeros and level of dispersion in the data generated with ZINB regression model, the ZINB model had convergence issues in the scenarios with small sample sizes. We found that at very low level of dispersion, 0.01, the model that produced the best fit both in terms of AIC and BIC was the ZIP regression model, regardless of the sample size, or proportion of structural zeros. This is not surprising, since, just as the NB distribution converges to the Poisson distribution as the dispersion parameter approaches zero, the ZINB distribution converges to the ZIP distribution as the dispersion parameter approaches zero [63]. In simulation scenarios with sample size greater than 50 and dispersion level fixed at 1, AIC suggested that the best model was ZINB (true model) regardless the proportion of structural zeros in the data. Research based on simulation studies on the use of AIC and BIC for model selection reported that BIC performed better in model selection in the case of large heterogeneity in data due to stronger penalty afforded [67]. In our study, according to the smallest BIC, the NB regression model fit the data better than ZINB in many of the scenarios, depending on the proportion of structural zeros and levels of overdispersion in the data. For instance, when the proportion of structural zeros exceeded 50% and the sample size was greater than 50, and with larger levels of overdispersion both BIC and AIC suggested that NB regression fits the data better. Similar findings were observed in a recent study conducted by Tlhaloganyang et al. [45] that showed that zero inflated models are not always necessary even if the data are characterized by both overdispersion and zero-inflation. Tlhaloganyang et al. [45] reported that the NB model provided a reasonable fit in all datasets when compared to ZIP and ZINB models in over-dispersed and zero inflated data. Similarly, Nekesa et al. (2019) reported results from simulation study and real data analysis of exposed infant diagnosis, showing the negative binomial emerging as the best performing model when fitting data with both structured and non-structured zeros under various settings. Tüzen et al. (2018) reported simulation scenarios, where the NB model outperformed other count models in the presence of outliers and/or excess zeros. Allison (2012) noted that some applications exist in which a compelling case could be made for a zero-inflated model and suggested the use of ZINB instead of ZIP when modeling zero-inflated count data. He stated that the zero-inflated negative binomial model may sometimes fit better than the conventional negative binomial, but for many applications it does not [68]. He recommended that in these cases, it’s important to test for the significance of the difference [68]. Lastly, Hilbe (2014) suggested that the model having substantially lower information test statistic should be preferred, other considerations being equal. In the situations when there is just a slight difference between which models fit statistics, the decision of which model to select should be based on context and how the models are to be interpreted [39].

The results from the empirical data analysis agreed with the findings based on the simulation study. The empirical data were not zero-inflated, and the data had overdispersion based on the Pearson dispersion statistic. The fitted NB regression model had the smallest AIC and BIC values followed by the ZINB regression model. This is in line with the findings of our simulation study, where NB was found to be the best model when dealing with overdispersed data. The Poisson and ZIP models underestimated the standard errors and overstated the significance of some covariates.

Since this study focused on regression analysis methods for count data, in this paper we do not fully discus the findings from the analysis of MIMIC-III data. However, the results from the NB regression analysis of the empirical data demonstrated, that health insurance type and admission type were significantly associated with the log transformed hospital LOS. Patients with elective admission had lower expected number of days of hospital stay compared to patients with urgent admission; and patients with Medicaid, Medicare, or private health insurance had longer expected number of days of hospital stay compared to self-pay patients, controlling for the effect of age and sex. Based on NB regression analysis, Soyiri et al. [69] reported significant associations of sex, age, admission type, ethnicity, week day of admission with asthma LOS in hospitals in London. Based on survey linear regression analysis of asthma Nationwide Inpatient Sample (NIS), collected between 2001 and 2010, Arora et al. [70] reported that white race and private insurance were significantly associated with longer and shorter LOS, respectively. However, the results across these studies are not fully comparable due to difference in variables’ definitions, covariates included in the regression models, type of regression methods, or potential difference in hospitals managements.

Different optimization procedures can produce different results and different rates of convergence. There is no perfect optimization procedure that finds the best solution within the most reasonable amount of time for all sets of data (SAS Institute, 2000). For the SAS programming language, the count regression models would typically be analyzed using PROC COUNTREG, PROC GENMOD, or PROC NLMIXED. The default optimization procedure here is usually the Quasi-Newton or Newton–Raphson. In this study, we used Quasi-Newton method that use iterative approximation and does not require computation of second order derivatives. Hence, it has the advantage of finding solutions quickly. However, this method does not consider the boundary constraints present in the zero-inflated data. Other optimization methods that are often used to fit zero inflated data (as in the case of R software) are the Nelder-Mead Simplex Optimization for ZIP regression models. The use of different optimization methods in the regression models across different software packages may explain some of the differences across studies conducted to evaluate regression methods.

AIC and BIC were used to determine the best model fitted to the data. When comparing the models, it is important to note that the best model will not be necessarily the one with the best fit. Rather, it will be the one that leads to correct inferences, interpretations, and decisions. Although one may not always know the exact model specification that will result in enhanced statistical conclusion, it is still possible to maintain a core principle that the ideal model should be simple and parsimonious [71].

Models’ performance measures were assessed based on the entire sample, but not based on models’ internal and/or external validation. The goal of our study was not to derive predictive modeling function for hospital LOS. The purpose of our study was to illustrate the choice of the count data regression model based on varying combinations of magnitudes of overdispersion, proportions of zeros, and sample sizes. In this research, we did not explore other count data regression models, such as Poisson and zero inflated Poisson inverse Gaussian, two-part Hurdle models, zero truncated model, mixture of a binomial and discretized gamma/beta distributions analysis, and others. Hurdle model is a modified count model in which the two processes generating the zeros and the positives are not constrained to be the same. The basic idea is that a binomial probability governs the binary outcome of whether a count variate has a zero or a positive realization. If the realization is positive, the “hurdle is crossed,” and the conditional distribution of the positives is governed by a truncated-at-zero count data model [38]. For example, the Hurdle model may be appropriate to analyze hospital length of stay if the data consisted of patients who were not hospitalized (i.e., zero days of hospital stay) and patients who were hospitalized. In this case, the probability of being hospitalized will be predicted by the logistic regression model and expected LOS will be predicted by the zero-truncated Poisson or NB regression models. Hospital LOS data are rarely zero-inflated. However, studies have reported rising in zero days hospital admissions (i.e., hospital stays of less than 24 h) in pediatric patients and in admissions with URTI/viral infection, gastroenteritis, croup, bronchiolitis, asthma, tonsillitis, non-specific abdominal pain, constipation, febrile convulsion, and rash diagnoses [69, 72, 73]. This may reflect a combination of factors including availability of more rapid assessment and effective treatment of acute presentation and declining hospital expertise and resources [72]. Further research should be conducted for scenarios of different data generating mechanisms in inpatient hospital LOS. Also, in this research we did not explore underdispersion. Even though is not as common to find underdispersed data in real life datasets, it would be interesting to evaluate the performance of count regression models in modeling such distributions. LOS can be analyzed as right-censored time-to-event data using survival analysis methods, where the event of interest is time to hospital discharge, or time to clinical stability, or time to death [74, 75]. If the interest is in estimating the probability of a patient reaching clinical stability or hospital discharge by a given day, Brock, et al. (2011) argued about right-censoring or disregarding the data for individuals who die prior the events. An alternative approach to analyze LOS time-to-event data with multiple events is to treat the events as a competing risk [75, 76]. Competing risk analyses extend survival analysis methods to situations with multiple possible events, where the occurrence of one either precludes the others or substantially alters the probability of other events [75, 76]. If additional measures such as, vital signs to monitor patients during hospitalizations or possible destinations after the first day of admission are obtained, Markov models can be utilized to capture the temporal sequences of events [77,78,79]. In this study we did not evaluate the extensions of the class of GLMs for analysis of count correlated data collected from multiple observations on individuals or count data that are clustered due to clustered survey data, family studies, or nested experimental designs. The generalized estimating equations (GEEs) [80] and the generalized linear mixed models (GLMMs) [81], also known as random effects models, multilevel, or hierarchical models, are used to deal with such clustered data and produce accurate regression coefficients and standard errors estimates. The technique of multilevel modeling incorporates cluster specific random effects to accounts for this within cluster dependency by partitioning the total data variance into between and within cluster variation [82]. For example, multilevel modeling will be appropriate if LOS data are collected from various participating hospitals with available hospital level variables and there is a systematic between-hospital variation in patient outcome; or LOS data are collected for the same subjects at multiple time points and there is a systematic between-subject variation in patient outcome. Whether a multilevel model is needed can be inferred from either the intraclass correlation coefficient (ICC) and/or the significance of the random effect variance component for the clustering variable in the null model. GEE methods account for correlation by incorporating predefined “working” correlation structures to describe the nature of within-clusters dependencies [80].

Concussions

The Poisson and ZIP regression models performed poorly in over-dispersed data. ZIP outperformed the Poisson, NB, and ZINB regression models when there is just zero-inflation but no overdispersion in the data. NB model provided the best fit in over-dispersed data and outperformed ZINB model in many cases of both zero-inflation and overdispersion. Just a slight difference existed between the fit statistics of NB and the more complex to fit and interpret ZINB model. The researcher should decide if a zero-inflated regression model is more appropriate to model the data. If the researcher believes there are two different data generating mechanism producing zeros, then the NB regression model may not capture the different characteristics of the two groups generating the zeros and, in this case, the ZINB regression model could provide greater flexibility when modeling the zeros. In addition, NB and ZINB regression models faced substantial convergence issues when incorrectly used to model equidispersed data. It is important to check for ovedispersion. Fitting incorrect models to overdispersed data leaded to incorrect regression coefficients estimates and overstated significance of some of the predictors.

Although the work presented here is based on the analysis of hospital LOS, the findings from the simulation study are generalizable to other count outcome variables. Our findings can guide in the selection from the studied generalized linear models in the development of hospital and public health analytical applications for the computation of risk-adjusted LOS.