WA Hospital Discharge Overdose 2010-2016

Assessing The ICD-9-CM/ICD-10-CM Transition with Time Series Analysis and Forecast.

WA Hospital Discharge Count trends

R packages used in this analysis.

See their Github or CRAN pages for further details

library(tidyverse)
library(tidyquant)
library(timetk)
library(sweep)
library(lubridate)
library(tibbletime)
library(plotly)
library(ggthemes)
library(forecast)
library(DT)

Table of the monthly counts of drug overdose by type of opioid. A few counts of less than 10 deaths have been suppressed from the table.

The data were fitted with a locally weighted smoothing (loess)

gg_ts <- drug_hosp_2010_16_m %>%
  ggplot(aes(x = dis_date, y = count, 
             color = drug_type)) +
  geom_point() +
  # geom_line() +
  geom_smooth(method = "loess") + 
  labs(title = "overdose monthly", x = "", 
       y = "od count") +
  facet_wrap(~ drug_type, nrow = 2, scale = "free_y") +
  expand_limits(y = 0) + 
  scale_color_tq() +
  theme_tufte() +
  theme(legend.position="none",
        axis.text = element_text(size = 15),
        strip.text = element_text(size = 20, 
                                  face = "bold"),
        panel.spacing = unit(1.5, "lines"))
gg_ts

Overdose Hospitalization Monthly Trends. Hovering over the graph provides more information

Time Series Decomposition

Preparing the data for time series analysis

Nesting the data by drug type

drug_hosp_2010_16_nest <- drug_hosp_2010_16_a %>%
  group_by(drug_type) %>%
  nest(.key = "drug_tab")

Create a time series objects for each drug type

drug_hosp_2010_16_ts <- drug_hosp_2010_16_nest %>%
  mutate(drug_ts = map(.x       = drug_tab, 
                       .f       = tk_ts, 
                       select = -dis_date,
                       start    = 2010,
                       freq     = 12))

decomposing and tidying the nested time series objects of each drug type

drug_hosp_2010_16_clas <- drug_hosp_2010_16_ts %>%
  mutate(drug_clas = map(.x = drug_ts,
                         .f = decompose,
                         type = "multiplicative"))

drug_hosp_2010_16_decomp <- drug_hosp_2010_16_clas %>%
  mutate(tidy_tab = map(drug_clas, sw_tidy_decomp)) %>%
  unnest(tidy_tab, .drop = TRUE)

Any drug

Additive decomposition

Recommended when the seasonal variation seems constant over time.

drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type) %>% 
  decompose() %>% autoplot()

Multiplicative decomposition

More appropriate when the seasonal variation increases over time.

drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type) %>% 
  decompose(type = "multiplicative") %>% autoplot()

Seasonal and Trend decomposition using Loess(STL) decomposition

it is more robust to outliers

drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>% select(n) %>% 
  unlist() %>% 
  ts(., start=2010, frequency = 12) %>% 
  stl(s.window = "periodic", robust = TRUE) %>% 
  autoplot()

The rest of the analysis will be based on the multiplicative assumption.

Any opioid

Any opioid (interactive version)

non_heroin_opioid

heroin

tidy decompose

Forecasting

See The forecast toolbox for more details.

In the graphs below the observed data are the black lines and the forcasted data are in blue with prediction intervals of 80% and 95% .

Any drug

drug_fc <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type)

drug_fc2 <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>%
  tk_ts(., start=2010, frequency = 12, 
        end =c(2015, 9), select = -drug_type)



 drug_fit2 <- snaive(drug_fc2,h=15)
 
 autoplot(drug_fc) +
 autolayer(drug_fit2, series="Seasonal naïve", 
           PI=TRUE, col = "skyblue", size = 2, alpha = .5) +
  xlab("Year") + ylab("count") +
  ggtitle("Forecasts for monthly drug overdose") +
  guides(colour=guide_legend(title="Forecast")) +
   theme_tufte()

Below are the actual forecast estimates for all drug and their prediction intervals, using seasonal naive.

	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Oct 2015	502	459.0190	544.9810	436.2662	567.7338
