1 Introduction

Health care spending has grown faster than other sectors of the United States economy for decades. Total healthcare expenditures accounted for 17 % of Gross Domestic Product (GDP) in 2015, a percentage far exceeding all other Organization for Economic Co-operation and Development (OECD) countries. This percentage is projected to continue rising over the next decade [1].

Given the prominence of healthcare in the U.S. economy, predictive models of medical costs, especially mean medical costs, are of great interest to policy makers. But medical costs are often right-skewed and heteroscedastic, making statistical analysis challenging. Logarithmic transformation has been often proposed, but requires modeling on the transformed scale, and heteroscedasticity may bias estimates derived from logged response variable models [2, 3]. It is therefore desirable to estimate logE(Y) instead of E(logY) to obtain the mean medical cost easily, and to consider robust statistical methods by avoiding unnecessary assumptions. For example, Chen, Liu, Zhang, and Shih [4] and Chen, Liu, Zhang, Shih, and Severini [5] developed models using log link and a fixed set of predictors, under assumptions only on the mean cost without additional restrictive assumptions on higher order moments.

Variable selection is critically important for cost prediction models because policy makers wish to know which variables “drive” medical costs, and favor simple models that are easy to interpret and estimate - usually models with fewer predictors.

Statistical methods of variable selections can be performed using either Bayesian or frequentist approaches. Bayesian approaches are more informative than frequentist counterparts, as they can provide posterior probability ratios for selected models, stating that the top model A is five times more likely to be the true model than the next model B, etc., which cannot be done by frequentist approaches. While the Bayesian approach of variable selection is appealing, it is computationally challenging. Specifically, Bayesian model selection methods such as BIC or Bayes Factor require assumptions on the true likelihood function. Since the posterior in a Bayesian approach is proportional to the product of the likelihood and the prior, it is seemingly impossible to conduct Bayesian analysis for robust methods without likelihood functions. To resolve this problem, we use an asymptotic likelihood based on the asymptotic normality of a parameter estimate instead of the original data and apply the robust sandwich estimate of the asymptotic variance, so that Bayesian inference based on this asymptotic normal likelihood function is valid under the assumption of the mean model. We employ a “quasi-Bayesian” approach from the econometrics literature (see, e. g. Chernozhukov and Hong [6], Inoue and Shintani [7]), that is less known in statistical applications. When multiplied by a prior density that either proposes a “spike” or a “slab” for the parameter components, we can form such an asymptotic posterior and perform Bayesian variable selection.

The method developed in this paper has four goals: (i) robustness to skewness and heteroscedasticity with only one assumption on the mean model; (ii) parsimonious variable selection; (iii) generation of posterior probabilities of candidate models; (iv) computational efficiency by avoiding either numerical or Monte Carlo integration, and limiting model search efforts. Our proposed method simplifies search for the model with the highest posterior probability and achieve computing efficiency based on ranking the Z-statistics, following Jiang and Liu [8]. In addition, with the Spike-or-Slab priors, the integration of the posterior computation is exact and no iteration is needed. Although recent work by Li and Jiang [9] addressed these three goals (i), (ii) and (iii) simultaneously, their method is computationally cumbersome because it involves a much more complicated posterior function that needs to be simulated by Markov chain Monte Carlo (MCMC) methods, which is no longer necessary in our model. Jiang and Liu [8] used the same idea of posterior inference based on parameter estimates as in this paper, but they did not do exact integration for the posterior probability of each candidate model as what we will do here. In addition, they did not use spike or slab prior densities to describe the prior belief of either small or large regression coefficients.

LASSO (Least Absolute Shrinkage and Selection Operator) and sslasso other commonly used variable selection methods, but neither achieves all our goals. LASSO is a frequentist approach and thus does not provide posterior probabilities of the candidate models. Ročková and George [10] proposed the spike-and-slab lasso and Tang, Shen, Zhang, and Yi [11] extended it to generalized linear models, termed sslasso, which could be applied to model the medical costs nonlinearly against the regressors. Their method is implemented in R package BhGLM. sslasso, as our method, is a Bayesian method with spike and slab priors. However, it does not use a robust estimate of the asymptotic variance, and therefore cannot guarantee that the reported posterior probabilities are asymptotically valid from the Bayesian perspective, when only the mean model is correctly specified. We will compare our method empirically to these alternative variable selection methods.

Our motivating example uses the Medical Expenditure Panel Survey (MEPS) dataset. We include 3,376 individuals aged over 65 from the MEPS 2014 Full Year Consolidated Data File. The mean medical cost in 2014 is $10,321 and the median is $4,394.50, illustrative of the right skewness. As a large scale survey of health, health services use and associated costs of noninstitutionalized civilians in the United States, the MEPS included many variables, and is thus well suited example on which to compare the abilities of variable selection methods to develop parsimonious cost prediction models.

The paper is organized as follows. Details about our proposed method are provided in Section 2. We compare our method’s estimation performance to other alternatives in Section 3. In Section 4, we demonstrate our method on the MEPS dataset. We discuss the results and draw conclusions in Section 5.

2 Methods

3 Simulation study

We conduct simulation studies to assess the performance of our proposed method (SorS). The response Y are generated from Gamma distribution where Gamma(a,b) density function is f(y)=ba/Γ(a)ya−1e−by. In our application study, we have many binary variables as covariates, so we only consider binary variables for X in the simulation. Correlated binary variables are generated using R package bindata. In all cases, we generated R = 100 datasets, each with sample size n = 1000 and p = 50 where μi=eXiβ. The true parameter coefficient has an intercept and four nonzero components (p0=4) as below:

1. Yi∼ Gamma(μi,1), that is, E(Yi)=μi, V(Yi)=μi.

   1. (Independent predictors) x1,…,xp’s are iid from Bernoulli(0.5).

   2. (Correlated predictors) x1,…,xp∼ Bernoulli(0.5) with corr(xj,xk)=0.5|j−k| for 1 ≤ j,k ≤ p.

2. Yi∼ Gamma(1,1/μi), that is, E(Yi)=μi, V(Yi)=μi2.

   1. (Independent predictors) x1,…,xp’s are iid from Bernoulli(0.5).

   2. (Correlated predictors) x1,…,xp∼ Bernoulli(0.5) with corr(xj,xk)=0.5|j−k| for 1 ≤ j,k ≤ p.

We compare SorS with three different methods: sslasso, LASSO and full model without variable selection. For estimation, we would like to compare RMSE =∑i=1n(Yi∗−Yˆ(Xi∗))2/n where Yˆ(Xi∗)=eXi∗βˆ and the test data (Xi∗,Yi∗), i = 1,…,n = 1000 are generated independently from the same distribution as the training set. In addition, we will also present the following criteria to evaluate the performance of these methods. The comparison results are shown in the following figures and table.

1. RMSE of coefficient estimates (bRMSE) =∑j=0pβj−βˆj2/(p+1).

2. Selected model size.

3. In the table,

   1. Coverage probability (Cov)=∑r=1RI(M∗⊂Mˆ(r))/R,

   2. True negative rate (TNR)=∑r=1R∑j=1pI(βˆj(r)=0,βj=0)/[R(p−p0)],

   3. False negative rate (FNR)=∑r=1R∑j=1pI(βˆj(r)=0,βj≠0)/[Rp0],

   4. Exact selection probability (Ext)=∑r=1RI(M∗=Mˆ(r))/R,

      where M∗ and Mˆ(r) denote a true model and a selected model at r th generated dataset, respectively.

1. Average accuracy rate of variable selection (Acc) =∑j=1pI(γˆj=γj)/p, where γj=I(βj≠0) and γˆj=I(βˆj≠0).

In Case 1 with moderate heteroscedasticity, variance is the same as mean. The median model size selected by LASSO is 15 as shown in Figure 1A, which implies a lot of false positives. For the case without heteroscedasticity (not shown here), LASSO works well. As is well known, LASSO is very fast and works better with independent predictors. The estimated coefficient βˆ from LASSO is re-fitted without shrinkage using only selected variables. Although sslasso is fast and still works well in the correlated setting, the weakness of sslasso is its incapability to deal with heteroscedasticity. As heteroscedasticity becomes more severe in Case 2 where variance is the square of the mean, SorS shows its advantage. Not only has SorS the smallest RMSE in Figure 2, but also SorS is the most parsimonious method and includes only a few false positives as seen in Figure 1. In Case 2 when heteroscedasticity is severe, its average accuracy is the highest as 0.97 and still can identify the true model more than half times in both independent and correlated covariate settings as shown in Table 1. In summary, SorS performs better than the other comparable methods in terms of selection (Figure 1), prediction (Figure 2) and estimation performance (see Figure 3). All method is implemented in R and the source code is available in https://github.com/GraceYoon/SorS.

Figure 1:

Boxplots of selected model sizes over 100 replications. The horizontal lines indicate the true model size, which is 4.

Figure 2:

Boxplots illustrating RMSE for test data in all cases.

