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

Background 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. 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 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. Results 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. Conclusions 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.


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 [4][5][6][7][8][9], with the majority of studies presenting cross-sectional 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 [11][12][13][14]; 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 data-bases, 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 length-of-stay); OLS (with both raw and log-scaled 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 skew-normal and skew-t multivariate distributions [19]; and finite mixture (FMM) models which consider the dependent variable as a mixture of distributions [20][21][22]; and (ii) determine the effect of (ICU) mortality outcome upon length of stay, allowing for "endogenous variable bias" [23] 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 [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 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 [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 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 [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 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) [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 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" [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 exp x 0 β À Á 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μ x i ð Þ ¼φ: exp x 0 iβ 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; 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. 3. A treatment effects model, log-dependent variable, via the Stata™ module "treatreg" as initially described by Cong and Drukker [43]. A treatment-effects model is a two-equation system estimator (non-linear probit [44] 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 [43]: the regression equation The outcome treatment-effect, estimated by the probit equation, is understood as being determined by an unobservable latent variable (z Ã j ), a linear function of covariates (w j ), and a random component (u j ): z Ã j ¼ w j γ þ 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 over-dispersion with the Poisson GLM was tested by calculating z ¼ and regressing z as a constant-only model; significance of z indicating over-dispersion [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 μ þ αμ 2 , 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 [5]. 5. EEE model allowing simultaneous estimation of flexible parametric link (Box-Cox function [48]) 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 h μ i ; θ 1 ; θ 2 ð Þ ð Þto The variance functions include the Poisson, gamma and inverse Gaussian variance functions as special cases [50]; as with linear OLS, where λ ¼ 1 and θ 2 ¼ 0 [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 skew-t and skew-normal distribution as implemented by Marchenko and Genton in the Stata™ user-written modules "skewtreg" and "skewnreg" respectively [19]. 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 [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 user-written Stata™ module "fmm", version 2.0.0, by Deb [54]. 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 [20][21][22]. The C-point finite mixture model is given by f ðy i jΘÞ ¼ P c j¼1 π i f j ðy i jθ j Þ, where the mixing probabilities π j are estimated with other parameters (Θ) [20][21][22]. 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 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) [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) 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 [41]. Distributional form of the raw-scale 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) user-written 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]: ; 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 [62]; forecast predictions are regressed on raw-scale length of stay and a significant difference of β predictions from unity suggests over-fitting [42]). (ii) 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 [63] 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: [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 [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 ("Died-in-ICU") 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]  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") [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 "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.

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.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   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 effectsand 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 R 2 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 (  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.

Discussion
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 [70], mixed models would appear to have been have been applied in a relatively small number of studies for length of stay analysis [16,[71][72][73][74][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 [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 (See figure on previous page.) Figure 5 Upper panel: Violin plots of ICU length of stay (raw-scale, length of stay < 16 days) and predicted length of stay (raw-scale, < 16 days; determination and validation data sets) for estimators: GLM Poisson-log, GLM negative binomial-log, GLM gamma-log, GLM inverse Gaussian-log. Lower panel: Hosmer-Lemeshow (F) test of residuals across deciles of linear-predictor, for estimators: GLM Poisson-log, GLM negative binomial-log, GLM gamma-log, GLM inverse Gaussian-log. X-axis, linear predictor; Y-axis, parameter estimates. determin., determination. Validat., validation. ". . .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 "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. [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 "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 [26]. 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 [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 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 u x ð Þ , have seen extensive application and assessment in the cost data literature, but not necessarily for length of stay. The recent introduction of log-skewnormal 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 [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 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 [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 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 squareroot (λ=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 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 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 [69]. Estimators on the raw scale (OLS and FMM) demonstrated consistent but suboptimal residual behaviour; and despite modest BIC values, the log-skew estimators had an overall poor residual performance.

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 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 rawscale residuals ( y i Àŷ i ) 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 [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 zero-truncated 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 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.