Nov 2015	491	448.0190	533.9810	425.2662	556.7338
Dec 2015	479	436.0190	521.9810	413.2662	544.7338
Jan 2016	483	440.0190	525.9810	417.2662	548.7338
Feb 2016	416	373.0190	458.9810	350.2662	481.7338
Mar 2016	526	483.0190	568.9810	460.2662	591.7338
Apr 2016	479	436.0190	521.9810	413.2662	544.7338
May 2016	510	467.0190	552.9810	444.2662	575.7338
Jun 2016	481	438.0190	523.9810	415.2662	546.7338
Jul 2016	534	491.0190	576.9810	468.2662	599.7338
Aug 2016	503	460.0190	545.9810	437.2662	568.7338
Sep 2016	477	434.0190	519.9810	411.2662	542.7338
Oct 2016	502	441.2157	562.7843	409.0384	594.9616
Nov 2016	491	430.2157	551.7843	398.0384	583.9616
Dec 2016	479	418.2157	539.7843	386.0384	571.9616

Residuals diagnostics

checkresiduals(snaive(drug_fc2)) 

    Ljung-Box test

data:  Residuals from Seasonal naive method
Q* = 18.011, df = 13.8, p-value = 0.1959

Model df: 0.   Total lags used: 13.8

For the overall drug overdose time-series forcast, The residuals, what is left after fitting the naive model, appear uncorrelated (p-value = 0.1959 base on the chi-square of Ljung-Box test) with a mean close to zero. However they didn’t seem to have a normal distribution, an assumption used when building the prediction intervals. So the PI may not be accurate. See below for residuals diagnostics of the opioids. All the residuals are calculated on the training sets

Any opioid

opioid_fc <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_opioid") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type)

opioid_fc2 <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_opioid") %>%
  tk_ts(., start=2010, frequency = 12, 
        end =c(2015, 9), select = -drug_type)



 opioid_fit2 <- snaive(opioid_fc2)
 
 autoplot(opioid_fc) +
 autolayer(opioid_fit2, series="Seasonal naïve", 
           PI=TRUE, col = "skyblue", size = 2, alpha = .5) +
   xlab("Year") + ylab("count") +
  ggtitle("Forecasts for monthly opioid overdose") +
  guides(colour=guide_legend(title="Forecast")) +
   theme_tufte()

For opioid, below are the actual forecast estimates and their prediction intervals, using seasonal naive.

	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Oct 2015	107	85.55209	128.4479	74.19825	139.8018
Nov 2015	108	86.55209	129.4479	75.19825	140.8018
Dec 2015	101	79.55209	122.4479	68.19825	133.8018
Jan 2016	120	98.55209	141.4479	87.19825	152.8018
Feb 2016	103	81.55209	124.4479	70.19825	135.8018
Mar 2016	133	111.55209	154.4479	100.19825	165.8018
Apr 2016	112	90.55209	133.4479	79.19825	144.8018
May 2016	112	90.55209	133.4479	79.19825	144.8018
Jun 2016	121	99.55209	142.4479	88.19825	153.8018
Jul 2016	112	90.55209	133.4479	79.19825	144.8018
Aug 2016	109	87.55209	130.4479	76.19825	141.8018
Sep 2016	111	89.55209	132.4479	78.19825	143.8018
Oct 2016	107	76.66807	137.3319	60.61132	153.3887
Nov 2016	108	77.66807	138.3319	61.61132	154.3887
Dec 2016	101	70.66807	131.3319	54.61132	147.3887

Residuals diagnostics

checkresiduals(snaive(opioid_fc2)) 

    Ljung-Box test

data:  Residuals from Seasonal naive method
Q* = 16.946, df = 13.8, p-value = 0.2471

Model df: 0.   Total lags used: 13.8

Non-Heroin Opioid

nhopi_fc <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "non_heroin_opioid") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type)

nhopi_fc2 <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "non_heroin_opioid") %>%
  tk_ts(., start=2010, frequency = 12, 
        end =c(2015, 9), select = -drug_type)



 nhopi_fit2 <- snaive(nhopi_fc2,h=15)
 
 autoplot(nhopi_fc) +
 autolayer(nhopi_fit2, series="Seasonal naïve", 
           PI=TRUE, col = "skyblue", size = 2, alpha = .5) +
  xlab("Year") + ylab("count") +
   
  ggtitle("Forecasts for monthly non-heroin opioid overdose") +
  guides(colour=guide_legend(title="Forecast")) +
   theme_tufte()

Below are the actual forecast estimates for non-heroin opioid and their prediction intervals, using seasonal naive.

	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Oct 2015	78	58.84797	97.15203	48.70950	107.2905
Nov 2015	84	64.84797	103.15203	54.70950	113.2905
Dec 2015	79	59.84797	98.15203	49.70950	108.2905
Jan 2016	93	73.84797	112.15203	63.70950	122.2905
Feb 2016	80	60.84797	99.15203	50.70950	109.2905
Mar 2016	102	82.84797	121.15203	72.70950	131.2905
Apr 2016	79	59.84797	98.15203	49.70950	108.2905
May 2016	82	62.84797	101.15203	52.70950	111.2905
Jun 2016	91	71.84797	110.15203	61.70950	120.2905
Jul 2016	88	68.84797	107.15203	58.70950	117.2905
Aug 2016	81	61.84797	100.15203	51.70950	110.2905
Sep 2016	86	66.84797	105.15203	56.70950	115.2905
Oct 2016	78	50.91494	105.08506	36.57698	119.4230
Nov 2016	84	56.91494	111.08506	42.57698	125.4230
Dec 2016	79	51.91494	106.08506	37.57698	120.4230

Residuals diagnostics

checkresiduals(snaive(nhopi_fc2)) 

    Ljung-Box test

data:  Residuals from Seasonal naive method
Q* = 21.546, df = 13.8, p-value = 0.08283

Model df: 0.   Total lags used: 13.8

Heroin

hero_fc <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "heroin") %>%
  tk_ts(., start=2010, frequency = 12, select = -drug_type)

hero_fc2 <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "heroin") %>%
  tk_ts(., start=2010, frequency = 12, 
        end =c(2015, 9), select = -drug_type)



 hero_fit2 <- snaive(hero_fc2)
 
 autoplot(hero_fc) +
 autolayer(hero_fit2, series="Seasonal naïve", 
           PI=TRUE, col = "skyblue", size = 2, alpha = .5) +
  xlab("Year") + ylab("count") +
  ggtitle("Forecasts for monthly Heroin overdose") +
  guides(colour=guide_legend(title="Forecast")) +
   theme_tufte()

Below are the actual forecast estimates for heroin and their prediction intervals, using seasonal naive.

	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Oct 2015	29	19.866798	38.1332	15.031973	42.96803
Nov 2015	24	14.866798	33.1332	10.031973	37.96803
Dec 2015	22	12.866798	31.1332	8.031973	35.96803
Jan 2016	27	17.866798	36.1332	13.031973	40.96803
Feb 2016	23	13.866798	32.1332	9.031973	36.96803
Mar 2016	31	21.866798	40.1332	17.031972	44.96803
Apr 2016	33	23.866798	42.1332	19.031972	46.96803
May 2016	30	20.866798	39.1332	16.031972	43.96803
Jun 2016	30	20.866798	39.1332	16.031972	43.96803
Jul 2016	24	14.866798	33.1332	10.031973	37.96803
Aug 2016	28	18.866798	37.1332	14.031973	41.96803
Sep 2016	25	15.866798	34.1332	11.031973	38.96803
Oct 2016	29	16.083702	41.9163	9.246226	48.75377
Nov 2016	24	11.083702	36.9163	4.246226	43.75377
Dec 2016	22	9.083702	34.9163	2.246226	41.75377
Jan 2017	27	14.083702	39.9163	7.246226	46.75377
Feb 2017	23	10.083702	35.9163	3.246226	42.75377
Mar 2017	31	18.083702	43.9163	11.246226	50.75377
Apr 2017	33	20.083702	45.9163	13.246226	52.75377
May 2017	30	17.083702	42.9163	10.246226	49.75377
Jun 2017	30	17.083702	42.9163	10.246226	49.75377
Jul 2017	24	11.083702	36.9163	4.246226	43.75377
Aug 2017	28	15.083702	40.9163	8.246226	47.75377
Sep 2017	25	12.083702	37.9163	5.246226	44.75377

