A review of statistical estimators for risk-adjusted length of stay: analysis of the Australian and new Zealand intensive care adult patient data-base, 2008–2009
© Moran and Solomon; licensee BioMed Central Ltd. 2012
Received: 14 October 2011
Accepted: 16 April 2012
Published: 16 May 2012
For the analysis of length-of-stay (LOS) data, which is characteristically right-skewed, a number of statistical estimators have been proposed as alternatives to the traditional ordinary least squares (OLS) regression with log dependent variable.
Using a cohort of patients identified in the Australian and New Zealand Intensive Care Society Adult Patient Database, 2008–2009, 12 different methods were used for estimation of intensive care (ICU) length of stay. These encompassed risk-adjusted regression analysis of firstly: log LOS using OLS, linear mixed model [LMM], treatment effects, skew-normal and skew-t models; and secondly: unmodified (raw) LOS via OLS, generalised linear models [GLMs] with log-link and 4 different distributions [Poisson, gamma, negative binomial and inverse-Gaussian], extended estimating equations [EEE] and a finite mixture model including a gamma distribution. A fixed covariate list and ICU-site clustering with robust variance were utilised for model fitting with split-sample determination (80%) and validation (20%) data sets, and model simulation was undertaken to establish over-fitting (Copas test). Indices of model specification using Bayesian information criterion [BIC: lower values preferred] and residual analysis as well as predictive performance (R2, concordance correlation coefficient (CCC), mean absolute error [MAE]) were established for each estimator.
The data-set consisted of 111663 patients from 131 ICUs; with mean(SD) age 60.6(18.8) years, 43.0% were female, 40.7% were mechanically ventilated and ICU mortality was 7.8%. ICU length-of-stay was 3.4(5.1) (median 1.8, range (0.17-60)) days and demonstrated marked kurtosis and right skew (29.4 and 4.4 respectively). BIC showed considerable spread, from a maximum of 509801 (OLS-raw scale) to a minimum of 210286 (LMM). R2 ranged from 0.22 (LMM) to 0.17 and the CCC from 0.334 (LMM) to 0.149, with MAE 2.2-2.4. Superior residual behaviour was established for the log-scale estimators. There was a general tendency for over-prediction (negative residuals) and for over-fitting, the exception being the GLM negative binomial estimator. The mean-variance function was best approximated by a quadratic function, consistent with log-scale estimation; the link function was estimated (EEE) as 0.152(0.019, 0.285), consistent with a fractional-root function.
For ICU length of stay, log-scale estimation, in particular the LMM, appeared to be the most consistently performing estimator(s). Neither the GLM variants nor the skew-regression estimators dominated.
Length of stay during an intensive care unit (ICU) or hospital admission is a function of diverse patient and organisational input variables . It is widely used as an indicator of performance  and is a determinant of costs, although resource allocation is also known to affect length of stay . Not surprisingly, ICU length of stay has been the subject of frequent analysis [4–9], with the majority of studies presenting cross-sectional analyses over a relatively short periods of months  to 1–2 years .
ICU patient length of stay (and costs) demonstrate skewed distribution and various statistical modelling strategies have been employed in analysis of such data [11–14]; albeit linear regression (ordinary least squares regression, OLS) of the logged dependent variable has demonstrated a remarkable persistence . Individual patient data, as accessed from ICU data-bases, have an intrinsic hierarchical structure (patients within ICUs) and due analytic consideration of this structure is also appropriate . Using such data from the Australian and New Zealand Intensive Care (ANZICS) adult patient database (APD) , calendar years 2008–2009, the purpose of this paper was to: (i) compare the predictive performance and model specification of conventional estimators for skewed data (ICU length-of-stay); OLS (with both raw and log-scaled dependent variable) and generalised linear models (GLMs ), with more innovative approaches: multilevel or hierarchical linear mixed models (LMM) incorporating random effects [11, 16]; extended generalised linear models (EEE) with flexible link and variance functions ; estimators utilising skew-normal and skew-t multivariate distributions ; and finite mixture (FMM) models which consider the dependent variable as a mixture of distributions [20–22]; and (ii) determine the effect of (ICU) mortality outcome upon length of stay, allowing for “endogenous variable bias”  by means of a treatment-effects model [24, 25] where specific allowance is made for the correlation of independent predictors and error terms. Under such conditions, model coefficients are biased .
The ANZICS adult patient database (APD) is a bi-national (Australia and New Zealand) voluntary data collection of individual ICU admissions which commenced in 1990. The data set requirements are specified in a data dictionary . The current data-base does not incorporate coronary care units. Data is collected at the individual ICU and uploaded to the central repository (ANZICS APD) for processing and quality assurance; which consists of a cycle of error and exception checks, site feed-back, resubmission and incorporation into a final reporting data-set.
The ANZICS APD was interrogated to define an appropriate patient set, over the time period 2008–2009. Physiological variables collected were the worst in the first 24 hours after ICU admission. All first ICU admissions to a particular hospital for the period 2008–2009 were selected. Exclusions: patients with an ICU length of stay ≤ 4 hours; and patients aged < 16 years of age. Specific attention was directed to the fidelity of severity of illness records; in particular the scoring of the Glasgow Coma Score (GCS). Records were used only when all three components of the GCS were provided. Records where all physiological variables were missing were excluded and for the remaining records, missing variables were replaced to the normal range and weighted accordingly . ICU length of stay, initially recorded in hours was transformed to days; patients with an ICU length of stay > 60 days were not considered, no formal trimming methods were employed .
Continuous variables were reported using mean (SD), except where otherwise indicated; categorical variables were analysed using the Chi-squared test. To avoid the confounding effect of calendar time, patient data over the two years 2008 and 2009 was pooled. Distributional form of continuous variables of interest was by means of histogram, hanging rootogram [30, 31] or violin plots . The “hangroot”, an alternative to the histogram, compares an empirical distribution with a theoretical distribution by means of “hanging” spikes from the theoretical distribution. The spike lengths represent bin counts, using the square root of the frequencies to stabilise the sampling variation across bins and make deviations in the tails, where the counts are small, more visible. Deviations are shown as deviations from a horizontal line (y = 0) instead of deviations from a curve (the density function) in order to facilitate identification of patterns of deviation. The violin plot is a modification of the box plot that adds plots of the estimated kernel density to the summary statistics displayed by box plots , incorporates a marker for the median of the data, a box indicating the interquartile range, and spikes extending to the upper- and lower-adjacent values; overlaid is a density plot. Kernel density estimates are modifications of the histogram (a “smoothed” histogram), where densities are the continuous analogues of proportions (derivatives of the cumulative distribution function, so that areas under the density function read off as probabilities) . Kurtosis (a measure of the heaviness of the tails of the distribution) and skewness were quantitatively expressed (a normal distribution having a skewness of 0 and a kurtosis of 3) according to standard formulae .
Continuous: age (in years), severity of illness score (Acute Physiology and Chronic Health Evaluation (APACHE) III score ); both centred for computation. The predictive effect of these variables was entered initially as both linear and quadratic.
Categorical: these were parameterized as indicator variables with the reference level (≡ 0) indicated in parentheses in the following list: gender (female); mechanical ventilation (not ventilated); ICU level, as defined in the ANZICS database, as Rural/Regional, Metropolitan, Tertiary and Private (Tertiary); State of origin; that is New Zealand and the States of the Commonwealth of Australia (New South Wales (NSW), the largest contributor); patient surgical status as post-elective surgery, post-emergency surgery and non-surgical (non-surgical); descriptors of ICU admission primary organ system dysfunction, these being a consolidation of the “diagnostic categories” of the APACHE algorithms: cardiovascular, gastrointestinal, metabolic, neurologic, respiratory, trauma, renal/genitourinary (cardiovascular); mean annual ICU volume across all ICUs (0 < mean volume; 1 > mean volume); Died-in-ICU (0 = Alive, 1 = Died)
Clinically meaningful combinations of variables and their interactions were assessed for effect; to preserve clinical transparency, higher order interactions were explored, but not entertained in the final model. The potential for multiple colinearity was tested using the variance inflation factor (VIF) and condition number (CN); where VIF < 10 and CN less than “30 or more”  were desirable.
OLS using length of stay as a raw and a logged dependent variable. For the log OLS the expected value, E(y), is proportional to and the variance is assumed to be proportional to the mean squared. Prediction is on the log-scale (the geometric mean) and re-transformation is dependent upon the distribution of the error term: if normally distributed [39, 40]; under homoscedasticity where is the estimated smearing factor and is usually between 1 and 4 [41, 42].
LMM (logged length of stay) using maximum likelihood for model estimates. Potential modifying covariates were computed as fixed effects; ICU-year units as random intercepts (or “levels”) and random coefficients (“slopes”: APACHE III and APACHE III squared, age and ventilation status) were incorporated into the model fit.
A treatment effects model, log-dependent variable, via the Stata™ module “treatreg” as initially described by Cong and Drukker . A treatment-effects model is a two-equation system estimator (non-linear probit  and linear OLS), in which the effect of an endogenous binary variable (in this case, “Died-in-ICU”) on the continuous dependent variable is estimated by maximum likelihood [24, 45]. Formally, the model is expressed in two equations : the regression equation , where z j is the endogenous dummy variable indicating treatment assignment (=1, if ; or not, ). The outcome treatment-effect, estimated by the probit equation, is understood as being determined by an unobservable latent variable (), a linear function of covariates (w j ), and a random component (u j ): . The covariates entered into the outcome treatment equation were: age mechanical ventilation status, APACHE III score and its square, patient surgical status and the interaction with ventilation status, descriptors of ICU admission primary organ system dysfunction. As these covariates were commensurate with those of the OLS equation, a question of model identification arises . However, as shown by Maddala, such a model is identified even if the error terms are not independent and the variable list is identical in the two equations .
GLMs  with a log link (relating the conditional mean to the covariates) and either a gamma, inverse Gaussian or Poisson family (specifying the relationship between the variance and the mean) [14, 42] using maximum likelihood. The possibility of over-dispersion with the Poisson GLM was tested by calculating and regressing z as a constant-only model; significance of z indicating over-dispersion . In the presence of over-dispersion, the Poisson GLM model prediction and performance measures were replicated using a GLM negative binomial (NB2) model with log link . For this model the variance function is given as , where μ is the mean and α the heterogeneity or over-dispersion parameter, estimated by maximum likelihood (in GLM Poisson regression, α=0). The negative binomial may be also be conceptualised as a “waiting time” distribution (the number of days to discharge), a generalisation of the geometric distribution .
EEE model allowing simultaneous estimation of flexible parametric link (Box-Cox function ) and variance functions, as implemented by Basu in the Stata™ user-written module “pglm”, using the raw dependent variable, scaled to the mean [18, 49]. Expected values, conditional on covariates, are an inverse Box-Cox function of the linear predictor: . The default power variance function was used, which sets the family of variance functions to . The variance functions include the Poisson, gamma and inverse Gaussian variance functions as special cases ; as with linear OLS, where . Estimation of regression and link parameters is via an extension of quasi-likelihood, and the variance parameters by additional estimating equations; that is, the EEE is a semi-parametric model. The link function transforms the mean of the outcome (not the outcome), overcoming the retransformation problem .
Regression utilising the skew-t and skew-normal distribution as implemented by Marchenko and Genton in the Stata™ user-written modules “skewtreg” and “skewnreg” respectively . These distributions, fitted to the transformed (log) dependent variable, are (equivalently) the log-skew-t and log-skew-normal distributions by virtue of the correspondence which connects the normal and log-normal distributions [52, 53]. The skew-normal distribution, beside location and scale parameters, has an additional shape parameter (α; if α = 0 ≡ normal distribution) regulating distributional asymmetry; the skew-t distribution allows for asymmetry and heavier tails via shape (α) and degrees-of-freedom (v) parameters . Regression and other model parameters are estimated by maximum likelihood.
FMM, with two mixing components of a gamma and negative binomial (NB1 and NB2) densities, using the user-written Stata™ module “fmm”, version 2.0.0, by Deb . In a FMM, the random variable is assumed to be drawn from a super-population that is an additive mixture of distinct subpopulations or classes (C) in proportions where [20–22]. The C-point finite mixture model is given by , where the mixing probabilities πj are estimated with other parameters (Θ) [20–22]. The θ j is a parameter vector characterizing the component density function fj(), the component density belonging to different parametric families (implemented in “fmm” are gamma, lognormal, negative binomial, normal, Poisson and Student-t). Estimation is by maximum likelihood.
user-written Stata™ code provided by Norton  and Jones ; the “GLM family” test (in GLMs there is accommodation of skewness, in particular, via the variance : Vary|x = α.[E(y|x)]; therefore test if γ = 0 (Gaussian), = 1 (Poisson), = 2 (gamma), and = 3 (inverse Gaussian)); and the “COPAS” test: testing for over-fitting via split sample (50:50) cross-validation (1000x) using a version of the test initially formulated by Copas ; forecast predictions are regressed on raw-scale length of stay and a significant difference of βpredictions from unity suggests over-fitting ).
With respect to the particular problem that raw-scale variance was a power function of the raw-scale mean function, the (modified [40, 59]) Park test  was implemented. Squared residuals from a provisional model (GLM or OLS log-scale) were regressed (using a GLM with robust variance estimate and log-link) on the predictions () from the same model, both log transformed: ; the coefficient on γ(see (i), above) indicating the GLM variance function.
a modified Hosmer-Lemeshow goodness-of-fit (F) test [50, 64], regressing the particular model residual distribution against deciles of the linear predictor (using OLS with robust variance). Graphical plots of the calibration of regression decile-parameter estimates against deciles of the linear predictor were undertaken; positive parameter estimates above the null-line (≡ 0) indicate under-prediction and negative estimates indicate over-prediction .
As the interpretation of β coefficients is not invariant to model link functions , the average marginal effect [65, 66] of ICU mortality upon length of stay (“Died-in-ICU”) was estimated using the “margins” routine of Stata™ and the specific computations provided by Basu  for the EEE estimator. A marginal (“partial” or “incremental” ) effect measures the effect on the conditional mean of the regression dependent variable (“y”, ≡ ICU length-of-stay) of a change in a regressor(s) of interest (in this case the binary variable “Died-in-ICU”) . The impact of progressive increases in patient severity of illness (APACHE III score) was elucidated by the “marginsplot” module of Stata™ , using the contrast across the 1 (≡ Died)/0 (≡ Alive) binary variable “Died-in-ICU”. Because of the technical problem of applying the Duan estimator directly to a contrast of the predicted logs in the log-scale estimators (exponentiating the difference and then multiplying by a function of the residuals would not produce an estimate directly comparable to the contrast produced after regression performed on raw-scale data), the illustrative plots were generated using raw-scale OLS and a log-link Gaussian family GLM.
Stata™ (Version 12 MP, 2011; Stata Corporation, College Station, Texas) statistical software was used.
Demographics of ICU-hospital level, ventilation and surgical status: mean(SD)
ICU length of stay (days)
APACHE III score
ICU length of stay (days)
APACHE III score
ICU length of stay (days)
APACHE III score
ICU length of stay (days)
APACHE III score
Performance comparisons for various estimators
Determination data-set (n = 86740)
Validation data-set (n = 21678)
OLS: raw scale
+ ve skew
OLS: log scale
+ ve skew
LMM (log scale)
+ ve skew
Treatment effects regression (log scale)
+ ve skew
GLM: family(Poisson), link(log)
+ ve skew
GLM: family(negbin), link(log)
+ ve skew
GLM: family(gamma), link(log)
+ ve skew#
+ ve skew#
family(inverse Gaussian), link(log)
+ ve skew#
Log-skew- t regression
FMM: raw scale
+ ve skew#
Specification tests for various estimators
Hosmer-Lemeshow (F) test
LP vs residuals
OLS raw scale @
OLS log scale @
LMM log scale @
Treatment effects model: log scale
GLM: family(Poisson), link(log) #
GLM: family(negative binomial), link(log) #
GLM: family(gamma), link(log) #
GLM: family(inverse Gaussian), link(log) #
EEE: raw scale
Skew- t regression: log scale
Skew normal regression: log scale
1.134 (1.113, 1.155
FMM: raw scale
β regression estimate and average marginal effect for “Died in ICU”
Average marginal effect
OLS raw scale
2.947 (2.368, 3.526)
2.486 (2.060, 2.912)
OLS log scale
0.404 (0.302, 0.507)
0.353 (0.268, 0.437)
LMM log scale
0.367 (0.318, 0.416)
0.320 (0.285, 0.355)
Treatment effects regression: log scale
0.297 (0.194, 0.401)
0.215 (0.125, 0.302)
GLM: family(Poisson), link(log)
0.600 (0.444, 0.676)
1.456 (1.064, 1.848)
GLM: family(negative binomial), link(log)
0.596 (0.476, 0.715)
1.679 (1.239, 2.119)
GLM: family(gamma), link(log)
0.619 (0.497, 0.740)
1.802 (1.336, 2.269)
GLM: family(inverse Gaussian), link(log)
0.757 (0.608, 0.906)
2.648 (1.914, 3.382)
EEE: raw scale
0.678 (0.539, 0.817)
1.507 (1.042, 1.972
Skew-t regression: log scale
0.376 (0.268, 0.483)
0.335 (0.250, 0.420)
Skew normal regression: log scale
0.314 (0.214, 0.413)
0.271 (0.195, 0.348)
FMM: raw scale, component 1
0.234 (0.121, 0.347)
0.130 (0.066, 0.195)
FMM: raw scale, component 2
0.648 (0.513, 0.784)
2.395 (1.871, 2.920)
GLM: family(Gaussian), link(log)
0.541 (0.421, 0.662)
1.285 (0.902, 1.667)
Methodological recommendations for non-normal data analysis have usually been subsumed under the rubric of cost analysis [12, 15, 69]. The distribution of positive health expenditures (and length of stay) exhibit skew, kurtosis and heteroscedasticity, with a non-constant and increasing variance, albeit particular cost problems, such as the accommodation of a probability mass at “zero”, may not have immediate relevance to length of stay analysis [51, 69]. This review has suggested that a LMM had better overall performance compared with 11 other estimators. However, a number of issues pertain to comparator studies of possible estimators; in particular data structure and the full assessment, or otherwise, of model specification.
The hierarchical nature of some clinical and/or administrative data sets and the requirement for appropriate analysis has been commented upon above. As opposed to the volume-outcome literature , mixed models would appear to have been have been applied in a relatively small number of studies for length of stay analysis [16, 71–75]. The use of administrative data [5, 16] may also be associated with integer based calendar–day recording which lacks the accuracy of fractional day based “exact” times [10, 71].
Because of the non-normality of length of stay, formal trimming [76, 77] or truncation [71, 77, 78] of the data, or deleting [5, 73, 79] “outliers” has been undertaken prior to possible data transformation. For instance, studies implementing the APACHE III  or IV  algorithms for predicting ICU length of stay have truncated the latter at 30 days, at a 1% data level; the same fraction as seen with data deletions . The current study used a deletion fraction of 0.001%. The motivation for such strategies are various , but the theoretical basis of these data revisions has been questioned  and bias in the estimation of the mean has also been suggested [29, 69]. More importantly, large values of length of stay are “true” values  and data reduction cannot necessarily be defended as “…represent[ing] unusual cases…”, or because of “…analytical problems…” . That trimming will cause decreases in residual variance in LMM was pointed out by Lee et al.  and was easily demonstrated in the current data set; where the residual variance was 0.646 for the full data, 0.600 at the 99th percentile trim (< 26 days) and 0.508 for the 95th percentile trim (< 12 days) of ICU length of stay. Model assessment must take note of data revisions.
Similarly, the repeated assertions [71, 72, 76, 77, 80, 81] that the objective of such data manipulations are model (in particular, OLS) requirements for normality of the dependent variable are problematic. As observed by Buntin and Zaslavsky “Plotting of the data with various transformations is thus a useful preliminary step, but ultimately the distribution of the residuals from the transformed model is critical because the model assumptions concern the distribution of the residuals, not of the data” . That the residuals “…usually have a similar distribution to the original data…”  does not absolve the analyst from model based residual interrogation; the simple reporting of overall measures such as R2 or agreement indices of observed versus predicted  are not satisfactory proxies. Moreover, inference on β-coefficients is biased under heteroscedasticity.
Effect of ICU death on length of stay.
The β coefficient for “Died-in-ICU”, as a main effect, ranged from +29.7% (treatment effects regression, log scale) to +300% (OLS, raw scale), consistent with the β estimates for the factor “Death” of +0.536 (age < 65 years) and +0.439 (age ≥ 65 years) reported by Angus et al. , using a GLM (geometric family, identity link). However, the results of the marginal analysis (Table 4) implied over-prediction of the β-coefficient of “Died in ICU” across the estimators, albeit the average marginal effect was considered, as opposed to the marginal effect at the mean or at a particular value [65, 66]. Furthermore, there was a progressive decrease in the contrast, died in ICU versus alive in ICU, of predicted mean ICU days with an increase in the APACHE III score (Figure 7); this effect being consistent with findings from our previous study  and those of Rapaport et al.  and Woods et al. . Of more interest was the demonstration, using the treatment effects regression model, that the covariate “Died-in-ICU” (which also has the status of an “outcome”) was associated with a significant correlation (ρ) between the error terms of the separate estimations (probit and linear regression) and was appropriately considered as endogenous . That ρ was positive at 0.065 also indicated that OLS overestimated the “treatment” effect (Died-in-ICU); that is, a face-value interpretation of the β mortality coefficient was problematic. Endogeneity may also be suspected in mortality models where length of stay  or mortality probability  are entered as predictive covariates. Such regression of a variable upon its components has been termed a “dubious practice” . This being said, the extension of 2-stage methods to non-linear models may not be without its own issues [15, 84, 85].
Estimators of length of stay.
The statistical basis of the various estimators used in the current paper has been canvassed in “Statistical methods”, above. These estimators, assessing the mean function , have seen extensive application and assessment in the cost data literature, but not necessarily for length of stay. The recent introduction of log-skew-normal and log-skew-t estimators has not afforded opportunity for extensive assessment. Within the ICU literature, the predominant estimators have been ordinary least squares regression, with  and without [7, 71, 78] log transformation; generalized linear model (geometric family with identity link)  and random intercept linear model . Lee and co-workers have explored alternative approaches suited to the analysis of hospital length of stay [74, 86].
It is perhaps not surprising that in a large hierarchically structured data set that the LMM appeared to dominate. The suggestion of lack of discriminatory power across alternative estimators in cost analysis, albeit with small sample size , also appears to be vindicated in the current analysis. However, the further suggestion that large sample size together with “…simple methods…”  and an identity link [12, 15] may be analytically sufficient seems not to be the case with our data, although the LMM and log-OLS models were relatively easy to fit and interpret. The estimated link function (λ=0.152(0.019, 0.285) supports neither an identity (λ=1) nor log link (λ=0), rather a fractional- but not a square-root (λ=0.5) function. Although the Park test and the EEE estimation of the variance function approximated 2, the GLM variants demonstrated poor goodness-of-fit (Figure 5) with the exception of the EEE (not with-standing the upper most decile of the linear predictor, Figure 6) where the link and variance functions were estimated, not fixed. The narrow confidence limits in the Copas test (Table 3) for the EEE estimator presumably reflects this re-estimation of link and variance functions in the cross-validation process. The residual behaviour of (i) the log scale estimators (Table 3) was generally satisfactory and empirical estimation suggested a quadratic mean-variance relationship which was consistent with the log-OLS estimator and (ii) the GLM models, was variable and with heavy tails. The use of gamma distribution in the presence of heavy tails has recently been cautioned . Estimators on the raw scale (OLS and FMM) demonstrated consistent but sub-optimal residual behaviour; and despite modest BIC values, the log-skew estimators had an overall poor residual performance.
In the current study, the range of R2 was modest (0.18-0.22; Table 2); consistent with that previously reported, 0.13  to 0.21 [71, 83], and the observations of Diehr et al. that R2 for cost utilisation data are usually ≤ 0.20 . Increases in R2 have been reported, not surprisingly, across groupings; for instance, ICU units (R2 = 0.78; ). In the current study, R2 across ICU units (n = 118) was 0.92 (LMM, determination set).
Critique of methodology
The hierarchical structure of data under analysis is not uncommon in clinical or administrative data sets, but may not be seen in cost utilisation studies where survey and administrative data predominate. Thus cross paradigm performance assessments may lack validity, although the skew and kurtosis indices (4.4 and 29.4 respectively, in the current data) at least were comparable (4.1 and 25.6 , 4.4 and 50 ). The significance of the residual tests (Table 3) and the tendency to over-fitting must be interpreted in the context of the large data sets  and the reported adverse effect of outliers and skewed data on these tests [15, 42]. We chose not to benchmark goodness of fit against the behaviour of raw residuals from the log-scale estimators. This preference was dictated by the application of homoscedastic retransformation from log-scale estimation which may not have been optimal, leading to the generation of raw-scale residuals () of uncertain status. Such uncertainty may be contrasted with the fixed relationship (via the inverse logit function) between the scale of estimation (log-odds) and the scale of interest (probability) in logistic regression where the Hosmer-Lemeshow test was first described. No generally applicable method of heteroscedastic retransformation recommends itself with large data sets and multiple predictor variables [39, 42, 51, 59] and was not undertaken. Our analysis suggests that percentage changes in days are important, not absolute changes; this follows from the fact that the data satisfy important optimal statistical properties on the log scale, and that we can therefore estimate and interpret effects on this scale with some confidence.
Comparisons between estimators with respect to the significance (or otherwise), the effect magnitude and form (for continuous variables) of other covariates entered into the regressions (the covariates were fixed [11, 59]) were not a focus of interest and were not reported. The relatively poor performance of the GLM model variants may be subject to enhancement by the use of generalised linear mixed models  which could incorporate the random effects structure of the LMM. Similarly, our use of the Poisson and negative binomial GLM fails to account for the structural exclusion of zero counts (both distributions include zeros) and use of zero-truncated count models may have been appropriate .
Model choice under the conditions of the current data set would appear to favour the LMM and log-OLS. The treatment effects model afforded extra explanation with respect to the covariate “Died-in-ICU” and had a good fit to the data. Mechanistic approaches to estimation are represented by the FMM which captures the bi-modal nature of the data (Figures 1 and 3), although the fit was not optimal; and the negative binomial GLM which represented the length of stay as “waiting times”. Neither the log-skew-normal nor log-skew-t estimator would appear to recommend themselves as the preferred analytic approach for the current data set.
- Burns LR, Wholey DR: The effects of patient, hospital, and physician characteristics on length of stay and mortality. Med Care. 1991, 29: 251-271. 10.1097/00005650-199103000-00007.View ArticlePubMedGoogle Scholar
- Thomas JW, Guire KE, Horvat GG: Is patient length of stay related to quality of care?. Hosp Health Serv Adm. 1997, 42: 489-507.PubMedGoogle Scholar
- Needham DM, Anderson G, Pink GH, McKillop I, Tomlinson GA, Detsky AS: A province-wide study of the association between hospital resource allocation and length of stay. Health Serv Manage Res. 2003, 16: 155-166. 10.1258/095148403322167915.View ArticlePubMedGoogle Scholar
- Afessa B, Keegan MT, Hubmayr RD, Naessens JM, Gajic O, Long KH, Peters SG: Evaluating the performance of an institution using an intensive care unit benchmark. Mayo Clin Proc. 2005, 80: 174-180. 10.4065/80.2.174.View ArticlePubMedGoogle Scholar
- Angus DC, Linde-Zwirble WT, Sirio CA, Rotondi AJ, Chelluri L, Newbold RC, Lave JR, Pinsky MR: The effect of managed care on ICU length of stay: implications for medicare. JAMA. 1996, 276: 1075-1082. 10.1001/jama.1996.03540130073033.View ArticlePubMedGoogle Scholar
- Becker RB, Zimmerman JE, Knaus WA, Wagner DP, Seneff MG, Draper EA, Higgins TL, Estafanous FG, Loop FD: The use of APACHE III to evaluate ICU length of stay, resource use, and mortality after coronary artery by-pass surgery. J Cardiovasc Surg (Torino). 1995, 36: 1-11.Google Scholar
- Knaus WA, Wagner DP, Zimmerman JE, Draper EA: Variations in mortality and length of stay in intensive care units. Ann Intern Med. 1993, 118: 753-761.View ArticlePubMedGoogle Scholar
- Rapoport J, Teres D, Zhao Y, Lemeshow S: Length of stay data as a guide to hospital economic performance for ICU patients. Med Care. 2003, 41: 386-397.PubMedGoogle Scholar
- Woods AW, MacKirdy FN, Livingston BM, Norrie J, Howie JC: Evaluation of predicted and actual length of stay in 22 Scottish intensive care units using the APACHE III system. Acute Physiology and Chronic Health Evaluation. Anaesthesia. 2000, 55: 1058-1065. 10.1046/j.1365-2044.2000.01552.x.View ArticlePubMedGoogle Scholar
- Marik PE, Hedman L: What's in a day? Determining intensive care unit length of stay. Crit Care Med. 2000, 28: 2090-2093. 10.1097/00003246-200006000-00071.View ArticlePubMedGoogle Scholar
- Austin PC, Rothwell DM, Tu JV: A comparison of statistical modeling strategies for analyzing length of stay after CABG surgery. Health Services & Outcomes Research Methodology. 2002, 3: 107-133. 10.1023/A:1024260023851.View ArticleGoogle Scholar
- Diehr P, Yanez D, Ash A, Hornbrook M, Lin DY: Methods for analyzing health care utilization and costs. Annu Rev Public Health. 1999, 20: 125-144. 10.1146/annurev.publhealth.20.1.125.View ArticlePubMedGoogle Scholar
- Manning WG, Basu A, Mullahy J: Generalized modeling approaches to risk adjustment of skewed outcomes data. J Health Econ. 2005, 24: 465-488. 10.1016/j.jhealeco.2004.09.011.View ArticlePubMedGoogle Scholar
- Moran JL, Solomon PJ, Peisach AR, Martin J: New models for old questions: Generalized Linear Models for cost prediction. J Eval Clin Pract. 2007, 13: 381-389. 10.1111/j.1365-2753.2006.00711.x.View ArticlePubMedGoogle Scholar
- Basu AP, Manning WGP: Issues for the next generation of health care cost analyses. Med Care. 2009, 47: S109-S114. 10.1097/MLR.0b013e31819c94a1.View ArticlePubMedGoogle Scholar
- Carey K: Hospital length of stay and cost: a multilevel modeling analysis. Health Services & Outcomes Research Methodology. 2002, 3: 41-56. 10.1023/A:1021530924455.View ArticleGoogle Scholar
- Stow PJ, Hart GK, Higlett T, George C, Herkes R, McWilliam D, Bellomo R: Development and implementation of a high-quality clinical database: the Australian and New Zealand intensive care society adult patient database. J Crit Care. 2006, 21: 133-141. 10.1016/j.jcrc.2005.11.010.View ArticlePubMedGoogle Scholar
- Basu A, Rathouz PJ: Estimating marginal and incremental effects on health outcomes using flexible link and variance function models. Biostat. 2005, 6: 93-109. 10.1093/biostatistics/kxh020.View ArticleGoogle Scholar
- Marchenko YV, Genton MG: A suite of commands for fitting the skew-normal and skew-t models. The Stata Journal. 2011, 10: 507-539.Google Scholar
- Conway KS, Deb P: Is prenatal care really ineffective? Or, is the `devil' in the distribution?. J Health Econ. 2005, 24: 489-513. 10.1016/j.jhealeco.2004.09.012.View ArticlePubMedGoogle Scholar
- Deb P, Holmes AM: Estimates of use and costs of behavioural health care: a comparison of standard and finite mixture models. Health Econ. 2000, 9: 475-489. 10.1002/1099-1050(200009)9:6<475::AID-HEC544>3.0.CO;2-H.View ArticlePubMedGoogle Scholar
- Singh CH, Ladusingh L: Inpatient length of stay: a finite mixture modeling analysis. Eur J Health Econ. 2010, 11: 119-126. 10.1007/s10198-009-0153-6.View ArticlePubMedGoogle Scholar
- Terza J: Estimating endogenous treatment effects in retrospective data analysis. Value in health: the journal of the International Society for Pharmacoeconomics and Outcomes Research. 1999, 2: 429-434. 10.1046/j.1524-4733.1999.26003.x.View ArticleGoogle Scholar
- Guo S, Fraser MW: Sample selection and related models. Propensity Score Analysis: Statistical Methods and Applications. Edited by: Guo S, Fraser MW. 2010, Thousand Oaks: SAGE Publications, Inc, 85-124.Google Scholar
- Wooldridge JM: Limited dependent variable models and sample selection corrections. Introductory econometrics: A modern approach. Edited by: Wooldridge JM, Mason OH. 2006, Thomson South-Western, 582-631. ThirdGoogle Scholar
- Hilbe JM: Handling endogeneity and latent class models. Negative Binomial Regression. 2011, Cambridge: Cambridge University Press, 407-446. 2ndView ArticleGoogle Scholar
- ANZICS: Adult Data Base. Data Dictionary Version 1.5. 2005, @ http://www.anzics.com.au/admc/files/data_dictionarypdf; Accessed, June 2006Google Scholar
- Wagner D, Knaus W, Bergner M: Statistical-Methods. Crit Care Med. 1989, 17: S194-S198.View ArticlePubMedGoogle Scholar
- Lee AH, Xiao J, Vemuri SR, Zhao Y: A discordancy test approach to identify outliers of length of hospital stay. Stat Med. 1998, 17: 2199-2206. 10.1002/(SICI)1097-0258(19981015)17:19<2199::AID-SIM917>3.0.CO;2-2.View ArticlePubMedGoogle Scholar
- Buis ML: HANGROOT: Stata module creating a hanging rootogram comparing an empirical distribution to the best fitting theoretical distribution. @ http://econpapers.repec.org/scripts/search/searchasp?ft=hangroot 2011, Accessed June 2011
- Wainer H: The suspended rootogram and other visual displays: an empirical validation. The American Statistician. 1974, 28: 143-145.Google Scholar
- Winter N, Nichols A: VIOPLOT: Stata module to produce violin plots. @ http://econpapers.repec.org/scripts/search/searchasp?ft=hangroot 2010, Accessed June 2010
- Hintze JL, Nelson RD: Violin plots: a box plot-density trace synergism. The American Statistician. 1998, 52: 181-184.Google Scholar
- Cox NJ: Speaking Stata: density probability plots. Stata Journal. 2005, 5: 259-273.Google Scholar
- Cox NJ: Speaking Stata: the limits of sample skewness and kurtosis. The Stata Journal. 2010, 10: 492-495.Google Scholar
- Knaus WA, Wagner DP, Draper EA, Zimmerman JE, Bergner M, Bastos PG, Sirio CA, Murphy DJ, Lotring T, Damiano A: The APACHE III prognostic system. Risk prediction of hospital mortality for critically ill hospitalized adults. Chest. 1991, 100: 1619-1636. 10.1378/chest.100.6.1619.View ArticlePubMedGoogle Scholar
- Belsley DA: Conditioning diagnostics, collinearity and weak data in regression. 1991, New York: John Wiley & SonsGoogle Scholar
- Rogers W: sg17: Regression standard errors in clustered samples. Stata Technical Bulletin Reprints. 1993, 3: 88-94.Google Scholar
- Manning WG: The logged dependent variable, heteroscedasticity, and the retransformation problem. J Health Econ. 1998, 17: 283-295. 10.1016/S0167-6296(98)00025-3.View ArticlePubMedGoogle Scholar
- Manning WG, Mullahy J: Estimating log models: to transform or not to transform?. J Health Econ. 2001, 20: 461-494. 10.1016/S0167-6296(01)00086-8.View ArticlePubMedGoogle Scholar
- Duan N: Smearing estimate: a nonparametric retransformation method. J Am Stat Assoc. 1983, 78: 605-610. 10.1080/01621459.1983.10478017.View ArticleGoogle Scholar
- Jones AM: Models for Health Care. Health Econometrics and Data Group Working Paper 10/01. 2010, @ york.ac.uk/res/herc/hedgwp, Accessed May 2011Google Scholar
- Cong R, Drukker DM: Treatment-effects model. Stata Technical Bulletin Reprints. 2000, 10: 159-169.Google Scholar
- Cheung YB: Adjustment for selection bias in cohort studies: an application of a probit model with selectivity to life course epidemiology. J Clin Epidemiol. 2001, 54: 1238-1243. 10.1016/S0895-4356(01)00403-6.View ArticlePubMedGoogle Scholar
- Maddala GS: Multivariate qualitative variables. Limited-Dependent and Qualitative Variables in Econometrics. Edited by: Maddala GS. 1983, Cambridgeshire: Cambridge University Press, 93-148.View ArticleGoogle Scholar
- Hardin J, Hilbe J: The Poisson family. Generalized linear models and extensions. 2007, College Station: TX: Stata Press, 183-198. 2ndGoogle Scholar
- Hilbe JM: Negative Binomial Regression. 2011, Cambridge, UK: Cambridge University Press, 2ndView ArticleGoogle Scholar
- Box GEP, Cox DR: An Analysis of Transformations. Journal of the Royal Statistical Society Series B-Statistical Methodology. 1964, 26: 211-252.Google Scholar
- Basu A: Extended beneralized linear models: Simulataneous estimation of flexible link and variance functions. Stata Journal. 2005, 5: 501-516.Google Scholar
- Basu A, Arondekar BV, Rathouz PJ: Scale of interest versus scale of estimation: comparing alternative estimators for the incremental costs of a comorbidity. Health Econ. 2006, 15: 1091-1107. 10.1002/hec.1099.View ArticlePubMedGoogle Scholar
- Hill SC, Miller GE: Health expenditure estimation and functional form: applications of the generalized gamma and extended estimating equations models. Health Econ. 2010, 19: 608-627.PubMedGoogle Scholar
- Azzalini A, dalCapello T, Kotz S: Log-skew-normal and log-skew-t distributions as models for family income data. Journal of Income Distribution. 2003, 11: 12-20.Google Scholar
- Walls WD: Modelling heavy tails and skewness in film returns. Applied Financial Economics. 2005, 15: 1181-1188. 10.1080/0960310050391040.View ArticleGoogle Scholar
- Deb P, FMM : Stata module to estimate finite mixture models. 2008, @ http://fmwww.bc.edu/repec/bocode/f/fmm 2008, Accessed November 2010Google Scholar
- Kuha J: AIC and BIC. Comparisons of assumptions and performance. Sociological Methods & Research. 2005, 33: 188-229.View ArticleGoogle Scholar
- Scott A, Wild C: Transformations and R2. The American Statistician. 1991, 45: 127-129.Google Scholar
- Thomas JW, Ashcraft ML: Measuring severity of illness: six severity systems and their ability to explain cost variations. Inquiry. 1991, 28: 39-55.PubMedGoogle Scholar
- Lin LI: A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989, 45: 255-268. 10.2307/2532051.View ArticlePubMedGoogle Scholar
- Buntin MB, Zaslavsky AM: Too much ado about two-part models and transformation? Comparing methods of modeling Medicare expenditures. J Health Econ. 2004, 23: 525-542. 10.1016/j.jhealeco.2003.10.005.View ArticlePubMedGoogle Scholar
- Norton EC: Stata code for modeling health care costs. @ http://www.uncedu/~enorton/ 2005, Accessed January 2011
- Jones AM: Modelling non-normal outcomes using linear models. @ http://melbourneinstitute.com/health/WorkshopApr202010.html 2010, Accessed May 2011
- Copas JB: Regression, prediction and shrinkage. Journal of the Royal Statistical Society Series B. 1983, 45: 311-354.Google Scholar
- Park RE: Estimation with heteroscedastic error terms. Econometrica. 1966, 34: 888-10.2307/1910108.View ArticleGoogle Scholar
- Hosmer DW, Taber S, Lemeshow S: The importance of assessing the fit of logistic regression models: a case study. Am J Public Health. 1991, 81: 1630-1635. 10.2105/AJPH.81.12.1630.View ArticlePubMedPubMed CentralGoogle Scholar
- Bartus T: Estimation of marginal effects using margeff. Stata Journal. 2011, 5: 309-329.Google Scholar
- Cameron CC, Trivedi PK: Nonlinear regression methods. Microeconomics using Stata. 2010, College Station: Stata Press, 319-362.Google Scholar
- Stata Corp: Reference Manual G-M: Release 12. 2011, College Station: Stata Corp LPGoogle Scholar
- Hardin J, Hilbe J: Analysis of fit. Generalized linear models and extensions. Edited by: Hardin J, Hilbe J. 2007, College Station, TX: Stata Press, 47-62. SecondGoogle Scholar
- Mihaylova B, Briggs A, O'Hagan A, Thompson SG: Review of statistical methods for analysing healthcare resources and costs. Health Econ. 2010, 20: 897-916.View ArticlePubMedPubMed CentralGoogle Scholar
- Panageas KS, Schrag D, Riedel E, Bach PB, Begg CB: The effect of clustering of outcomes on the association of procedure volume and surgical outcomes. Ann Intern Med. 2003, 139: 658-665.View ArticlePubMedGoogle Scholar
- Kramer AA, Zimmerman JE: The relationship between hospital and intensive care unit length of stay. Crit Care Med. 2011, 39: 1015-1022. 10.1097/CCM.0b013e31820eabab.View ArticlePubMedGoogle Scholar
- Lee AH, Gracey M, Wang K, Yau KKW: A robustified modeling approach to analyze pediatric length of stay. Ann Epidemiol. 2005, 15: 673-677. 10.1016/j.annepidem.2004.10.001.View ArticlePubMedGoogle Scholar
- Leung KM, Elashoff RM, Rees KS, Hasan MM, Legorreta AP: Hospital- and patient-related characteristics determining maternity length of stay: a hierarchical linear model approach. Am J Public Health. 1998, 88: 377-381. 10.2105/AJPH.88.3.377.View ArticlePubMedPubMed CentralGoogle Scholar
- Yau KKW, Lee AH, Ng ASK: Finite mixture regression model with random effects: application to neonatal hospital length of stay. Computational Statistics & Data Analysis. 2003, 41: 359-366. 10.1016/S0167-9473(02)00180-9.View ArticleGoogle Scholar
- Yau KKW, Wang K, Lee AH: Zero-inflated negative binomial mixed regression modeling of over-dispersed count data with extra zeros. Biom J. 2003, 45: 437-452. 10.1002/bimj.200390024.View ArticleGoogle Scholar
- Lee AH, Fung WK, Fu B: Analyzing hospital length of stay - Mean or median regression?. Med Care. 2003, 41: 681-686.PubMedGoogle Scholar
- Liu VM, Kipnis PP, Gould MKM, Escobar GJM: Length of stay predictions: improvements through the use of automated laboratory and comorbidity variables. Med Care. 2010, 48: 739-744. 10.1097/MLR.0b013e3181e359f3.View ArticlePubMedGoogle Scholar
- Rosenberg AL, Zimmerman JE, Alzola C, Draper EA, Knaus WA: Intensive care unit length of stay: recent changes and future challenges. Crit Care Med. 2000, 28: 3465-3473. 10.1097/00003246-200010000-00016.View ArticlePubMedGoogle Scholar
- Huntley DA, Cho DW, Christman J, Csernansky JG: Predicting length of stay in an acute psychiatric hospital. Psychiatr Serv. 1998, 49: 1049-1053.View ArticlePubMedGoogle Scholar
- Marazzi A, Paccaud F, Ruffieux C, Beguin C: Fitting the distributions of length of stay by parametric models. Med Care. 1998, 36: 915-927. 10.1097/00005650-199806000-00014.View ArticlePubMedGoogle Scholar
- Nathanson BH: Making progress with the egress. Would P.T. Barnum make a good hospital administrator?. Crit Care Med. 2011, 39: 1208-1209. 10.1097/CCM.0b013e318211fa81.View ArticlePubMedGoogle Scholar
- Moran JL, Bristow P, Solomon PJ, George C, Hart GK: for the Australian and New Zealand Intensive Care Society Database Management Committee (ADMC): Mortality and length-of-stay outcomes, 1993–2003, in the binational Australian and New Zealand intensive care adult patient database. Crit Care Med. 2008, 36: 46-61. 10.1097/01.CCM.0000295313.08084.58.View ArticlePubMedGoogle Scholar
- Render ML, Kim HM, Deddens J, Sivaganesin S, Welsh DE, Bickel K, Freyberg R, Timmons S, Johnston J, Connors AF, Wagner D, Hofer TP: Variation in outcomes in veterans affairs intensive care units with a computerized severity measure. Crit Care Med. 2005, 33: 930-939. 10.1097/01.CCM.0000162497.86229.E9.View ArticlePubMedGoogle Scholar
- Cai B, Small DS, Have TRT: Two-stage instrumental variable methods for estimating the causal odds ratio: analysis of bias. Stat Med. 2011, 30: 1809-1824. 10.1002/sim.4241.View ArticlePubMedGoogle Scholar
- Terza JV, Basu A, Rathouz PJ: Two-stage residual inclusion estimation: addressing endogeneity in health econometric modeling. J Health Econ. 2008, 27: 531-543. 10.1016/j.jhealeco.2007.09.009.View ArticlePubMedGoogle Scholar
- Lee AH, Wang K, Yau KK, Somerford PJ: Truncated negative binomial mixed regression modelling of ischaemic stroke hospitalizations. Stat Med. 2003, 22: 1129-1139. 10.1002/sim.1419.View ArticlePubMedGoogle Scholar
- Rowan KM, Kerr JH, Major E, McPherson K, Short A, Vessey MP: Intensive Care Society's APACHE II study in Britain and Ireland–I: variations in case mix of adult admissions to general intensive care units and impact on outcome. BMJ. 1993, 307: 972-977. 10.1136/bmj.307.6910.972.View ArticlePubMedPubMed CentralGoogle Scholar
- Rabe-Hesketh S, Skrondal A: Generalized Linear Mixed Models. International Encyclopedia of Education. Edited by: Peterson P, Baker E, McGaw B, Barry M. 2010, Oxford: Elsevier, 171-177. ThirdthView ArticleGoogle Scholar
- The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2288/12/68/prepub