Skip to main content

Statistical models for the analysis of skewed healthcare cost data: a simulation study

Abstract

Skewed data is the main issue in statistical models in healthcare costs. Data transformation is a conventional method to decrease skewness, but there are some disadvantages. Some recent studies have employed generalized linear models (GLMs) and Cox proportional hazard regression as alternative estimators.

The aim of this study was to investigate how well these alternative estimators perform in terms of bias and precision when the data are skewed. The primary outcome was an estimation of population means of healthcare costs and the secondary outcome was the impact of a covariate on healthcare cost. Alternative estimators, such as ordinary least squares (OLS) for Ln(y) or Log(y), Gamma, Weibull and Cox proportional hazard regression models, were compared using Monte Carlo simulation under different situations, which were generated from skewed distributions.

We found that there was not one best model across all generated conditions. However, GLMs, especially the Gamma regression model, behaved well in the estimation of population means of healthcare costs. The results showed that the Cox proportional hazard model exhibited a poor estimation of population means of healthcare costs and the β1 even under proportional hazard data. Approximately results are consistent by increasing the sample size. However, increasing the sample size could improve the performance of the OLS-based model.

Background

Statistical models are often used in many healthcare economics and policy studies. The main issues in such studies are the estimation of mean population healthcare costs and finding the best relationship between costs and covariates through regression modeling [1]. However, these cannot be implemented by simple statistical models as the healthcare costs data have specific characterizations [2]. Healthcare costs data demonstrate the substantial positive skewness and are sometimes characterized by the use of large resources with zero cost [3]. These specifications of data impose a number of difficulties in using standard statistical analysis, such as implementing linear regression causes unreliable results [2].

Two-part models based on mixture models are performed when excess zeroes are present in data [3]. Further, logarithmic (or other) transformations are commonly used to decrease the skewness and drive them close to normal distribution, in order to implement linear regression models. The logarithmic transformation with ordinary least squares (OLS) regression is a very common approach in applied economics. However, it also presents several drawbacks. One of these drawbacks is that the predictions are not robust enough to detect the heteroscedasticity in the transformed scale [1,4]. The general consensus is that estimating the mean cost using a logarithmic regression model leads to biased estimation [2,4-6].

An alternative approach is using nonlinear regression models, of which exponential conditional mean (ECM) models in generalized linear models (GLMs) are examples [7]. Generally, GLMs extend the linear modeling framework to allow response variables that are not normally distributed. In healthcare studies, generalized linear modeling through log-link function avoids the weakness and problems of OLS regression. In addition, the Cox proportional hazards model has been a controversial issue for healthcare data modeling. It has been used as a special flexible model for skewed healthcare data in many studies [8,9].

In recent years, extensive research efforts have been done to propose suitable regression methods for the analysis of skewed healthcare data [1,3,10,11]. Many studies also set out a clear framework for comparing these methods from a variety of aspects [5,6,12,13]. Moreover, a few have provided prominent reviews of the statistical methods for analyzing healthcare data [2,7].

However, there is no comparative study that investigates the different methods using different sample sizes. This paper was conducted to compare the proposed statistical models in the available literature using different sample sizes. We specifically focused on comparing proposed statistical models for positive skewed healthcare costs, but not zero mass problems. It was developed based on a Monte Carlo simulation to find appropriate methods to get the unbiased and precise estimates of the mean costs. This aspect is particularly important in the literature [5,13]. Furthermore, in this paper, the coefficient estimations of covariates are also evaluated in our simulations using different sample sizes.

Methods

Let y i denote healthcare expenditures for person i, and x i denote the set of covariates, including the intercept. We estimated the following models.

Ordinary least square based on log transformation

It is common to use linear regression models for the log of costs in healthcare expenditures. Logarithmic transformation is most commonly used to decrease skewness and to make the distribution more symmetric and closer to normality. The log regression model is as follows:

$$ \ln \left({y}_i\right)={x}_i\beta +{\varepsilon}_i $$

Where it was assumed that E() = 0 and E(ε) = 0, since predicting costs on the original scale is primary objective so:

$$ {y}_i= \exp \left({x}_i\beta +{\varepsilon}_i\right) $$
$$ E\left(\left.{y}_i\right|{x}_i\right)=E\left( \exp \left({\varepsilon}_i\right)\Big|{x}_i\right) \exp \left({x}_i\beta \right) $$

If the error term is \( N\left(0,{\sigma}_{\upvarepsilon}^2\right) \) distribution, it is a log-normal case, and then:

$$ E\left(\left.{y}_i\right|{x}_i\right)= \exp \left({x}_i\beta +0.5{\sigma}_{\upvarepsilon}^2\right) $$

However, if the error term is not normally distributed, but is homoscedastic, then the smearing estimator is applied.

Generalized linear models

GLMs are a broad class of statistical models for relating non-normal dependent variables to linear combinations of predictor variables. An invertible link function (g (.)) converts the expectation of the response variable, E (Y i ), to the linear predictor:

$$ g\left(E\left({y}_i\right)\right)=g\left({\mu}_i\right)={x}_i\beta $$

The ECM model is a special type of GLM with log-link function, and can be viewed as a nonlinear regression model:

$$ E\left(\left.{y}_i\right|{x}_i\right)= \exp \left({x}_i\beta \right) $$

Weibull and Gamma regression models are assumed as two special types of ECM model; β values were estimated here using quasi-maximum likelihood estimation. The exponential distribution was considered to be a special case of the Weibull and Gamma regression models when the shape parameter was equal to 1.

Cox proportional hazard model

The Cox proportional hazard model is based on hazard and survival functions, instead of ECM or direct estimation of E (y|x). It is a semi-parametric model because it does not specify the baseline hazard function:

$$ h\left(\left.{y}_i\right|{x}_i\right)={h}_0(y) \exp \left({x}_i\beta \right) $$

Where h 0(y) is the baseline hazard, estimated using the Breslow method. The main issue in this model, which should be considered, is the proportional hazard assumption. This means that the hazard ratio of two individuals is independent of time [14]. Note that the interpretation of the estimated coefficients in this model is based on hazard ratio rather than the covariate effect on the mean.

Comparing model performance

The interested estimations are evaluated as follows:

Two statistics were calculated to evaluate the quality of cost predictions using above mentioned models. The first was the mean prediction error (MPE), which measures the bias and predictive accuracy, and the second was the mean absolute prediction error (MAPE):

$$ MPE=\frac{1}{n}{\displaystyle \sum_{i=1}^n\left({y}_i-{\widehat{y}}_i\right)} $$
$$ MAPE=\frac{1}{n}{\displaystyle \sum_{i=1}^n\left|{y}_i-{\widehat{y}}_i\kern0.1em \right|} $$

Actually, MPE indicates how the mean of predicted healthcare expenditures from a particular model compares with the mean of healthcare costs. Models with lower values of MPE have smaller biases than models with higher values. However, MAPE indicates how values of individual predicted healthcare expenditures from a particular model compare with values of actual healthcare expenditures in the sample [6].

Mean square of error (MSE) and 95% confidence interval of the estimate of β1 coefficient were calculated to evaluate the accuracy and precision of the estimated parameter. A more precise estimator should be closer to the true value. A Goodness of fit test provided by Hosmer-Lemeshow test and the Akaike information criteria (AIC) used as an aid to choosing between competing models. Lower values of the AIC indicate the preferred model criterion were also used to evaluate. The mean of the residuals across deciles of x was also plotted in order to assess a systematic bias in the predictions.

Simulation study

To compare the performance of the alternative models, a Monte Carlo simulation was used to show how each estimator behaves under different conditions of skewness that are common in healthcare expenditure studies.

To determine the effect of the level of skewness on the estimated outcome, some skewed probability density function (pdf), such as log-normal, Gamma and Weibull distribution, was used as a data-generating mechanism. To assess how the sample size affects the estimations, 10,000 times batch samples (n = 25, 50, 100, 500 and 1,000) were examined by comparing all models mentioned in the previous section. All generated data were standardized according to Basu et al., in which β 0 was considered as intercept, estimated assuming E(y) = 1.

Log-normal data generation

The true model assumed is as follows:

$$ \ln (y)={\beta}_0+{\beta}_1+\varepsilon $$

Where x is uniform (0, 1), εN(0, σ 2), in which σ 2 = 0.5, 1.0, 1.5, and β 1 = 1 were used. β 0 was estimated based on E(y) =1:

$$ E\left(y\Big|\;x\right)= \exp \left({\beta}_0+{\beta}_1x+0.5{\sigma}^2\right) $$

The skewness of log-normal models is an increasing function of the variance as follows:

$$ \left( \exp \left({\sigma}^2\right)+2\right){\left( exp\left({\sigma}^2\right)-1\right)}^{0.5} $$

We considered σ 2 = 0.5, 1, 1.5 and 2.

Gamma data generation

The pdf of Gamma distribution is:

$$ f(y)=\frac{1}{\varGamma \left(\alpha \right){b}^{\alpha }}{y}^{\alpha -1}{e}^{-y/b} $$