Residuals diagnostics

checkresiduals(snaive(hero_fc2,h=15)) 

    Ljung-Box test

data:  Residuals from Seasonal naive method
Q* = 24.666, df = 13.8, p-value = 0.03521

Model df: 0.   Total lags used: 13.8

Concluding Remarks

• The forecasting model has been developed on the first period then tested on the second period.

• The increase in opioid hospitalization discharges in the months following the adoption of ICD-10 CM was in contrast to the downwards trend in the preceding months.

• The seasonal naive uses the last known monthly data in the training period for the forecasted data to preserve the seasonality. Its prediction intervals widen as the forecast goes further in the future.

• In fact the subsequent 2017 months show the return of the expected decrease in opioid hospitalization discharges. The 2016 bump was probably due to the transition from ICD-9-CM to ICD-10-CM.

Appendix

Comparing forecast Methods

- mean, seasonal naive and naive methods:

drug_fit1 <- meanf(drug_fc2, h = 15)
drug_fit2 <- snaive(drug_fc2, h = 15)
drug_fit3 <- rwf(drug_fc2, h = 15)

- Visualizing the different methods

autoplot(drug_fc) +
  autolayer(drug_fit1, series="Mean", PI=FALSE) +
  autolayer(drug_fit3, series="Naïve", PI=FALSE) +
  autolayer(drug_fit2, series="Seasonal naïve", PI=FALSE) +
  xlab("Month") + ylab("Count") +
  ggtitle("Forecasts for monthly drug overdose") +
  guides(colour=guide_legend(title="Forecast")) +
  theme_tufte()

- Accuracy measures comparing the three methods

drug_fc3 <- drug_hosp_2010_16_a %>% 
  filter(drug_type == "any_drug") %>%
  tk_ts(., start=c(2015, 10), frequency = 12, 
        select = -drug_type)

accuracy(drug_fit1, drug_fc3)

                        ME     RMSE      MAE        MPE     MAPE
Training set -1.729489e-14 35.42222 28.28859 -0.5152459 5.691566
Test set      2.063768e+01 29.70751 26.78261  3.7226984 4.984015
                  MASE       ACF1 Theil's U
Training set 0.9814059 0.28905030        NA
Test set     0.9291593 0.07697424  1.046967

accuracy(drug_fit2, drug_fc3)

                    ME     RMSE      MAE       MPE     MAPE     MASE
Training set -9.105263 33.53736 28.82456 -2.037625 5.824856 1.000000
Test set     41.133333 53.15951 42.20000  7.587379 7.808221 1.464029
                    ACF1 Theil's U
Training set  0.06606783        NA
Test set     -0.13692265   1.90383

accuracy(drug_fit3, drug_fc3)

                     ME     RMSE      MAE       MPE     MAPE     MASE
Training set -0.7352941 42.30526 32.88235 -0.508113  6.52251 1.140775
Test set     54.3333333 58.38436 54.33333 10.075066 10.07507 1.884967
                    ACF1 Theil's U
Training set -0.55189829        NA
Test set      0.07697424  2.026206

The function accuracy provides multiple measures of accuracy of the model fit: mean error (ME), root mean squared error (RMSE), mean absolute error (MAE), mean percentage error (MPE), mean absolute percentage error (MAPE), mean absolute scaled error (MASE) and the first-order autocorrelation coefficient (ACF1).

The seasonal naive had larger errors compared to the other methods. The errors are the deviation from the observed data in October 2015 through the end of 2016, but the seasonal naive reflects better the precedent period.