3.1 Missing Observations
About 2% of quantitative observations were missing, so simple imputation using the mean was employed. This is conservative, as it tends to hide results that might be statistically relevant by reinforcing mean values. For the categorical variables, all but ownership were fully complete. There were only 14 missing observations for this variable, and these were imputed with the mode.
3.2. Descriptive Statistics-Quantitative Data
Descriptive statistics for the quantitative data are provided in Table 2. The average hospital observation during any given year had 1,635.62 observations of DRG 291, 292, and 293 (median of 383). That same average hospital had 146.39 staffed beds (median of 87), 6,996.58 discharges (median of 2,825), 6,348.76 surgeries (median of 4,490), and 34,181.38 acute days (median of 14,051). The average hospital had positive income (in millions) of $17.232 (median of $2.044), significant cash-on-hand ($20.28, median of $1.99), and positive equity. The typical hospital had 1005.53 employees (median of 437) with 231.37 affiliated physicians (median of 104) and was reimbursed 45% by Medicare (median of 42%). Only 9% reimbursement was from Medicaid (median of 6%).
Table 2. Descriptive statistics for the study (dollars in millions)
n=40,523 hospital observations
|
Mean
|
SD
|
Median
|
Min
|
Max
|
Quantitative Technical Variables
|
Number DRGs
|
1,635.62
|
3330.63
|
383.00
|
11
|
57,461
|
Staffed Beds
|
146.39
|
171.98
|
87
|
2
|
2,753
|
Affiliated Physicians
|
231.37
|
352.66
|
104
|
1
|
4328
|
Employees
|
1,005.53
|
1679.14
|
437
|
4
|
26,491
|
Percent Medicare
|
0.45
|
0.19
|
0.42
|
0
|
.98
|
Percent Medicaid
|
0.09
|
.09
|
.06
|
0
|
.87
|
Workload Variables
|
Discharges
|
6,996.58
|
9,881.39
|
2,825
|
1
|
129,339
|
ER Visits
|
32,865.58
|
33,893.77
|
25,236
|
0
|
543,457
|
Surgeries
|
6,348.76
|
7,965.82
|
4,490
|
0
|
130,741
|
Acute Days
|
34,181.38
|
51,725.31
|
14,051
|
5
|
701,074
|
Financial Variables
|
Net Income ($ in M)
|
$17.23
|
$117.65
|
$2.04
|
-$1.21
|
$3,31
|
Cash on Hand ($ in M)
|
$20.28
|
$120.24
|
$1.99
|
-$2.51
|
$3.88
|
Profit Margin
|
-0.03
|
1.25
|
-0.02
|
-15.45
|
62.07
|
Equity ($ in M)
|
$174.11
|
$625.76
|
$33.94
|
-$3.25
|
$10.24
|
Year over year, both DRGs and rates of DRGs per 1000 population increased as illustrated in Figure 2. The significance of the DRG increase is the financial consideration. The significance of the rate of DRG increase is the epidemiological consideration. If the DRG rate is considered a proxy for incidence rate, then there is either a significant increase, a coding issue, or something else. These considerations are found in the discussion section. One might expect the DRG rate graph to remain horizontal (static). Independent variables remained relatively constant year-over-year likely due to repeated measures on the same facilities.
Figure 2. Number and rates of DRGs as a function of year
3.3. Descriptive Statistics-Categorical Data
California, Texas, and Florida had the largest number of diagnoses for all years and year-over-year, largely due to population size, with averages of {1,669,210; 1,631,021; 1,490,983}, respectively. When adjusted per 1000 population, the District of Columbia, West Virginia, and Delaware dominated the with total rates per 1,000 population of {109.47, 102.86, 94.15}, respectively. Utah, Hawaii, and Colorado had the smallest average rates, {26.17, 28.87, 34.74}, respectively. Appendix Aillustrates the rates by state / territory.
Of the hospital observations, 6700 were rural (42%) while 9279 were urban (58%). Most of the hospitals (8342 or 52%) were voluntary non-profits with 29% (4641) proprietary and 18.7% (2996) governmental. The vast majority (11,914 or 75%) had no affiliation with a medical school and were short-term care facilities (9,604 or 60%). Nearly no hospitals were classified as Department of Defense (DoD) or children’s hospitals. Figure 3 depicts the categorical breakout by year.
Figure 3. Categorical variables by year
3.4. Descriptive Statistics-Financial Estimates
In FY 2008, the Centers for Medicare and Medicaid (CMS) estimated that heart failure DRGs 291, 292, and 293 national average total costs per case were {$10.235, $6.882, $5.038} thousand, respectively. By FY 2012, CMS increased those estimates to {$11.437, $7.841, $5.400} thousand, respectively. In four years, the accumulation rates (1 plus the inflation rate) were 1.139, 1.117, and 1.072 for the DRGs in ascending order. Using these accumulation rates, estimates for 2016, 2017, and 2018 were generated. Table 3 shows these extrapolated estimates.
Table 3. Estimated total costs for heart failure by DRG in thousands, linear extrapolation method
DRG
|
2016
|
2017
|
2018
|
DRG 291
|
$12,780
|
$13,155
|
$13,243
|
DRG 292
|
$8,934
|
$9,245
|
$9,257
|
DRG 293
|
$5,788
|
$5,891
|
$5,998
|
Another method for estimating these costs involved the use of the Federal Reserve Bank of Saint Louis (FRED) producer price index for general medical and surgical hospitals [30]. The annual accumulation rates for 2013 through 2018 were estimated as {1.022, 1.012, 1,007, 1.013, 1.018, 1.023}, respectively. Applying these to the 2012 total costs from CMS results in Table 4 estimates for 2016 through 2018.
Table 4. Estimated total costs for heart failure by DRG in thousands, medical inflation rate method
DRG
|
2016
|
2017
|
2018
|
DRG 291
|
$12,058
|
$12,273
|
$12,582
|
DRG 292
|
$8,267
|
$8,414
|
$8,626
|
DRG 293
|
$5,693
|
$5,795
|
$5,491
|
Both estimates are fairly close. To estimate costs, we used both of these tables separately as upper and lower bounds. Since these total costs represent only CMS costs, the actual financial burden across all payers is likely underestimated as commercial third-party insurers can reimburse up to 90% more than Medicare for the same diagnosis [31]. Figure 4 illustrates the number of DRGs by year, while Figure 5 shows the associated aggregate cost estimates.
Figure 4. Number of DRGs by type (left axis) and cost estimates by DRG type and total, 2016 through 2018
Figure 5. Associated cost estimates in billions (total and by DRG) per year
In Figure 4, it is clear that DRG 291, the DRG with the highest average reimbursement rate per case, has increased nonlinearly, while DRG2 292 has seen a small drop, and DRG 293 is flat. In Figure 5, the total cost estimates for 2018 are nearly $66 billion more than 2016 on average. DRG 291, the most expensive DRG, has seen reimbursement increases of $92 billion on average. Reasons for such an increase are explored in the discussion section.
3.4. Descriptive Statistics-Correlational Analysis
Hierarchical clustered correlation analysis of quantitative variables (Figure 6) illustrate tight relationships among many variables. Hierarchical clustered correlation analysis clusters variables based on distance measures (e.g., Euclidean), so that those which are most highly correlated are close in location. These variables are then placed into a correlation plot or correlogram. Figure 7 illustrates that discharges, acute days, and staffed beds are most closely associated with the number of diagnoses, our primary variable of interest.
Figure 6. Hierarchical clustered correlation of quantitative variables
Analysis of the relationship between some categorical variables and the number of diagnoses also proved interesting. Notched boxplots by year and medical school affiliation reveal that a major affiliation experiences a larger number of diagnoses at the .05 level a result that is to be expected. (See Figure 6). Further, voluntary not-for-profits see a larger number of diagnoses (Figure 8).
Figure 7. Number of diagnoses by year by medical school affiliation
Figure 8. Number of diagnoses by year by type of hospital
3.5 Geospatial Analysis
A descriptive analysis of heart failure over time using geographical informat systems was conducted to evaluate regional differences. Primarily, we were interested in rates per 1000 in the population. Populations over time were based on Census Bureau estimates [5].
While DRG rates per 1000 were not constant over time, the concentrations were fairly consistent. There is a clear bifurcation in the center of the United States separating high and low rates. That bifurcation suggests a clear West-East difference, favoring the West Coast. Washington, D.C. has had (on average) the highest diagnoses of heart failure followed by West Virginia, Alabama, Mississippi, Michigan, Louisiana, Kentucky, and North Dakota. Of interest is that previous studies indicate these states are also plagued by the opioid crisis [23]. Figures 9, 10, and 11show the diagnoses per 1000 by year 2016, 2017, and 2018, respectively.
Figure 9. DRG rates per 1000, 2016
Figure 10. DRG rates per 1000, 2017
Figure 11. DRG rates per 1000, 2018
From 2016 through 2018, the average rate of diagnoses per 1,000 population increased for nearly all states. A Friedman rank sum test (paired, non-parametric ANOVA) of rates by state by year revealed significantly different rates by year by state (c22=70.941, p<.001). Figure 12 illustrates the changes by year and by state.
Figure 12. Diagnoses per 1,000 by year by state
When 2018 data are aggregated at the state / territory level, the DRGs per 1000 paint a slightly different picture, with high-intensities in Washington D.C., West Virginia, Delaware, Mississippi, Kentucky, North Dakota, Michigan, and Missouri (listed in descending order.) Further, superimposing obesity prevalence intensity from the Centers for Disease Control and Prevention (CDC) shows significant correlation between obesity and DRGs per 1000 [32]. A Spearman’s test for correlation of obesity prevalence and 2018 DRGs per 1000 was statistically significant with rho=.689, S=6,867.7, p<.001. Figure 13 provides the map of obesity and DRG rates.
Figure 13. Map of DRG Rates / 1000 versus obesity prevalence
3.6. Explanatory Models
Explanatory models were sought to explain the number of diagnoses by facility. The importance of these models is that we might estimate demand based on workload, technical, financial, and geospatial-temporal variables. A discussion of data preparation and analysis follows.
3.6.1. Box-Cox Multivariate Transformations
To meet required regression assumptions, multivariate transformation using Box-Cox methods was conducted on location-transformed variables. The location transform was necessary to ensure that all variables were positive definite. Box-Cox methods search for the optimal power transform of all variables simultaneously such that the assumption of multivariate normality cannot be rejected. A logarithmic transform is defined as the power of zero. In order to ensure that all possible transformations are feasible, the data must be positive definite. Thus each variable that was non-positive definite had the absolute value of the minumum added to each observations plus .01. Doing so ensured a positive definite location transform. These transformations are necessary only for non-tree models. (Tree models are location-scale invariant.) To prevent bias from being induced into the unknown test set, the transformations are completed only on the training set. The optimal powers found from the optimization associated with Box-Cox multivariate analysis are then applied to the test set. See Table 4 for the optimal powers.
Table 4. Optimal power transformations, Box-Cox methods
Variable
|
Power
|
Diagnoses
|
-.0673
|
Discharges
|
0.290
|
ER_Visits
|
0.417
|
Surgeries
|
0.364
|
Acute_Days
|
0.238
|
Net_Income
|
0.541
|
Operating_Profit_Margin
|
0.530
|
Cash_on_Hand
|
0.919
|
Equity
|
0.344
|
Staffed_Beds
|
0.170
|
Affiliated_Physicians
|
0.226
|
Employees
|
0.086
|
Medicare
|
0.788
|
Medicaid
|
0.032
|
3.6.2. Regression Models
Using the postive definite, Box-Cox transformed data, a regression model was fit hierarchially using the following blocks (in order): technical, workload, financial, geo-spatial. The multivariate transformation assumes that at least some independent variables cannot be fully observed or that we have incomplete observations on variables that might be fully observed. Thus, the transformations from the Box-Cox methods attempt to achieve multivariate normality rather than univariate normality. Hierarchical models attempt to fit obvious (known) variable blocks first followed by those of mmost interest. In our case, all blocks were statistically relevant to the analysis (see Table 5).
Table 5. Hierarchical analysis suggests all blocks are important
Block
|
Res.Df
|
RSS
|
Df
|
Sum of Sq
|
F
|
Pr(>F)
|
Technical
|
32404
|
53.933
|
|
|
|
|
Workload
|
32400
|
44.398
|
4
|
9.595
|
1853.223
|
<.0001
|
Financial
|
32396
|
43.882
|
4
|
0.516
|
99.609
|
<.0001
|
Geospatial
|
32343
|
41.864
|
53
|
2.018
|
29.416
|
<.0001
|
Linear regression on the training set resulted in a reasonable fit that accounted for R2=.750 or 75% of the sum of squared variability. No collinearity problems were present after transformation. Performance on the training set was insightful; however, the proof of model explanatory power rests in the training set estimates of the test set values. Applying the parameter estimates generated from the training set to the test set resulted in an R2 of .749, barely any loss.
Given the model’s ability to predict, the linear regression model was re-run on the entirety of the dataset after re-estimating the Box-Cost transformations, transformations which were only slight different in magnitude than those produced by the training set. The results again produce R2=0.749. The actual versus predicted plot is shown in Figure 14.
Figure 14. Plot of actual test-set data versus predictions from the training set
Further, we evaluated the coefficients and directions of those coefficients for forecasting the transformed dependent variable based on the number of variables included in the model. The top 1- variables in the regression model with their associated R2 are shown in Table 6. Discharges, Medicare percentage, and hospital type are the primary variables of interest.
Table 6. Most important variables and associated parameter estimates based on number of variables in the model
# Variables:
|
1
|
2
|
3
|
4
|
5
|
6
|
7
|
8
|
Intercept
|
0.790
|
0.775
|
0.762
|
0.760
|
0.809
|
0.831
|
0.831
|
0.831
|
Discharges
|
-0.011
|
-0.011
|
-0.011
|
-0.009
|
-0.010
|
-0.010
|
-0.010
|
-0.010
|
DRG 293
|
|
0.056
|
0.069
|
0.070
|
0.070
|
0.071
|
0.071
|
0.071
|
DRG 292
|
|
|
0.026
|
0.026
|
0.026
|
0.027
|
0.027
|
0.027
|
Acute-Care Hospital
|
|
|
|
-0.025
|
-0.033
|
-0.051
|
-0.050
|
-0.048
|
Medicare
|
|
|
|
|
-0.065
|
-0.068
|
-0.066
|
-0.065
|
Critical Access
|
|
|
|
|
|
-0.025
|
-0.025
|
-0.025
|
Major Med. School
|
|
|
|
|
|
|
0.014
|
0.013
|
Voluntary Non-Profit
|
|
|
|
|
|
|
|
-0.006
|
R^2
|
0.545
|
0.680
|
0.703
|
0.719
|
0.733
|
0.741
|
0.744
|
0.745
|
Outside of simple linear regression, we explored constrained regression techniques (lasso, ridge, and Elastic Net). Lasso regression was inferior to linear regression in terms of R2 (.651 vs. .750) on the test set. Ridge regression and Elastic Net were also unable to beat linear regression in terms of effects with both R2 nearly identical to that of lasso, {.651, .681}, respectively. In this case, the linear regression model was not overfit.
3.6.3. Tree Ensemble Models
Several tree regressor models were built and compared on an 80% training set. There was no need to use the transformed data for tree models, as these are location / scale invariant. These tree models included a bagging regressor (BR), a random forest regressor (RFR), an extra trees regressor (ETR), a gradient boosted regressor (GBR), and an extreme-gradient boosted model (XGR). Tree models are atheoretic, as each tree developed may be different from the previous one. When ensembled, variable importances emerge that determine which items are most important to determining how to classify or regress in a nonlinear fashion (piecewise). The number of trees used for each estimator was tuned along with the maximum depth of the trees (number of branches). A pseudo-random number ensured that any model improvements were not due to the random number stream. The results of these models on the unseen test set are shown in Table 6. Most importantly, all of these models account for more variance than regression models. The models predict at 97.1% and above in terms of variability capture. (See Table 7.)
Table 7. Coefficients of determination for the five tree ensemble models
Model
|
R2
|
Extreme Gradient Boosting
|
0.975
|
Extra Trees
|
0.975
|
Gradient Boosting
|
0.974
|
Random Forest
|
0.971
|
Bagging
|
0.970
|
Because of the tight congruence of these models, we ensembled the estimates of the number of DRGs forecast by each to produce importance statistics. The variables of most importance includes discharges, Medicare percentage, acute days, affiliated physicians, staffed beds, employees psychiatric hospital status, ER visits, medical school affiliation status, Puerto Rico status and surgeries. See Table 8.
Table 8. The variables and the importance factors associated with them, averaged over all tree ensemble models
Variable
|
Importance
|
Discharges
|
0.331
|
DRG__DRG291
|
0.287
|
Medicare
|
0.063
|
Acute_Days
|
0.035
|
Year__Y16
|
0.031
|
ER_Visits
|
0.029
|
Hospital_Type__Short Term Acute Care Hospital
|
0.026
|
State__VA
|
0.023
|
Affiliated_Physicians
|
0.023
|
State__MD
|
0.017
|
Hospital_Type__Rehabilitation Hospital
|
0.017
|
Staffed_Beds
|
0.016
|
State__FL
|
0.015
|
Surgeries
|
0.014
|
DRG__DRG292
|
0.011
|
State__OK
|
0.011
|
Medical_School_Affiliation__Major
|
0.010
|
State__IL
|
0.010
|
State__NE
|
0.010
|
State__OH
|
0.010
|
DRG__DRG293
|
0.010
|
Employees
|
0.009
|
State__NC
|
0.009
|
State__MI
|
0.009
|
Equity
|
0.008
|
Medicaid
|
0.007
|
Cash_on_Hand
|
0.007
|
Profit_Margin
|
0.006
|
Net_Income
|
0.006
|
State__CT
|
0.003
|
Year__Y18
|
0.002
|
Year__Y17
|
0.002
|
When comparing the regression models with the ensembled forests, we see that the first two terms are congruent (discharges and Medicare percent). Interestingly, no financial models are in the top 10 effect sizes of the regression or tree models. Facility technical and workload variables are the most important determinants of heart failure. In the tree models, there were piecewise linear effects identified for states that were not seen in the regression models.