Where b = exp(β 0 + β 1 x) and α are the scale and shape parameters, respectively. The mean is equal to αb and the skewness is a decreasing function of the shape parameter, as follows:

$$ \frac{2}{\sqrt{\alpha }} $$

Where x is uniform (0, 1), β 1 = 1 and β 0 was estimated so that E(y) = 1. Also, we used the assumption that α = 0.5, 1, 2 and 4 in the data generating process.

Weibull data generation

Weibull data generation is considered as a function of the data-generating mechanism, which has proportional hazard properties. It was used to generate proportional hazard data, since the Cox proportional hazards model is based on this assumption:

$$ f(y)=\frac{\alpha }{b}{\left(\frac{y}{b}\right)}^{\alpha +1}{e}^{{\left(-y/b\right)}^{\alpha }} $$

Where b = exp(β 0 + β 1 x) and α are the scale and shape parameters, respectively. The mean is equal to \( b\varGamma \left(1+\frac{1}{\alpha}\right) \) and the skewness is also a decreasing function of the shape parameter, as follows:

$$ {b}^3\varGamma \left(1+\frac{3}{\alpha}\right)-3\varGamma \left(1+\frac{1}{\alpha}\right)\varGamma \left(1+\frac{2}{\alpha}\right)+2{\left(\varGamma \left(1+\frac{1}{\alpha}\right)\right)}^3 $$

Shape parameter was considered as 0.5, 1 and 5 in this scenario. The proportional hazards assumption was evaluated in all of the simulations.

Results

Mean, standard deviation, skewness and kurtosis for the y in various data-generating mechanisms are presented in Table 1. Based on this result, the log-normal and Weibull models provided greater skewness than the Gamma model. It should be noted that the skewness of data from the log-normal and Gamma models increased monotonically as the sample size increased.

Table 1 Simple statistics of y

The results in Tables 2, 3, 4, 5 and 6 were based on 10,000 times batch replication, in sample sizes of 25, 50, 100, 500 and 1,000, respectively. These tables show the results of the estimates of population means and β1 for each model under the various data-generating processes. Minimum deviance (MPE) and absolute deviance (MAPE) of predicting the value of the response variable (health-care costs) considered as adequacy of methods.

Table 2 Alternative estimator results for log-normal, gamma and weibull distributions for n=25
Table 3 Alternative estimator results for log-normal, gamma and weibull distributions for n=50
Table 4 Alternative estimator results for log-normal, gamma and weibull distributions for n=100
Table 5 Alternative estimator results for log-normal, gamma and weibull distributions for n=500
Table 6 Alternative estimator results for log-normal, gamma and weibull distributions for n=1000

Generally, entire models exhibited lower MPE by declining skewness and increasing sample size. However, the Gamma regression model had the smallest biases across all data-generating processes. Moreover, our results indicated that its ability to predict the expenditures in a small sample size was as good as for large sample sizes. Furthermore, OLS for Ln(y) and Weibull regression models showed a lower bias than the Cox proportional hazard model, even in proportional hazard data-generating process (Figure 1).

Figure 1
figure 1

Mean residual from different estimators across deciles of ‘X’ for log-normal data (n=25) with variance a: 0.5, b: 1.0, c: 1.5, d: 2.0.

Figure 2
figure 2

Mean residual from different estimators across deciles of ‘X’ for Gamma data (n=25) with shape parameter a: 0.5, b: 1, c: 2, d: 4.

Figure 3
figure 3

Mean residual from different estimators across deciles of ‘X’ for Weibull data (n=25) with shape parameter a: 0.5, b: 1, c: 5.

Figure 4
figure 4

Mean residual from different estimators across deciles of ‘X’ for log-normal data (n=50) with variance a: 0.5, b: 1.0, c: 1.5, d: 2.0.

Figure 5
figure 5

Mean residual from different estimators across deciles of ‘X’ for Gamma data (n=50) with shape parameter a: 0.5, b: 1, c: 2, d: 4.

Figure 6
figure 6

Mean residual from different estimators across deciles of ‘X’ for Weibull data (n=50) with shape parameter a: 0.5, b: 1, c: 5.

Figure 7
figure 7

Mean residual from different estimators across deciles of ‘X’ for log-normal data (n=100) with variance a: 0.5, b: 1.0, c: 1.5, d: 2.0.

Figure 8
figure 8

Mean residual from different estimators across deciles of ‘X’ for Gamma data (n=100) with shape parameter a: 0.5, b: 1, c: 2, d: 4.

Figure 9
figure 9

Mean residual from different estimators across deciles of ‘X’ for Weibull data (n=100) with shape parameter a: 0.5, b: 1, c: 5.

