Abstract
Background
For the analysis of lengthofstay (LOS) data, which is characteristically rightskewed, a number of statistical estimators have been proposed as alternatives to the traditional ordinary least squares (OLS) regression with log dependent variable.
Methods
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 riskadjusted regression analysis of firstly: log LOS using OLS, linear mixed model [LMM], treatment effects, skewnormal and skewt models; and secondly: unmodified (raw) LOS via OLS, generalised linear models [GLMs] with loglink and 4 different distributions [Poisson, gamma, negative binomial and inverseGaussian], extended estimating equations [EEE] and a finite mixture model including a gamma distribution. A fixed covariate list and ICUsite clustering with robust variance were utilised for model fitting with splitsample determination (80%) and validation (20%) data sets, and model simulation was undertaken to establish overfitting (Copas test). Indices of model specification using Bayesian information criterion [BIC: lower values preferred] and residual analysis as well as predictive performance (R^{2}, concordance correlation coefficient (CCC), mean absolute error [MAE]) were established for each estimator.
Results
The dataset 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 lengthofstay was 3.4(5.1) (median 1.8, range (0.1760)) days and demonstrated marked kurtosis and right skew (29.4 and 4.4 respectively). BIC showed considerable spread, from a maximum of 509801 (OLSraw scale) to a minimum of 210286 (LMM). R^{2} ranged from 0.22 (LMM) to 0.17 and the CCC from 0.334 (LMM) to 0.149, with MAE 2.22.4. Superior residual behaviour was established for the logscale estimators. There was a general tendency for overprediction (negative residuals) and for overfitting, the exception being the GLM negative binomial estimator. The meanvariance function was best approximated by a quadratic function, consistent with logscale estimation; the link function was estimated (EEE) as 0.152(0.019, 0.285), consistent with a fractionalroot function.
Conclusions
For ICU length of stay, logscale estimation, in particular the LMM, appeared to be the most consistently performing estimator(s). Neither the GLM variants nor the skewregression estimators dominated.
Background
Length of stay during an intensive care unit (ICU) or hospital admission is a function of diverse patient and organisational input variables [1]. It is widely used as an indicator of performance [2] and is a determinant of costs, although resource allocation is also known to affect length of stay [3]. Not surprisingly, ICU length of stay has been the subject of frequent analysis [49], with the majority of studies presenting crosssectional analyses over a relatively short periods of months [10] to 1–2 years [9].
ICU patient length of stay (and costs) demonstrate skewed distribution and various statistical modelling strategies have been employed in analysis of such data [1114]; albeit linear regression (ordinary least squares regression, OLS) of the logged dependent variable has demonstrated a remarkable persistence [15]. Individual patient data, as accessed from ICU databases, have an intrinsic hierarchical structure (patients within ICUs) and due analytic consideration of this structure is also appropriate [16]. Using such data from the Australian and New Zealand Intensive Care (ANZICS) adult patient database (APD) [17], 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 lengthofstay); OLS (with both raw and logscaled dependent variable) and generalised linear models (GLMs [14]), 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 [18]; estimators utilising skewnormal and skewt multivariate distributions [19]; and finite mixture (FMM) models which consider the dependent variable as a mixture of distributions [2022]; and (ii) determine the effect of (ICU) mortality outcome upon length of stay, allowing for “endogenous variable bias” [23] by means of a treatmenteffects model [24,25] where specific allowance is made for the correlation of independent predictors and error terms. Under such conditions, model coefficients are biased [26].
Methods
The ANZICS adult patient database (APD) is a binational (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 [27]. The current database 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 feedback, resubmission and incorporation into a final reporting dataset.
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 [28]. 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 [29].
Statistical methods
Continuous variables were reported using mean (SD), except where otherwise indicated; categorical variables were analysed using the Chisquared 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 [32]. 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 [33], incorporates a marker for the median of the data, a box indicating the interquartile range, and spikes extending to the upper and loweradjacent 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) [34]. 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 [35].
Predictor variables considered were
a) Continuous: age (in years), severity of illness score (Acute Physiology and Chronic Health Evaluation (APACHE) III score [36]); both centred for computation. The predictive effect of these variables was entered initially as both linear and quadratic.
b) 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 postelective surgery, postemergency surgery and nonsurgical (nonsurgical); 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); DiedinICU (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” [37] were desirable.
ICU length of stay was modelled (with a fixed covariate list, and clustering on ICU site with robust variance [38]) via the estimators;
1. 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 logscale (the geometric mean) and retransformation 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].
2. LMM (logged length of stay) using maximum likelihood for model estimates. Potential modifying covariates were computed as fixed effects; ICUyear 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.
3. A treatment effects model, logdependent variable, via the Stata™ module “treatreg” as initially described by Cong and Drukker [43]. A treatmenteffects model is a twoequation system estimator (nonlinear probit [44] and linear OLS), in which the effect of an endogenous binary variable (in this case, “DiedinICU”) on the continuous dependent variable is estimated by maximum likelihood [24,45]. Formally, the model is expressed in two equations [43]: the regression equation , where z_{j} is the endogenous dummy variable indicating treatment assignment (=1, if ; or not, ). The outcome treatmenteffect, 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 [25]. 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 [45].
4. GLMs [14] 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 overdispersion with the Poisson GLM was tested by calculating and regressing z as a constantonly model; significance of z indicating overdispersion [46]. In the presence of overdispersion, the Poisson GLM model prediction and performance measures were replicated using a GLM negative binomial (NB2) model with log link [47]. For this model the variance function is given as , where μ is the mean and α the heterogeneity or overdispersion 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 [5].
5. EEE model allowing simultaneous estimation of flexible parametric link (BoxCox function [48]) and variance functions, as implemented by Basu in the Stata™ userwritten module “pglm”, using the raw dependent variable, scaled to the mean [18,49]. Expected values, conditional on covariates, are an inverse BoxCox 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 [50]; as with linear OLS, where [51]. Estimation of regression and link parameters is via an extension of quasilikelihood, and the variance parameters by additional estimating equations; that is, the EEE is a semiparametric model. The link function transforms the mean of the outcome (not the outcome), overcoming the retransformation problem [50].
6. Regression utilising the skewt and skewnormal distribution as implemented by Marchenko and Genton in the Stata™ userwritten modules “skewtreg” and “skewnreg” respectively [19]. These distributions, fitted to the transformed (log) dependent variable, are (equivalently) the logskewt and logskewnormal distributions by virtue of the correspondence which connects the normal and lognormal distributions [52,53]. The skewnormal distribution, beside location and scale parameters, has an additional shape parameter (α; if α = 0 ≡ normal distribution) regulating distributional asymmetry; the skewt distribution allows for asymmetry and heavier tails via shape (α) and degreesoffreedom (v) parameters [19]. Regression and other model parameters are estimated by maximum likelihood.
7. FMM, with two mixing components of a gamma and negative binomial (NB1 and NB2) densities, using the userwritten Stata™ module “fmm”, version 2.0.0, by Deb [54]. In a FMM, the random variable is assumed to be drawn from a superpopulation that is an additive mixture of distinct subpopulations or classes (C) in proportions where [2022]. The Cpoint finite mixture model is given by , where the mixing probabilities π_{j} are estimated with other parameters (Θ) [2022]. The θ_{j} is a parameter vector characterizing the component density function f_{j}(), the component density belonging to different parametric families (implemented in “fmm” are gamma, lognormal, negative binomial, normal, Poisson and Studentt). Estimation is by maximum likelihood.
Model adequacy was gauged by (i) progressive reduction in AIC (Akaike Information Criterion, for nested models) and BIC (Bayesian Information Criterion, for nonnested models) [55], both of which are penalised (with respect to number of observations and model parameters) likelihood methods for model determination (ii) likelihood ratio tests where appropriate and (iii) normality and lack of heteroscedasticity of residuals based upon graphical analysis [14]. Residuals were generated from the fitted model according to the scale of estimation [50]. The final model was developed on a determination set (80% of data) and validated on a validation set (20% of data); the random samples being stratified by (the 2) calendaryears. Predicted values were generated with continuous covariates centred and categorical covariates held at the reference category, as above, and in the presence of a logged dependent variable (length of stay) backtransformation to the original metric (days) utilised Duan’s smearing estimate [41]. Distributional form of the rawscale predictions (determination and validation) were displayed using violin plots. Final model performance was assessed by R^{2} (on the “day” scale [56]) computed as the square of the correlation between predicted and observed length of stay [57]; Lin’s concordance correlation coefficient (between raw and predicted raw length of stay [58]); mean absolute error (MAE); and root mean squared error (RMSE) [12,14]. Supplementary diagnostic analyses [15,40,50,51,59] were also performed, using:
(i) userwritten Stata™ code provided by Norton [60] and Jones [61]; the “GLM family” test (in GLMs there is accommodation of skewness, in particular, via the variance [42]: Varyx = α.[E(yx)]; therefore test if γ = 0 (Gaussian), = 1 (Poisson), = 2 (gamma), and = 3 (inverse Gaussian)); and the “COPAS” test: testing for overfitting via split sample (50:50) crossvalidation (1000x) using a version of the test initially formulated by Copas [62]; forecast predictions are regressed on rawscale length of stay and a significant difference of β_{predictions} from unity suggests overfitting [42]).
(ii) With respect to the particular problem that rawscale variance was a power function of the rawscale mean function, the (modified [40,59]) Park test [63] was implemented. Squared residuals from a provisional model (GLM or OLS logscale) were regressed (using a GLM with robust variance estimate and loglink) on the predictions () from the same model, both log transformed: ; the coefficient on γ(see (i), above) indicating the GLM variance function.
(iii) a modified HosmerLemeshow goodnessoffit (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 decileparameter estimates against deciles of the linear predictor were undertaken; positive parameter estimates above the nullline (≡ 0) indicate underprediction and negative estimates indicate overprediction [50].
As the interpretation of β coefficients is not invariant to model link functions [18], the average marginal effect [65,66] of ICU mortality upon length of stay (“DiedinICU”) was estimated using the “margins” routine of Stata™ and the specific computations provided by Basu [49] for the EEE estimator. A marginal (“partial” or “incremental” [18]) effect measures the effect on the conditional mean of the regression dependent variable (“y”, ≡ ICU lengthofstay) of a change in a regressor(s) of interest (in this case the binary variable “DiedinICU”) [66]. The impact of progressive increases in patient severity of illness (APACHE III score) was elucidated by the “marginsplot” module of Stata™ [67], using the contrast across the 1 (≡ Died)/0 (≡ Alive) binary variable “DiedinICU”. Because of the technical problem of applying the Duan estimator directly to a contrast of the predicted logs in the logscale 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 rawscale data), the illustrative plots were generated using rawscale OLS and a loglink Gaussian family GLM.
Stata™ (Version 12 MP, 2011; Stata Corporation, College Station, Texas) statistical software was used.
Results
The database for the period 2008–2009 contained records of 114798 patients with exclusion (4.20%) of patients having: ICU length of stay < 4 hours, incomplete GCS or APACHE III records. The final data set had missing ICU length of stay in 0.03% and 141 patients (0.001%) with ICU length of stay > 60 days (mean 81.98(24.37), range 60.05199.67) not considered in analysis. Patients (n = 111663) admitted to 131 ICUs in the same number of hospitals, were of mean(SD) age 60.6(18.8) years, APACHE III score 52.5(28.9); 43.0% were female, 40.7% were mechanically ventilated (on the day of admission) and ICU mortality was 7.8%. ICU length of stay was 3.4(5.1) (median 1.8, interquartile range 2.8 (0.933.7)) days. Summary statistics of length of stay, APACHE III score, age and gender by ventilation status and patient surgical type are shown in Table 1. Overall ICU length of stay and the log transform are displayed in histogram form with normal curve overlay in Figure 1 (left and right panels respectively); raw scale length of stay demonstrated marked kurtosis and right skew (29.4 and 4.4 respectively) whilst the log transform demonstrated a more symmetrical distribution (kurtosis 3.1 and skewness 0.5). Comparison of the empirical distribution of ICU length of stay and its’ log and square root transformation showed poor approximation to standard theoretical distributions (normal, lognormal and gamma; Figure 2).
Table 1. Demographics of ICUhospital level, ventilation and surgical status: mean(SD)
Figure 1. Histograms of ICU length of stay with normal curve overlay: left panel, raw scale; right panel, log transform.
Figure 2. Hanging rootograms displaying approximations of the empirical distribution of ICU length of stay, and transformations (log and squareroot), to various theoretical distributions. Shaded areas a tip of spiles are 95% CI.
Comparative performance of the estimators for ICU length of stay is seen in Tables 2 and 3; the regression models had 52 fixed parameters, including the constant term. Over the range of measures assessing predictive performance and model specification (including goodnessoffit), the LMM demonstrated a consistently better performance compared with the other 11 estimators; albeit all estimators revealed some problematic aspects of residual analyses. BIC for all the logscale estimators (LMM, treatment effects and logskewregressions) were considerably less than for estimation on the raw scale (Table 2). The final model variable set displayed a condition number of 19.4 and a mean VIF of 3.76; the best R^{2} and concordance correlation coefficient was 0.22 and 0.334 respectively.
Table 2. Performance comparisons for various estimators
Table 3. Specification tests for various estimators
The outcome discrimination (DiedinICU) of the probit subequation of the treatment effects regression was excellent for both determination and validation sets (ROC curve area 0.92). The null hypothesis, that the correlation (ρ=0.065) between the error terms of the separate estimations is zero, was rejected at P = 0.0001. The Poisson GLM exhibited overdispersion (P value for “z” = 0.0001), the negative binomial GLM demonstrating an expected decrease in BIC; other performance measures were comparable between these two models. However, overfitting, as assessed by the Copas test, was demonstrated in all other models except for the negative binomial GLM (Table 3). The “GLM family test” of the variance function of GLMs (Varyx = α.[E(yx)]γ [60]) rejected tests of γ for 0, 1, 2 and 3 (P ≤ 0.0001). The modified Park test generated estimates for γ of: 0.348(0.333, 0.364) for logscale OLS; 1.659(1.609, 1.708) for GLM Poissonlog; 2.106(2.005, 2.207) for GLM gammalog; and 2.352(2.088, 2.616) for GLM inverse Gaussianlog. The EEE link parameter (λ) estimate was 0.152(0.019, 0.285) and the variance function (θ_{2}) was 2.115(2.013, 2.218). The logskewnormal and logskewt regressions reported α values of 1.073(0.933, 1.213) and 1.490(1.323, 1.659) respectively, suggesting positiveskew; the degreesoffreedom for the logskewt regression was 10.53(7.97, 13.91) implying heavierthannormal tails for the conditional distribution of (log) ICU lengthofstay. The 2 component FMM gamma model generated predicted mean ICU days of 0.56(0.32), range 0.06, 3.64 (component 1) and 3.69(2.45), range 0.49, 19.1 (component 2) as illustrated in Figure 3. The FMM negative binomial models demonstrated nonconvergence.
Figure 3. Histograms of predicted mean ICU days for component 1 (left panel) and component 2 (right panel) of the FMM.
Figures 4, 5, 6 show (i) upper panels: violinplot distributional form of predictions and (ii) lower panels: HosmerLemeshowtest parameter estimates, for determination and validation sets for each of the estimators, grouped by increasing scalar values of the HosmerLemeshow parameter estimates. The graphical displays show a progressive lackoffit (HosmerLemeshow test estimates) across Figures 4, 5, 6 in concert with corresponding over and underprediction as revealed by the violin plots. The LMM, logscale OLS and treatment effects regression were the best calibrated; the EEE showed good calibration across the linear predictor deciles except for the upper deciles where substantial discordance was evident.
Figure 4. Upper panel: Violin plots of ICU length of stay (rawscale, length of stay < 16 days) and predicted length of stay (rawscale, < 16 days; determination and validation data sets) for estimators: OLS (rawscale), logscale OLS, LMM, and treatment effects regression. Lower panel: HosmerLemeshow (F) test of residuals across deciles of linearpredictor, for estimators: OLS (rawscale), logscale OLS, LMM, and treatment effects regression. Xaxis, linear predictor; Yaxis, parameter estimates.
Figure 5. Upper panel: Violin plots of ICU length of stay (rawscale, length of stay < 16 days) and predicted length of stay (rawscale, < 16 days; determination and validation data sets) for estimators: GLM Poissonlog, GLM negative binomiallog, GLM gammalog, GLM inverse Gaussianlog. Lower panel: HosmerLemeshow (F) test of residuals across deciles of linearpredictor, for estimators: GLM Poissonlog, GLM negative binomiallog, GLM gammalog, GLM inverse Gaussianlog. Xaxis, linear predictor; Yaxis, parameter estimates. determin., determination. Validat., validation.
Figure 6. Upper panel: Violin plots of ICU length of stay (rawscale, length of stay < 16 days) and predicted length of stay (raw scale, < 16 days; determination and validation data sets) for estimators: EEE, logskewt regression, logskewnormal regression and FMM. Lower panel: HosmerLemeshow (F) test of residuals across deciles of linearpredictor, for estimators: EEE, logskewt regression, logskewnormal regression and FMM. Xaxis, linear predictor; Yaxis, parameter estimates.
The β regression estimates for the binary variable “DiedinICU” are seen in Table 4, final column, and range from a low of 29% (treatment effects regression) to 295% (rawscale OLS).The average marginal effects are also tabulated, and were generally decreased in magnitude compared with the β estimates. Figure 7 illustrates the contrasts of predictive margins for “DiedinICU” across increments of APACHE III score for rawscale OLS and loglink Gaussian family GLM (see also Table 4). Similar decreases (not shown) of the contrasts of predictive margins across the APACHE III score range were evident for the (log) linear predictor for estimators using logICU length of stay.
Discussion
Methodological recommendations for nonnormal 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 nonconstant 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 volumeoutcome literature [70], mixed models would appear to have been have been applied in a relatively small number of studies for length of stay analysis [16,7175]. 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 nonnormality 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 [78] or IV [71] 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 [73]. The current study used a deletion fraction of 0.001%. The motivation for such strategies are various [80], but the theoretical basis of these data revisions has been questioned [29] 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 [69] and data reduction cannot necessarily be defended as “…represent[ing] unusual cases…”, or because of “…analytical problems…” [73]. That trimming will cause decreases in residual variance in LMM was pointed out by Lee et al. [72] and was easily demonstrated in the current data set; where the residual variance was 0.646 for the full data, 0.600 at the 99^{th} percentile trim (< 26 days) and 0.508 for the 95^{th} 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” [59]. That the residuals “…usually have a similar distribution to the original data…” [12] does not absolve the analyst from model based residual interrogation; the simple reporting of overall measures such as R^{2}[77] or agreement indices of observed versus predicted [71] are not satisfactory proxies. Moreover, inference on βcoefficients is biased under heteroscedasticity.
Effect of ICU death on length of stay.
The β coefficient for “DiedinICU”, 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. [5], using a GLM (geometric family, identity link). However, the results of the marginal analysis (Table 4) implied overprediction 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 [82] and those of Rapaport et al. [8] and Woods et al. [9]. Of more interest was the demonstration, using the treatment effects regression model, that the covariate “DiedinICU” (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 [26]. That ρ was positive at 0.065 also indicated that OLS overestimated the “treatment” effect (DiedinICU); that is, a facevalue interpretation of the β mortality coefficient was problematic. Endogeneity may also be suspected in mortality models where length of stay [7] or mortality probability [83] are entered as predictive covariates. Such regression of a variable upon its components has been termed a “dubious practice” [15]. This being said, the extension of 2stage methods to nonlinear 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 logskewnormal and logskewt estimators has not afforded opportunity for extensive assessment. Within the ICU literature, the predominant estimators have been ordinary least squares regression, with [8] and without [7,71,78] log transformation; generalized linear model (geometric family with identity link) [5] and random intercept linear model [83]. Lee and coworkers 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 [15], also appears to be vindicated in the current analysis. However, the further suggestion that large sample size together with “…simple methods…” [69] and an identity link [12,15] may be analytically sufficient seems not to be the case with our data, although the LMM and logOLS 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 squareroot (λ=0.5) function. Although the Park test and the EEE estimation of the variance function approximated 2, the GLM variants demonstrated poor goodnessoffit (Figure 5) with the exception of the EEE (not withstanding 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 reestimation of link and variance functions in the crossvalidation process. The residual behaviour of (i) the log scale estimators (Table 3) was generally satisfactory and empirical estimation suggested a quadratic meanvariance relationship which was consistent with the logOLS 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 [69]. Estimators on the raw scale (OLS and FMM) demonstrated consistent but suboptimal residual behaviour; and despite modest BIC values, the logskew estimators had an overall poor residual performance.
In the current study, the range of R^{2} was modest (0.180.22; Table 2); consistent with that previously reported, 0.13 [78] to 0.21 [71,83], and the observations of Diehr et al. that R^{2} for cost utilisation data are usually ≤ 0.20 [12]. Increases in R^{2} have been reported, not surprisingly, across groupings; for instance, ICU units (R^{2} = 0.78; [7]). In the current study, R^{2} 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 [42], 4.4 and 50 [50]). The significance of the residual tests (Table 3) and the tendency to overfitting must be interpreted in the context of the large data sets [87] 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 logscale estimators. This preference was dictated by the application of homoscedastic retransformation from logscale estimation which may not have been optimal, leading to the generation of rawscale residuals () of uncertain status. Such uncertainty may be contrasted with the fixed relationship (via the inverse logit function) between the scale of estimation (logodds) and the scale of interest (probability) in logistic regression where the HosmerLemeshow 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 [88] 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 zerotruncated count models may have been appropriate [47].
Conclusion
Model choice under the conditions of the current data set would appear to favour the LMM and logOLS. The treatment effects model afforded extra explanation with respect to the covariate “DiedinICU” and had a good fit to the data. Mechanistic approaches to estimation are represented by the FMM which captures the bimodal 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 logskewnormal nor logskewt estimator would appear to recommend themselves as the preferred analytic approach for the current data set.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
The study was conceived, designed, (data)analysed, written and critically revised jointly by both authors (JLM, PJS). Both authors read and approved the final manuscript.
References

Burns LR, Wholey DR: The effects of patient, hospital, and physician characteristics on length of stay and mortality.
Med Care 1991, 29:251271. PubMed Abstract  Publisher Full Text

Thomas JW, Guire KE, Horvat GG: Is patient length of stay related to quality of care?
Hosp Health Serv Adm 1997, 42:489507. PubMed Abstract

Needham DM, Anderson G, Pink GH, McKillop I, Tomlinson GA, Detsky AS: A provincewide study of the association between hospital resource allocation and length of stay.
Health Serv Manage Res 2003, 16:155166. PubMed Abstract  Publisher Full Text

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:174180. PubMed Abstract  Publisher Full Text

Angus DC, LindeZwirble WT, Sirio CA, Rotondi AJ, Chelluri L, Newbold RC III, Lave JR, Pinsky MR: The effect of managed care on ICU length of stay: implications for medicare.
JAMA 1996, 276:10751082. PubMed Abstract  Publisher Full Text

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 bypass surgery.

Knaus WA, Wagner DP, Zimmerman JE, Draper EA: Variations in mortality and length of stay in intensive care units.
Ann Intern Med 1993, 118:753761. PubMed Abstract  Publisher Full Text

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:386397. PubMed Abstract  Publisher Full Text

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:10581065. PubMed Abstract  Publisher Full Text

Marik PE, Hedman L: What's in a day? Determining intensive care unit length of stay.
Crit Care Med 2000, 28:20902093. PubMed Abstract  Publisher Full Text

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:107133. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Diehr P, Yanez D, Ash A, Hornbrook M, Lin DY: Methods for analyzing health care utilization and costs.
Annu Rev Public Health 1999, 20:125144. PubMed Abstract  Publisher Full Text

Manning WG, Basu A, Mullahy J: Generalized modeling approaches to risk adjustment of skewed outcomes data.
J Health Econ 2005, 24:465488. PubMed Abstract  Publisher Full Text

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:381389. PubMed Abstract  Publisher Full Text

Basu AP, Manning WGP: Issues for the next generation of health care cost analyses.
Med Care 2009, 47:S109S114. PubMed Abstract  Publisher Full Text

Carey K: Hospital length of stay and cost: a multilevel modeling analysis.
Health Services & Outcomes Research Methodology 2002, 3:4156. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stow PJ, Hart GK, Higlett T, George C, Herkes R, McWilliam D, Bellomo R: Development and implementation of a highquality clinical database: the Australian and New Zealand intensive care society adult patient database.
J Crit Care 2006, 21:133141. PubMed Abstract  Publisher Full Text

Basu A, Rathouz PJ: Estimating marginal and incremental effects on health outcomes using flexible link and variance function models.
Biostat 2005, 6:93109. Publisher Full Text

Marchenko YV, Genton MG: A suite of commands for fitting the skewnormal and skewt models.

Conway KS, Deb P: Is prenatal care really ineffective? Or, is the `devil' in the distribution?
J Health Econ 2005, 24:489513. PubMed Abstract  Publisher Full Text

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:475489. PubMed Abstract  Publisher Full Text

Singh CH, Ladusingh L: Inpatient length of stay: a finite mixture modeling analysis.
Eur J Health Econ 2010, 11:119126. PubMed Abstract  Publisher Full Text

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:429434. Publisher Full Text

Guo S, Fraser MW: Sample selection and related models. In Propensity Score Analysis: Statistical Methods and Applications. Edited by Guo S, Fraser MW. Thousand Oaks: SAGE Publications, Inc; 2010:85124.

Wooldridge JM: Limited dependent variable models and sample selection corrections. In Introductory econometrics: A modern approach. Third edition. Edited by Wooldridge JM, Mason OH. Thomson SouthWestern; 2006:582631.

Hilbe JM: Handling endogeneity and latent class models. In Negative Binomial Regression. 2nd edition. Cambridge: Cambridge University Press; 2011:407446.

ANZICS: Adult Data Base. Data Dictionary Version 1.5.
2005.
@ http://www anzics com au/admc/files/data_dictionary pdf webcite; Accessed, June 2006

Wagner D, Knaus W, Bergner M: StatisticalMethods.
Crit Care Med 1989, 17:S194S198. PubMed Abstract

Lee AH, Xiao J, Vemuri SR, Zhao Y: A discordancy test approach to identify outliers of length of hospital stay.
Stat Med 1998, 17:21992206. PubMed Abstract  Publisher Full Text

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/search asp?ft=hangroot webcite 2011, Accessed June 2011

Wainer H: The suspended rootogram and other visual displays: an empirical validation.

Winter N, Nichols A: VIOPLOT: Stata module to produce violin plots. @ http://econpapers repec org/scripts/search/search asp?ft=hangroot webcite 2010, Accessed June 2010

Hintze JL, Nelson RD: Violin plots: a box plotdensity trace synergism.

Cox NJ: Speaking Stata: the limits of sample skewness and kurtosis.

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:16191636. PubMed Abstract  Publisher Full Text

Belsley DA: Conditioning diagnostics, collinearity and weak data in regression. New York: John Wiley & Sons; 1991.

Rogers W: sg17: Regression standard errors in clustered samples.

Manning WG: The logged dependent variable, heteroscedasticity, and the retransformation problem.
J Health Econ 1998, 17:283295. PubMed Abstract  Publisher Full Text

Manning WG, Mullahy J: Estimating log models: to transform or not to transform?
J Health Econ 2001, 20:461494. PubMed Abstract  Publisher Full Text

Duan N: Smearing estimate: a nonparametric retransformation method.
J Am Stat Assoc 1983, 78:605610. Publisher Full Text

Jones AM: Models for Health Care.
Health Econometrics and Data Group Working Paper 10/01 2010.
@ york.ac.uk/res/herc/hedgwp, Accessed May 2011

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:12381243. PubMed Abstract  Publisher Full Text

Maddala GS: Multivariate qualitative variables. In LimitedDependent and Qualitative Variables in Econometrics. Edited by Maddala GS. Cambridgeshire: Cambridge University Press; 1983:93148.

Hardin J, Hilbe J: The Poisson family. In Generalized linear models and extensions. 2nd edition. College Station: TX: Stata Press; 2007:183198.

Hilbe JM: Negative Binomial Regression. 2nd edition. Cambridge, UK: Cambridge University Press; 2011.

Box GEP, Cox DR: An Analysis of Transformations.
Journal of the Royal Statistical Society Series BStatistical Methodology 1964, 26:211252.

Basu A: Extended beneralized linear models: Simulataneous estimation of flexible link and variance functions.

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:10911107. PubMed Abstract  Publisher Full Text

Hill SC, Miller GE: Health expenditure estimation and functional form: applications of the generalized gamma and extended estimating equations models.
Health Econ 2010, 19:608627. PubMed Abstract  Publisher Full Text

Azzalini A, dalCapello T, Kotz S: Logskewnormal and logskewt distributions as models for family income data.

Walls WD: Modelling heavy tails and skewness in film returns.
Applied Financial Economics 2005, 15:11811188. Publisher Full Text

Stata module to estimate finite mixture models. 2008.
@ http://fmwww bc edu/repec/bocode/f/fmm webcite 2008, Accessed November 2010

Kuha J: AIC and BIC. Comparisons of assumptions and performance.
Sociological Methods & Research 2005, 33:188229. PubMed Abstract  Publisher Full Text

Thomas JW, Ashcraft ML: Measuring severity of illness: six severity systems and their ability to explain cost variations.
Inquiry 1991, 28:3955. PubMed Abstract

Lin LI: A concordance correlation coefficient to evaluate reproducibility.
Biometrics 1989, 45:255268. PubMed Abstract  Publisher Full Text

Buntin MB, Zaslavsky AM: Too much ado about twopart models and transformation? Comparing methods of modeling Medicare expenditures.
J Health Econ 2004, 23:525542. PubMed Abstract  Publisher Full Text

Norton EC: Stata code for modeling health care costs. @ http://www unc edu/~enorton/ webcite 2005, Accessed January 2011

Jones AM: Modelling nonnormal outcomes using linear models. @ http://melbourneinstitute com/health/WorkshopApr202010 html webcite 2010, Accessed May 2011

Copas JB: Regression, prediction and shrinkage.
Journal of the Royal Statistical Society Series B 1983, 45:311354.

Park RE: Estimation with heteroscedastic error terms.
Econometrica 1966, 34:888. Publisher Full Text

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:16301635. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cameron CC, Trivedi PK: Nonlinear regression methods. In Microeconomics using Stata. College Station: Stata Press; 2010:319362.

Stata Corp: Reference Manual GM: Release 12. College Station: Stata Corp LP; 2011.

Hardin J, Hilbe J: Analysis of fit. In Generalized linear models and extensions. Second edition. Edited by Hardin J, Hilbe J. College Station, TX: Stata Press; 2007:4762.

Mihaylova B, Briggs A, O'Hagan A, Thompson SG: Review of statistical methods for analysing healthcare resources and costs.
Health Econ 2010, 20:897916. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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:658665. PubMed Abstract  Publisher Full Text

Kramer AA, Zimmerman JE: The relationship between hospital and intensive care unit length of stay.
Crit Care Med 2011, 39:10151022. PubMed Abstract  Publisher Full Text

Lee AH, Gracey M, Wang K, Yau KKW: A robustified modeling approach to analyze pediatric length of stay.
Ann Epidemiol 2005, 15:673677. PubMed Abstract  Publisher Full Text

Leung KM, Elashoff RM, Rees KS, Hasan MM, Legorreta AP: Hospital and patientrelated characteristics determining maternity length of stay: a hierarchical linear model approach.
Am J Public Health 1998, 88:377381. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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:359366. PubMed Abstract  Publisher Full Text

Yau KKW, Wang K, Lee AH: Zeroinflated negative binomial mixed regression modeling of overdispersed count data with extra zeros.
Biom J 2003, 45:437452. Publisher Full Text

Lee AH, Fung WK, Fu B: Analyzing hospital length of stay  Mean or median regression?
Med Care 2003, 41:681686. PubMed Abstract  Publisher Full Text

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:739744. PubMed Abstract  Publisher Full Text

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:34653473. PubMed Abstract  Publisher Full Text

Huntley DA, Cho DW, Christman J, Csernansky JG: Predicting length of stay in an acute psychiatric hospital.
Psychiatr Serv 1998, 49:10491053. PubMed Abstract  Publisher Full Text

Marazzi A, Paccaud F, Ruffieux C, Beguin C: Fitting the distributions of length of stay by parametric models.
Med Care 1998, 36:915927. PubMed Abstract  Publisher Full Text

Nathanson BH: Making progress with the egress. Would P.T. Barnum make a good hospital administrator?
Crit Care Med 2011, 39:12081209. PubMed Abstract  Publisher Full Text

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 lengthofstay outcomes, 1993–2003, in the binational Australian and New Zealand intensive care adult patient database.
Crit Care Med 2008, 36:4661. PubMed Abstract  Publisher Full Text

Render ML, Kim HM, Deddens J, Sivaganesin S, Welsh DE, Bickel K, Freyberg R, Timmons S, Johnston J, Connors AF Jr, Wagner D, Hofer TP: Variation in outcomes in veterans affairs intensive care units with a computerized severity measure.
Crit Care Med 2005, 33:930939. PubMed Abstract  Publisher Full Text

Cai B, Small DS, Have TRT: Twostage instrumental variable methods for estimating the causal odds ratio: analysis of bias.
Stat Med 2011, 30:18091824. PubMed Abstract  Publisher Full Text

Terza JV, Basu A, Rathouz PJ: Twostage residual inclusion estimation: addressing endogeneity in health econometric modeling.
J Health Econ 2008, 27:531543. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lee AH, Wang K, Yau KK, Somerford PJ: Truncated negative binomial mixed regression modelling of ischaemic stroke hospitalizations.
Stat Med 2003, 22:11291139. PubMed Abstract  Publisher Full Text

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:972977. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

RabeHesketh S, Skrondal A: Generalized Linear Mixed Models. In International Encyclopedia of Education. Thirdth edition. Edited by Peterson P, Baker E, McGaw B, Barry M. Oxford: Elsevier; 2010:171177.
Prepublication history
The prepublication history for this paper can be accessed here: