 Research
 Open Access
 Published:
A comparison of statistical methods for modeling count data with an application to hospital length of stay
BMC Medical Research Methodology volume 22, Article number: 211 (2022)
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 postmedical 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 nonnegative values, typically right skewed, and often exhibiting excessive zeros. In this study, we compared the performance of the Poisson, negative binomial (NB), zeroinflated Poisson (ZIP), and zeroinflated 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 zeroinflation 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 zeroinflation 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 zeroinflation only and NB models for overdispersed data or data with combinations of zeroinflation 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 zeroinflation 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/COVID19 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 hospitalacquired 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 postmedical 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 sociodemographic 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 rightskewed 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 nontransformed or logarithmtransformed 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 logarithmtransformed 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 logtransformations 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 nonnormality 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 logodds 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), zeroinflated Poisson (ZIP), and zeroinflated 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 zeroinflation. 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 zeroinflated 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 zeroinflation 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 zeroinflation. 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 integervalued 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 kdimensional 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 Poissondistributed with:
and the mean parameter (i.e., the mean number of events per period) is given by:
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):
Poisson regression model is also called loglinear model because the logarithm of the conditional mean is linear in the parameters:
The marginal effect of a predictor variable \({X}_{j}\) is given by:
and the interpretation of this effect is that a oneunit 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)\)
Reallife count data often exhibit two (related) characteristics: overdispersion and zeroinflation. Overdispersion refers to an excess of variability in the data (i.e., the variance exceeds the mean), while zeroinflation 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:
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}\):
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:
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:
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.
Zeroinflated count models
Zeroinflated 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:
where \(\mathrm{g}({y}_{i})\) follows either Poisson or NB distributions, defined in (1) and (8), respectively, and therefore the zeroinflated count models are called either zeroinflated Poisson (ZIP) or zeroinflated 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 zeroinflated covariates and \(\gamma\) is a (q + 1)dimensional vector of zeroinflated regression coefficients to be estimated. The function F is called zeroinflated 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:
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:
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 nonnormal 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_{1} and x_{2} were generated from the following four distributions: Poisson, NB, ZIP, and ZINB. Variable x_{1} was continuous and generated from a normal distribution with mean µ = 57.3 and a variance σ^{2} = 306.25 representing the distribution of variable age observed in the Medical Information Mart for Intensive Care (MIMICIII) dataset for patients with an asthma diagnosis [54,55,56]. The binary variable x_{2} was generated from Bernoulli distribution with probability of success p = 0.43, representing the distribution of variable sex in the MIMICIII dataset for patients with an asthma diagnosis. The values of the population regression coefficient β_{0}, β_{1}, and β_{2} were prespecified and obtained by fitting a NB regression model for the outcome variable hospital LOS in the same MIMICIII 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.
Models evaluation
Poisson, NB, ZIP, and ZINB regression models with dependent variable y and independent variables x_{1} and x_{2} were fitted on the generated data under each of the simulation scenarios and replications using the maximum likelihood estimation (MLE) method [57]. The QuasiNewton 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 widelyused 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 (MIMICIII) [54,55,56]. MIMICIII is a large, singlecenter 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 (ICD9) 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 selfpay. 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 Chisquare 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 twosided 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 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 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.
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 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.
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.
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 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.
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.
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 sixtyfour 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%.
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.
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%.
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 ICD9 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].
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 zeroinflated 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 secondbest 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 (zeromodel), none of the variables in both ZIP and ZINB regression models had significant contribution to the structural zerogenerating process of LOS.
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 reallife 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 zeroinflation 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 zeroinflation 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 QuasiNewton 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 zeroinflation. Tlhaloganyang et al. [45] reported that the NB model provided a reasonable fit in all datasets when compared to ZIP and ZINB models in overdispersed 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 nonstructured 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 zeroinflated model and suggested the use of ZINB instead of ZIP when modeling zeroinflated count data. He stated that the zeroinflated 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 zeroinflated, 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 MIMICIII 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 selfpay 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 QuasiNewton or Newton–Raphson. In this study, we used QuasiNewton 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 zeroinflated data. Other optimization methods that are often used to fit zero inflated data (as in the case of R software) are the NelderMead 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, twopart 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 truncatedatzero 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 zerotruncated Poisson or NB regression models. Hospital LOS data are rarely zeroinflated. 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, nonspecific 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 rightcensored timetoevent 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 rightcensoring or disregarding the data for individuals who die prior the events. An alternative approach to analyze LOS timetoevent 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 betweenhospital variation in patient outcome; or LOS data are collected for the same subjects at multiple time points and there is a systematic betweensubject 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 withinclusters dependencies [80].
Concussions
The Poisson and ZIP regression models performed poorly in overdispersed data. ZIP outperformed the Poisson, NB, and ZINB regression models when there is just zeroinflation but no overdispersion in the data. NB model provided the best fit in overdispersed data and outperformed ZINB model in many cases of both zeroinflation 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 zeroinflated 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 riskadjusted LOS.
Availability of data and materials
The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.
Abbreviations
 AIC:

Akaike’s information criteria
 BIC:

Bayesian information criteria
 COVID19:

Corona virus disease 2019
 GEE:

Generalized estimating equation
 GLM:

Generalized linear model
 GLMM:

Generalized linear mixed model
 ICC:

Intraclass correlation coefficient
 ICD9:

International classification of diseases, ninth revision
 LOS:

Length of stay
 MAE:

Mean absolute error
 MIMIC:

Medical Information mart for Intensive care
 MLE:

Maximum likelihood estimation
 NB:

Negative binomial
 NIS:

Nationwide inpatient sample
 SAS:

Statistical analysis system
 ZINB:

Zeroinflated negative binomial
 ZIP:

Zeroinflated Poisson
References
Thomas JW, Guire KE, Horvat GG. Is patient length of stay related to quality of care? United States. 1997;42:489–507.
Taheri PA, Butz DA, Greenfield LJ. Length of stay has minimal impact on the cost of hospital admission. United States. 2000;191:123–30.
Kossovsky MP, Sarasin FP, Chopard P, LouisSimonet M, Sigaud P, Perneger TV, et al. Relationship between hospital length of stay and quality of care in patients with congestive heart failure. England. 2002;11:219–23.
Khalifa M. Reducing Length of Stay by Enhancing Patients’ Discharge: A Practical Approach to Improve Hospital Efficiency. Netherlands. 2017;238:157–60.
Baek H, Cho M, Kim S, Hwang H, Song M, Yoo S. Analysis of length of hospital stay using electronic health records: A statistical and data mining approach. PLoS One. 2018;13(4):e0195901. Available from: https://doi.org/10.1371/journal.pone.0195901.
Giraldi G, Montesano M, Sandorfi F, Iachini M, Orsi GB. Excess length of hospital stay due to healthcare acquired infections: methodologies evaluation. Italy. 2019;31:507–16.
Rees EM, Nightingale ES, Jafari Y, Waterlow NR, Clifford S, Pearson CAB, et al. COVID19 length of hospital stay: a systematic review and data synthesis. BMC Med. 2020;18:270. Available from: https://doi.org/10.1186/s12916020017263.
Systematic, datadriven approach lowers length of stay and improves care coordination [Internet]. 2018. Available from: https://www.healthcatalyst.com/success_stories/reducinglengthofstaymemorialhospitalatgulfport. Accessed 16 Feb 2022.
Freeman WJ, Weiss AJ, Heslin KC. Overview of U.S. Hospital Stays in 2016: Variation by Geographic Region. 2018. In: Healthcare Cost and Utilization Project (HCUP) Statistical Briefs [Internet]. Rockville: Agency for Healthcare Research and Quality (US); 2006. Statistical Brief #246.
Pickering BW, Dong Y, Ahmed A, Giri J, Kilickaya O, Gupta A, et al. The implementation of clinician designed, humancentered electronic medical record viewer in the intensive care unit: a pilot stepwedge cluster randomized trial. Int J Med Inform. 2015;84:299–307 Ireland: Elsevier Ireland Ltd.
Lingsma HF, Bottle A, Middleton S, Kievit J, Steyerberg EW, Marangvan de Mheen PJ. Evaluation of hospital outcomes: the relation between lengthofstay, readmission, and mortality in a large international administrative database. BMC Health Serv Res. 2018;18(1):116. Available from: https://doi.org/10.1186/s1291301829161.
BetancourtGarcia MM, Vatcheva K, Gupta PK, Martinez RD, McCormick JB, FisherHoch SP, et al. The effect of Hispanic ethnicity on surgical outcomes: An analysis of the NSQIP database. Am J Surg. 2019;217:618–33 United States: Elsevier Inc.
Almashrafi A, Elmontsri M, Aylin P. Systematic review of factors influencing length of stay in ICU after adult cardiac surgery. BMC Health Serv Res. 2016;16:318. Available from: https://doi.org/10.1186/s1291301615913.
Rosenthal MJ, Fajardo M, Gilmore S, Morley JE, Naliboff BD. Hospitalization and Mortality of Diabetes in Older Adults: A 3year prospective study. Diabetes Care. 1998;21:231–5. Available from: https://doi.org/10.2337/diacare.21.2.231.
Aro S, Kangas T, Reunanen A, Salinto M, Koivisto V. Hospital use among diabetic patients and the general population. United States. 1994;17:1320–9.
Bo S, Ciccone G, Grassi G, Gancia R, Rosato R, Merletti F, et al. Patients with type 2 diabetes had higher rates of hospitalization than the general population. United States. 2004;57:1196–201.
Carter EM, Potts HW. Predicting length of stay from an electronic patient record system: a primary total knee replacement example. BMC Med Inform Decis Mak. 2014;14:26. Available from: https://doi.org/10.1186/147269471426.
Comino EJ, Harris MF, Islam MD, Tran DT, Jalaludin B, Jorm L, Flack J, Haas M. Impact of diabetes on hospital admission and length of stay among a general population aged 45 year or more: a record linkage study. BMC Health Serv Res. 2015;15:12. Available from: https://doi.org/10.1186/s1291301406662.
Feng CX, Li L. Modeling Zero Inflation and Overdispersion in the Length of Hospital Stay for Patients with Ischaemic Heart Disease. In: Chen DG, Chen J, Lu X, Yi GY, Yu H, editors. Singapore: Springer Singapore; 2016. p. 35–53. Available from: https://doi.org/10.1007/9789811025945_3.
Cheng SW, Wang CY, Ko Y. Costs and Length of Stay of Hospitalizations due to DiabetesRelated Complications. J Diabetes Res. 2019;2019:2363292. Available from: https://doi.org/10.1155/2019/2363292.
Donnan PT, Leese GP, Morris AD. Diabetes Audit and Research in Tayside SMUC. Hospitalizations for people with type 1 and type 2 diabetes compared with the nondiabetic population of Tayside, Scotland: a retrospective cohort study of resource use. Diabetes Care. 2000;23:1774–9 United States.
Priyadi A, Permana H, Muhtadi A, Sumiwi SA, Sinuraya RK, Suwantika AA. CostEffectiveness Analysis of Type 2 Diabetes Mellitus (T2DM) Treatment in Patients with Complications of Kidney and Peripheral Vascular Diseases in Indonesia. Healthcare (Basel). 2021;9(2):211. Available from: https://doi.org/10.3390/healthcare9020211.
Dictionary SM. length of stay. (n.d.) [Internet]. 2011. Available from: https://www.who.int/data/gho/indicatormetadataregistry/imrdetails/2541. Accessed 16 Feb 2022.
Shaaban AN, Peleteiro B, Martins MRO. Statistical models for analyzing count data: predictors of length of stay among HIV patients in Portugal using a multilevel model. BMC Health Serv Res. 2021;21:372. Available from: https://doi.org/10.1186/s12913021063891.
World Health Organization. The Global Health Observatory. Explore a world of health data. Length of stay for inpatient shortterm treatment, days [Internet]. Available from: https://www.who.int/data/gho/indicatormetadataregistry/imrdetails/2541. Accessed 16 Feb 2022.
Bert F, Kakaa O, Corradi A, Mascaro A, Roggero S, Corsi D, Scarmozzino A, Siliquini R. Predicting Length of Stay and Discharge Destination for Surgical Patients: A Cohort Study. Int J Environ Res Public Health. 2020;17(24):9490. Available from: https://doi.org/10.3390/ijerph17249490.
Garrison SR, Schneider KE, Singh M, Pogodzinski J. Preoperative physical therapy results in shorter length of stay and discharge disposition following total knee arthroplasty: a retrospective study. J Rehabil Med Clin Commun. 2019;2:1000017. Available from: https://doi.org/10.2340/200307111000017.
Lim ATP. Methods for analyzing hospital length of stay with application to inpatients dying in Southern Thailand. Glob J Health Sci. 2009;1(1):27. Available from: https://doi.org/10.5539/gjhs.v1n1p27.
Mullahy J. Much ado about two: reconsidering retransformation and the twopart model in health econometrics. Netherlands. 1998;17:247–81.
Manning WG. The logged dependent variable, heteroscedasticity, and the retransformation problem. Netherlands. 1998;17:283–95.
Gardner W, Mulvey EP, Shaw EC. Regression analyses of counts and rates: Poisson, overdispersed Poisson, and negative binomial models. Psychol Bull. 1995;118:392–404 US: American Psychological Association.
O’Hara RB, Kotze DJ. Do not logtransform count data. Methods Ecol Evol. 2010;1:118–22. Available from: https://doi.org/10.1111/j.2041210X.2010.00021.x.
MartínFernández JA, Hron K, Templ M, Filzmoser P, PalareaAlbaladejo J. Modelbased replacement of rounded zeros in compositional data: Classical and robust approaches. Computational Statistics & Data Analysis. 2012;56:2688–704. Available from: https://www.sciencedirect.com/science/article/pii/S0167947312000941.
Bryk AS, Raudenbush SW, Congdon RT. HLM: hierarchical linear and nonlinear modeling with the HLM2L and HLM3L programs. Chicago: Scientific Software International; 1996.
Huang JQ, Hooper PM, Marrie TJ. Factors associated with length of stay in hospital for suspected communityacquired pneumonia. Egypt. 2006;13:317–24.
Sroka CJ, Nagaraja HN. Odds ratios from logistic, geometric, Poisson, and negative binomial regression models. BMC Med Res Methodol. 2018;18:112. Available from: https://doi.org/10.1186/s1287401805689.
Lambert D. ZeroInflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics. 1992;34:1–14. Available from: https://www.tandfonline.com/doi/abs/10.1080/00401706.1992.10485228.
Cameron AC, Trivedi P. Regression Analysis of Count Data, 2nd edition, 2013. Econometric Society Monograph No.53. Cambridge University Press; 1998.
Hilbe JM. Modeling Count Data [Internet]. Cambridge: Cambridge University Press; 2014. Available from: https://www.cambridge.org/core/books/modelingcountdata/BFEB3985905CA70523D9F98DA8E64D08.
Agresti A. Foundations of linear and generalized linear models. Wiley; 2015.
Greene W. Accounting for excess zeros and sample selection in poisson and negative binomial regression models [Internet]. In: Leonard N. Stern School of Business, Department of Economics. New York University; 1994.
Slymen DJ, Ayala GX, Arredondo EM, Elder JP. A demonstration of modeling count data with an application to physical activity. Epidemiol Perspect Innov. 2006;3:3. Available from: https://doi.org/10.1186/1742557333.
Lee JH, Han G, Fulp WJ, Giuliano AR. Analysis of overdispersed count data: application to the Human Papillomavirus Infection in Men (HIM) Study. Epidemiol Infect. 2012;140(6):108794. Available from: https://doi.org/10.1017/S095026881100166X.
Tüzen F, Erbaş S, Olmuş H. A simulation study for count data models under varying degrees of outliers and zeros. Communications in Statistics  Simulation and Computation. 2020;49:1078–88. Available from: https://doi.org/10.1080/03610918.2018.1498886.
Tlhaloganyang BP TK. Are zero inflated distributions compulsory in the presence of zeroinflation? Int J Innov Sci Res Tech. 2020;5:1274–7.
Poisson SD. Recherches sur la probabilité des jugements en matière criminelle et en matière civile, précédées des règles générales du calcul des probabilités. ImprimeurLibraire pour les Mathématiques, Paris: Bachelier; 1837.
Yang Z, Hardin JW, Addy CL. Score Tests for ZeroInflation in Overdispersed Count Data. Communications in Statistics  Theory and Methods. 2010;39:200830. Available from: https://doi.org/10.1080/03610920902948228.
Perga F Pierre de, and Apollonius, of. Varia Opera Mathematica. olosæ: apud Johannem Pech. 1679. Available from: https://doi.org/10.5479/sil.128299.39088002705879.
Greenwood OR, Yule GU. An Inquiry into the Nature of Frequency Distributions Representative of Multiple Happenings with Particular Reference to the Occurrence of Multiple Attacks of Disease or of Repeated Accidents. J R Stat Soc. 1920;83:255–79. Available from: https://doi.org/10.1111/j.23972335.1920.tb00606.x John Wiley & Sons, Ltd.
Eggenberger F, Pólya G. über die Statistik verketteter Vorgange. Z Angew Math Mech. 1923;3:27989. https://doi.org/10.1002/zamm.19230030407.
He H, Tang W, Wang W, CritsChristoph P. Structural zeroes and zeroinflated models. Shanghai Arch Psychiatry. 2014;26:236–42 China.
Nelder JA, Wedderburn RWM. Generalized Linear Models. J R Stat Soc Ser A. 1972;135:37084. Available from: https://doi.org/10.2307/2344614.
McCullagh PNJ. Generalized Linear Models. 2nd ed. London: Chapman and Hall; 1989.
Johnson AE, Pollard TJ, Shen L, Lehman LW, Feng M, Ghassemi M, et al. MIMICIII, a freely accessible critical care database. England. 2016;3:160035.
Johnson A, Pollard T, Mark R. MIMICIII Clinical Database Demo (version1.4). 2019. Available from: https://doi.org/10.13026/C2HM2Q.
Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PC, Mark RG, et al. PhysioBank, PhysioToolkit, and PhysioNet. American Heart Association. 2000;101:e215–20. Available from: https://doi.org/10.1161/01.CIR.101.23.e215.
Enders CK. Maximum Likelihood Estimation. In Encyclopedia of statistics in behavioral science. American Cancer Society; 2005. https://doi.org/10.1002/0470013192.bsa200.
Beaujean AA, Grant MB. Tutorial on using regression models with count outcomes using R. Practical Assessment, Research, and Evaluation, Vol. 21, Article 2. 2016. Available from: https://doi.org/10.7275/pj8ch254.
Akaike H. A New Look at the Statistical Model Identification. In: Parzen E, Tanabe K, Kitagawa G, editors. Selected Papers of Hirotugu Akaike. Springer Series in Statistics. 1974. Available from: https://doi.org/10.1007/9781461216940_16.
Schwarz G. Estimating the Dimension of a Model. Ann Statist. 1978;6(2):4614. Available from: https://doi.org/10.1214/aos/1176344136.
ZZeileis A, Kleiber C, Jackman S. Regression Models for Count Data in R Journal of Statistical Software. 2008;27(8):1–25. Available from: https://doi.org/10.18637/jss.v027.i08.
Nekesa F, Odhiambo C, Chaba L. Comparative assessment of zeroinflated models with application to HIV exposed infants data. Open J Stat. 2019;9:664–85. Available from: https://doi.org/10.4236/ojs.2019.96043.
Minami M, LennertCody CE, Gao W, RománVerdesoto M. Modeling shark bycatch: The zeroinflated negative binomial regression model with smoothing. Fish Res. 2007;84:210–21 Available from: (https://www.sciencedirect.com/science/article/pii/S0165783606003821).
Saffari SE, Adnan R, Greene. Handling of Overdispersion of Count Data via Truncation using Poisson Regression Model. Journal of Computer Science & Computational Mathematics. 2011;1(1). Available from: https://doi.org/10.20967/jcscm.2011.01.001.
Sawyer R. Sample Size and the Accuracy of Predictions Made from Multiple Regression Equations. Am Educ Res J. 1982;7:91–104. Available from: https://doi.org/10.3102/10769986007002091.
Tang W, Lu N, Chen T, Wang W, Gunzler DD, Han Y, et al. On performance of parametric and distributionfree models for zeroinflated and overdispersed count responses. Stat Med. 2015;34:3235–45 England: John Wiley & Sons, Ltd.
Brewer MJ, Butler A, Cooksley SL. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. Methods Ecol Evol. 2016;7:679–92. Available from: https://doi.org/10.1111/2041210X.12541.
Allison PD. Logistic Regression Using SAS®: Theory and Application. 2nd ed. Cary, NC: SAS Institute Inc.; 2012.
Soyiri IN, Reidpath DD, Sarran C. Asthma length of stay in hospitals in London 2001–2006: demographic, diagnostic and temporal factors. PLoS One. 2011;6(11):e27184. Available from: https://doi.org/10.1371/journal.pone.0027184.
Arora S, Kaur P, Panaich SS, Sagar H, Levine D. Asthma Exacerbations, Length of Stay and Hospitalization Costs: Insights from the Nationwide Inpatient Sample. J Allergy Clin Immunol. 2015;135(2):AB241.
Blumer A, Ehrenfeucht A, Haussler D, Warmuth MK. Occam’s Razor. Inf Process Lett. 1987;24:377–80. Available from: https://www.sciencedirect.com/science/article/pii/0020019087901141.
AlMahtot M, BarwiseMunro R, Wilson P, Turner S. Changing characteristics of hospital admissions but not the children admitteda whole population study between 2000 and 2013. Eur J Pediatr. 2018;177(3):381–8. Available from: https://doi.org/10.1007/s004310173064z.
Turner S, Raja EA. The association between opening a short stay pediatric assessment unit and trends in short stay hospital admissions. BMC Health Serv Res. 2021;21(1):523. Available from: https://doi.org/10.1186/s1291302106541x.
Thiruvengadam G, Lakshmi M, Ramanujam R. A Study of Factors Affecting the Length of Hospital Stay of COVID19 Patients by CoxProportional Hazard Model in a South Indian Tertiary Care Hospital. J Prim Care Community Health. 2021. Available from: https://doi.org/10.1177/21501327211000231.
Brock GN, Barnes C, Ramirez JA, Myers J. How to handle mortality when investigating length of hospital stay and time to clinical stability. BMC Med Res Methodol. 2011;11:144. Available from: https://doi.org/10.1186/1471228811144.
Taylor SL, Sen S, Greenhalgh DG, Lawless M, Curri T, Palmieri TL. A competing risk analysis for hospital length of stay in patients with burns. JAMA Surg. 2015;150(5):450–6. Available from: https://doi.org/10.1001/jamasurg.2014.3490.
Perez A, Chan W, Dennis RJ. Predicting the Length of Stay of Patients Admitted for Intensive Care Using a First Step Analysis. Health Serv Outcomes Res Methodol. 2006;6(3–4):127–38. Available from: https://doi.org/10.1007/s1074200600099.
Sotoodeh M, Ho JC. Improving length of stay prediction using a hidden Markov model. AMIA Jt Summits Transl Sci Proc. 2019;2019:425–34 Available from: (https://pubmed.ncbi.nlm.nih.gov/31258996).
Xie H, Chaussalet TJ, Millard PH. A continuous time Markov model for the length of stay of elderly people in institutional longterm care. J R Statist SocA. 2005;168(1):51–61. Available from: https://doi.org/10.1111/j.1467985X.2004.00335.x.
Zeger SL, Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics. 1986;42(1):12130.
Wolfinger R, O’connell M. Generalized linear mixed models a pseudolikelihood approach. J Stat Comput Simul. 1993;48:233–43. Available from: https://doi.org/10.1080/00949659308811554.
Goldstein H, Browne W, Rasbash J. Partitioning variation in multilevel models. Underst Stat. 2002;1:223–32. Available from: https://doi.org/10.1207/S15328031US0104_02.
Acknowledgements
Not applicable.
Funding
None.
Author information
Authors and Affiliations
Contributions
Research idea and study design: KPV, GAF; statistical analysis GAF, KPV; interpretation: GAF, KPV; wrote the manuscript: GAF, KPV; reviewed/edited manuscript: KPV. Each author: provided intellectual content; contributed significantly to the preparation and/or revision of the manuscript; and approved the final version of the manuscript. GAF takes responsibility for the integrity of the data and the accuracy of the data analysis. The authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
None.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Fernandez, G.A., Vatcheva, K.P. A comparison of statistical methods for modeling count data with an application to hospital length of stay. BMC Med Res Methodol 22, 211 (2022). https://doi.org/10.1186/s12874022016858
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12874022016858
Keywords
 Count data
 Poisson regression
 Negative binomial regression
 Zeroinflated Poisson regression
 Zeroinflated negative binomial regression
 Simulation study