Figure 10
figure 10

Mean residual from different estimators across deciles of ‘X’ for log-normal data (n=500) with variance a: 0.5, b: 1.0, c: 1.5, d: 2.0.

Figure 11
figure 11

Mean residual from different estimators across deciles of ‘X’ for Gamma data (n=500) with shape parameter a: 0.5, b: 1, c: 2, d: 4.

Figure 12
figure 12

Mean residual from different estimators across deciles of ‘X’ for Weibull data (n=500) with shape parameter a: 0.5, b: 1, c: 5.

Figure 13
figure 13

Mean residual from different estimators across deciles of ‘X’ for log-normal data (n=1000) with variance a: 0.5, b: 1.0, c: 1.5, d: 2.0.

Figure 14
figure 14

Mean residual from different estimators across deciles of ‘X’ for Gamma data (n=1000) with shape parameter a: 0.5, b: 1, c: 2, d: 4.

Figure 15
figure 15

Mean residual from different estimators across deciles of ‘X’ for Weibull data (n=1000) with shape parameter a: 0.5, b: 1, c: 5.

In addition, evaluating MAPE as an accuracy measure showed that Gamma and Weibull regression models have almost equal MAPE values. In many conditions, such as the log-normal model with σ2 = 1.5, 2, the Gamma model with shape equal to 0.5 (monotonically declining pdf) and the Weibull model with shape equal to 0.5 (linearly increasing hazard), as higher skewness data-generating mechanisms, the MAPE from Weibull regression model was always lower than Gamma regression model.

Interestingly, as the sample size increased, the MAPE of OLS for Ln(y) became very similar to that of the Gamma regression model. However, the MAPE of all models had an insignificant upward trend as the sample size increased.

Since there was also a concern about consistency and precision in the estimates of β1 coefficients, MSE and 95% simulation intervals were investigated. All three regression Gamma and Weibull and OLS for Ln(y) models provided approximately similar MSEs of β1 as data generated using log normal. However, the Gamma regression model showed minimum MSE values. We also found that MSE decreased by reducing skewness and increasing sample size. For the Weibull-generated data, Gamma and Weibull regression models exhibited similar and minimum values of MSE. Under all data-generating mechanisms, 95% simulation intervals were closer to true values in all three regression models. Surprisingly, the Cox proportional hazard model revealed maximum MSE and less accurate 95% simulation intervals, even within proportional hazards data-generating scenario.

Comparison goodness of fit tests (Hosmer-Lemeshow test and AIC criterion) revealed that, under a different range of data conditions, Gamma and Weibull regression models were better behaved. Finally, investigation of the pattern of the residuals as a function of X, which have been implemented by the mean of the residuals across deciles of X, showed more bias for the Cox proportional hazard model across all generated data and sample sizes (see Figures 2-15).

Discussion

Although there are many substantial studies addressing the statistical issues in healthcare cost analysis over the last few decades, it is still an important issue that needs further evaluation. In this paper, we assessed the performance of various well-known statistical regression-based models in healthcare expenditure analysis, through different sample sizes and data-generating processes, using a Monte Carlo simulation. Each model was evaluated on 10,000 batch random samples, with 25, 50, 100, 500 and 1,000 sample sizes. Other studies were performed using just one large sample size (such as 10,000) [5,10], while we know the sample size is an important issue in healthcare studies and in precision of model-based estimators.

We found that, by considering the different interest points of various research and various data conditions, different model-based estimators could be used. Indeed, no estimator is considered to be the best across all ranges of data-generating processes. In addition, the accuracy of the results was almost the same in all sample sizes.

However, the GLMs estimated population means more precisely in all data-generating processes and sample sizes. In this respect, our results are consistent with other studies [2,5,6,10]. Comparative studies between log models were evaluated on 1,000 random samples, with a sample size of 10,000. They found almost identical results in estimating the slope β1, but the GLMs were substantially more precise than OLS-based model [5]. In this paper, as the sample size increased, the precision of estimating the mean population and the β1 using an OLS-based model became closer to that of GLMs.

Based on our result, the Gamma regression model provided more accurate estimates of population mean. In other studies, which compare log and Cox proportional hazard models, the Gamma regression model was introduced as the reasonable model across all of the simulation processes [13]. They have also found that the Cox proportional hazard model exhibited good performance when data were generated by distribution with a proportional hazards assumption [13]. In this paper, a Weibull distribution was selected as the proportional hazard data-generating mechanism. In addition, investigating proportional hazards assumption detected that gamma generation process also has produced data with proportional hazard properties but the Cox proportional hazard model showed a poor result within these data generation process. We also found that the Cox proportional hazard model behaved poorly in other data generation scenarios.

