- Research article
- Open Access
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
BMC Medical Research Methodology volume 12, Article number: 68 (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 .
Predictor variables considered were
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.
ICU length of stay was modelled (with a fixed covariate list, and clustering on ICU site with robust variance ) via the estimators;
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.
Model adequacy was gauged by (i) progressive reduction in AIC (Akaike Information Criterion, for nested models) and BIC (Bayesian Information Criterion, for non-nested models) , 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 . Residuals were generated from the fitted model according to the scale of estimation . 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) calendar-years. 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) back-transformation to the original metric (days) utilised Duan’s smearing estimate . Distributional form of the raw-scale predictions (determination and validation) were displayed using violin plots. Final model performance was assessed by R2 (on the “day” scale ) computed as the square of the correlation between predicted and observed length of stay ; Lin’s concordance correlation coefficient (between raw and predicted raw length of stay ); mean absolute error (MAE); and root mean squared error (RMSE) [12, 14]. Supplementary diagnostic analyses [15, 40, 50, 51, 59] were also performed, using:
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.
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.05-199.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, inter-quartile range 2.8 (0.93-3.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).
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 goodness-of-fit), 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 log-scale estimators (LMM, treatment effects- and log-skew-regressions) 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 R2 and concordance correlation coefficient was 0.22 and 0.334 respectively.
The outcome discrimination (Died-in-ICU) of the probit sub-equation 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, over-fitting, 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 (Vary|x = α.[E(y|x)]γ ) 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 log-scale OLS; 1.659(1.609, 1.708) for GLM Poisson-log; 2.106(2.005, 2.207) for GLM gamma-log; and 2.352(2.088, 2.616) for GLM inverse Gaussian-log. 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 log-skew-normal and log-skew-t regressions reported α values of 1.073(0.933, 1.213) and 1.490(1.323, 1.659) respectively, suggesting positive-skew; the degrees-of-freedom for the log-skew-t regression was 10.53(7.97, 13.91) implying heavier-than-normal tails for the conditional distribution of (log) ICU length-of-stay. 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 non-convergence.
Figures 4, 5, 6 show (i) upper panels: violin-plot distributional form of predictions and (ii) lower panels: Hosmer-Lemeshow-test parameter estimates, for determination and validation sets for each of the estimators, grouped by increasing scalar values of the Hosmer-Lemeshow parameter estimates. The graphical displays show a progressive lack-of-fit (Hosmer-Lemeshow test estimates) across Figures 4, 5, 6 in concert with corresponding over- and under-prediction as revealed by the violin plots. The LMM, log-scale 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.
The β regression estimates for the binary variable “Died-in-ICU” are seen in Table 4, final column, and range from a low of 29% (treatment effects regression) to 295% (raw-scale 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 “Died-in-ICU” across increments of APACHE III score for raw-scale OLS and log-link 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 log-ICU length of stay.
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.
Thomas JW, Guire KE, Horvat GG: Is patient length of stay related to quality of care?. Hosp Health Serv Adm. 1997, 42: 489-507.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Basu AP, Manning WGP: Issues for the next generation of health care cost analyses. Med Care. 2009, 47: S109-S114. 10.1097/MLR.0b013e31819c94a1.
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.
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.
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.
Marchenko YV, Genton MG: A suite of commands for fitting the skew-normal and skew-t models. The Stata Journal. 2011, 10: 507-539.
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.
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.
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.
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.
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.
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. Third
Hilbe JM: Handling endogeneity and latent class models. Negative Binomial Regression. 2011, Cambridge: Cambridge University Press, 407-446. 2nd
ANZICS: Adult Data Base. Data Dictionary Version 1.5. 2005, @ http://www.anzics.com.au/admc/files/data_dictionarypdf; Accessed, June 2006
Wagner D, Knaus W, Bergner M: Statistical-Methods. Crit Care Med. 1989, 17: S194-S198.
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.
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.
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.
Cox NJ: Speaking Stata: density probability plots. Stata Journal. 2005, 5: 259-273.
Cox NJ: Speaking Stata: the limits of sample skewness and kurtosis. The Stata Journal. 2010, 10: 492-495.
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.
Belsley DA: Conditioning diagnostics, collinearity and weak data in regression. 1991, New York: John Wiley & Sons
Rogers W: sg17: Regression standard errors in clustered samples. Stata Technical Bulletin Reprints. 1993, 3: 88-94.
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.
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.
Duan N: Smearing estimate: a nonparametric retransformation method. J Am Stat Assoc. 1983, 78: 605-610. 10.1080/01621459.1983.10478017.
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
Cong R, Drukker DM: Treatment-effects model. Stata Technical Bulletin Reprints. 2000, 10: 159-169.
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.
Maddala GS: Multivariate qualitative variables. Limited-Dependent and Qualitative Variables in Econometrics. Edited by: Maddala GS. 1983, Cambridgeshire: Cambridge University Press, 93-148.
Hardin J, Hilbe J: The Poisson family. Generalized linear models and extensions. 2007, College Station: TX: Stata Press, 183-198. 2nd
Hilbe JM: Negative Binomial Regression. 2011, Cambridge, UK: Cambridge University Press, 2nd
Box GEP, Cox DR: An Analysis of Transformations. Journal of the Royal Statistical Society Series B-Statistical Methodology. 1964, 26: 211-252.
Basu A: Extended beneralized linear models: Simulataneous estimation of flexible link and variance functions. Stata Journal. 2005, 5: 501-516.
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.
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.
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.
Walls WD: Modelling heavy tails and skewness in film returns. Applied Financial Economics. 2005, 15: 1181-1188. 10.1080/0960310050391040.
Deb P, FMM : Stata module to estimate finite mixture models. 2008, @ http://fmwww.bc.edu/repec/bocode/f/fmm 2008, Accessed November 2010
Kuha J: AIC and BIC. Comparisons of assumptions and performance. Sociological Methods & Research. 2005, 33: 188-229.
Scott A, Wild C: Transformations and R2. The American Statistician. 1991, 45: 127-129.
Thomas JW, Ashcraft ML: Measuring severity of illness: six severity systems and their ability to explain cost variations. Inquiry. 1991, 28: 39-55.
Lin LI: A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989, 45: 255-268. 10.2307/2532051.
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.
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.
Park RE: Estimation with heteroscedastic error terms. Econometrica. 1966, 34: 888-10.2307/1910108.
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.
Bartus T: Estimation of marginal effects using margeff. Stata Journal. 2011, 5: 309-329.
Cameron CC, Trivedi PK: Nonlinear regression methods. Microeconomics using Stata. 2010, College Station: Stata Press, 319-362.
Stata Corp: Reference Manual G-M: Release 12. 2011, College Station: Stata Corp LP
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. Second
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.
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.
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.
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.
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.
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.
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.
Lee AH, Fung WK, Fu B: Analyzing hospital length of stay - Mean or median regression?. Med Care. 2003, 41: 681-686.
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.
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.
Huntley DA, Cho DW, Christman J, Csernansky JG: Predicting length of stay in an acute psychiatric hospital. Psychiatr Serv. 1998, 49: 1049-1053.
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.
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.
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.
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.
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.
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.
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.
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.
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. Thirdth
The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2288/12/68/prepub
The authors declare that they have no competing interests.
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.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Moran, J.L., Solomon, P.J. & the ANZICS Centre for Outcome and Resource Evaluation (CORE) of the Australian and New Zealand Intensive Care Society (ANZICS). 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. BMC Med Res Methodol 12, 68 (2012). https://doi.org/10.1186/1471-2288-12-68
- Intensive Care Unit
- Ordinary Little Square
- Intensive Care Unit Admission
- Variance Function
- Intensive Care Unit Mortality