Figure 3:

RMSE of coefficient estimates in each case from 100 replications.

Our proposed method is based on estimating a Poisson estimating equation with the log-linear mean specification. Our method is robust against violation of the Poisson-type variance assumption, as is shown in our simulation that allows the variance to be unequal to the mean. As is typical in the biostatistical applications of the estimating equations, the proposed method is NOT intended to be robust against mean specification. For example, for a true model Y=exp(X2) on X∼Unif(−1,1), fitting a Poisson model with log-linear mean in X would obtain a constant fit, suggesting (incorrectly) that E(Y|X)=constant, independent of X.

Specifically, for medical cost data, generalized linear models (GLM) with log link have been used most often to describe the skewed and heteroscedastic medical cost data in the original scale [3, 14, 15, 16, 17, 18]. In this paper, we have not considered other type of link functions. However, more flexible link functions, e. g. Box-Cox transformation [19], could be considered in the future work. Another option is to add nonlinear functions for covariate effects, e. g. through polynomial functions or splines [4], to accommodate the curvature and account for mean misspecification.

4 Application to MEPS data

We construct an analytical sample from the 2014 MEPS Full Year Consolidated Data File, including every subject ≥ 65 years old in the survey year, and extract a total of 33 variables to account for the impact of each individual’s disease(s), symptom(s) or any related conditions on medical costs. These variables are listed in Appendix B. Most are self-explanatory. EDUCAT, the highest degree of education in 2014, is coded as an integer from 1 to 3, for < 13 years, 13–16 years, or > 16 years, respectively. POVCAT is coded as: 1: poor, 2: nearly poor, 3: low income, 4: middle income, 5: high income. Other continuous variables are standardized in our analyses. Correlations of covariates range from −0.55 to 0.6. The final study sample includes 3,376 individuals. As noted, the response variable in our model, annual medical costs in U.S. dollars, is highly right-skewed; 3.6 % of subjects had no medical cost, which our method accommodates since the naive Poisson likelihood supports 0 values.

We apply Spike-or-Slab normal priors for the parameters a uniform prior for model selection. Posterior probabilities are obtained from the exact integration as described in Section 2.2. To assess the performance of various variable selection methods, we randomly sample half the observations to form a training set and use the remaining as test data to compute the RMSE, repeating this 100 times. The results are summarized in Figure 4. RMSE values for SorS method are close to those from other methods, but selected model sizes are less than half of others.

Figure 4:

Prediction performance and selected model size by different methods in MEPS data analysis.

Applying the SorS method to all the data, the model with 10 variables has the largest posterior probability. The candidate model with 11 variables has the second highest posterior probability and is 21.21 % as likely as the model with 10 variables. Using 10 variables in the highest posterior probability model, the estimated coefficients are shown in Table 2.

Results from the best performing model (with 10 variables) show that hospitalization (HOSPEXP) and emergency room visit (EMERG) increase the medical cost by a large percentage (3.4 times and 1.3 times, respectively). Heart and blood vessel disease (CORHRT and STRK), body movement disorder (ANYLMT), cancer and diabetes are all significantly associated with annual medical costs. Gender and race variables are not selected in our model.

Expected medical costs are higher for the higher educated: more educated individuals may be more likely to have regular checkups and better able and prone to spend more for treatment. Physical Composite Scores (PCS), a quality of life measure, is inversely associated with medical costs.

5 Discussion

In this work, we propose a new quasi-Bayesian variable selection method with Spike-or-Slab (SorS) priors to simultaneously accommodate severe skewness and heteroscedasticity of the medical cost data without any transformation. We integrate four features: log-linear mean model; spike and slab prior; ranking by Z-statistics; and posterior computation without MCMC. The novelty of our work lies in the innovative integration and application to skewed and heteroscedastic data. Our proposed method is applicable to many situations where the likelihood function is not available. In addition, our method uses the robust sandwich formula for constructing the posterior, and is therefore more robust in the sense that the asymptotic validity of the reported posterior probabilities depends only on the correct specification of the mean model without any additional assumptions on the variance structure or density function.

This method however is limited to low-dimensional data. Our asymptotic likelihood is based on asymptotic normality of the parameter estimates which relies on the standard large sample theory and ranking based on Z-statistics method cannot be implemented when the dimension of covariates is higher than the sample size. In high-dimensional case, sslasso works and is computationally feasible, but it performs poorly when the heteroscedasticity level is severe. It will be an interesting area of future research to consider robust Bayes variable selection for higher dimensional data with p > n.