Our study has some limitations, including the fact that our focus was on generating skewed data, while kurtosis may have affected the results. Furthermore, the study was limited to fixed covariates.

Conclusions

Selecting the best model is dependent on the interest point of research, which could be the estimated mean of the population or covariate effects. There is no best model among all data conditions. It seems that the GLMs, especially the Gamma regression model, behave well regarding the estimation of population means of healthcare costs in most of the conditions. The results are consistent among all sample sizes; however, increasing sample size leads to improvement in the performance of the OLS-based model.

Based on estimation of the β1, GLMs seems to provide plausible estimations and as the sample size increased, estimated the β1 more precisely in all data-generating processes. Under all data generation, process even proportional hazard data generation scenarios the Cox proportional hazard model provided a poor estimation of mean population and the β1.

Abbreviations

Prob. H.L.:

Hosmer–Lemeshow test

signif.:

At the 5% level

References

  1. Gregori D, Petrinco M, Bo S, Desideri A, Merletti F, Pagano E. Regression models for analyzing costs and their determinants in health care: an introductory review. Int J Qual Health Care. 2011;23(3):331–41.

    Article  PubMed  Google Scholar 

  2. Mihaylova B, Briggs A, O’Hagan A, Thompson SG. Review of statistical methods for analysis healthcare resources and costs. Health Econ. 2011;20:897–916.

    Article  PubMed Central  PubMed  Google Scholar 

  3. Gilleskie DB, Mroz TA. A flexible approach for estimating the effects of covariates on health expenditures. J Health Econ. 2004;23:391–418.

    Article  PubMed  Google Scholar 

  4. Cantoni E, Ronchetti E. A robust approach for skewed and heavy-tailed outcomes in the analysis of health care expenditures. J Health Econ. 2006;25(2):198–213.

    Article  PubMed  Google Scholar 

  5. Manninga WG, Mullahy J. Estimating log models: to transform or not to transform? J Health Econ. 2001;20:461–94.

    Article  Google Scholar 

  6. Deb P, Burgess JF. A quasi-experimental comparison of econometric models for health care expenditures. New York: Hunter College Department of Economics; 2003. Report No.: 212.

  7. Basu A, Manning WG. Issues for the next generation of health care cost analyses. Med Care. 2009;47(7):S109–14.

    Article  PubMed  Google Scholar 

  8. Ravangard R, Arab M, Rashidian A, Akbarisari A, Zare A, Zeraati H. Comparison of the results of cox proportional hazards model and parametric models in the study of length of stay in a tertiary teaching hospital in Tehran. Iran Acta Med Iran. 2011;49(10):650–8.

    Google Scholar 

  9. Jain AK, Strawderman RL. Flexible hazard regression modeling for medical cost data. Biostatistics. 2002;3(1):101–18.

    Article  PubMed  Google Scholar 

  10. Manning WG, Basu A, Mullahy J. Generalized modeling approaches to risk adjustment of skewed outcomes data. J Health Econ. 2005;24:465–88.

    Article  PubMed  Google Scholar 

  11. Faddy M, Graves N, Pettitt A. Modeling length of stay in hospital and other right skewed data: comparison of phase-type, gamma and log-normal distributions. Value Health. 2009;12(2):309–14.

    Article  PubMed  Google Scholar 

  12. Griswold M, Parmigiani G, Potosky A, Lipscomb J. Analyzing health care costs: a comparison of statistical methods motivated by Medicare colorectal cancer charges. Biostatistics. 2004;1(1):1–23.

    Google Scholar 

  13. Basu A, Manning WG, Mullahy J. Comparing alternative models: log vs Cox proportional hazard? Health Econ. 2004;13:749–65.

    Article  PubMed  Google Scholar 

  14. Lee ET, Wang JW. Statistical methods for survival data analysis. 3rd ed. New Jersey: John Wiley & Sons; 2003.

    Book  Google Scholar 

Download references

Acknowledgements

This study is part of biostatistics MS degree thesis of Fatemeh Pourmotahari and it was supported by the Ahvaz Jundishapur University of Medical Sciences.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Amal Saki Malehi.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

ASM contributed to the study design, wrote and revised the manuscript. FP analyzed the data and drafted the manuscript. KAA revised the manuscript. All authors read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0), which permits use, duplication, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Malehi, A.S., Pourmotahari, F. & Angali, K.A. Statistical models for the analysis of skewed healthcare cost data: a simulation study. Health Econ Rev 5, 11 (2015). https://doi.org/10.1186/s13561-015-0045-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13561-015-0045-7

Keywords