Overview
The fiphde (forecasting influenza in support of public health decision making) package provides utilities for forecasting influenza hospitalizations in the United States. fiphde includes functions for retrieving hospitalization time series data from the HHS Protect system at HealthData.gov, preparing raw data for forecasting, fitting time series and count regression models to create probabilistic forecasts for influenza hospitalizations at state and national levels, visualizing and evaluating forecasts, and formatting forecasts for submission to FluSight.
fiphde rhymes with “fifty,” as in the 50 states in the US.
Usage
The fiphde package retrieves current data from HHS and CDC APIs and fits models and forecasts using this data. This vignette uses data current to May 28, 2022 (MMWR epidemiological week 21 of 2022). Running the code here as written will produce different results depending on when you run the code, as new data is constantly being added and historical data is constantly being revised.
To get started, load the packages that are used in this vignette.
Data retrieval
Prior to fitting any forecasts we need to first retrieve data from
the HealthData.gov
COVID-19 Reported Patient Impact and Hospital Capacity by State
Timeseries API. Running get_hdgov_hosp(limitcols=TRUE)
will initiate an API call with an argument to return only a selection of
fields relevant to flu hospitalization reporting.
hosp <- get_hdgov_hosp(limitcols = TRUE)
hosp
#> # A tibble: 43,967 × 14
#> state date flu.admits flu.admits.cov flu.deaths flu.deaths.cov flu.icu
#> <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 AL 2019-12-31 NA 0 NA 0 NA
#> 2 HI 2019-12-31 NA 0 NA 0 NA
#> 3 IN 2019-12-31 NA 0 NA 0 NA
#> 4 LA 2019-12-31 NA 0 NA 0 NA
#> 5 MN 2019-12-31 NA 0 NA 0 NA
#> 6 MT 2019-12-31 NA 0 NA 0 NA
#> 7 NC 2019-12-31 NA 0 NA 0 NA
#> 8 PR 2019-12-31 NA 0 NA 0 NA
#> 9 TX 2019-12-31 NA 0 NA 0 NA
#> 10 AL 2020-01-01 NA 0 NA 0 NA
#> # ℹ 43,957 more rows
#> # ℹ 7 more variables: flu.icu.cov <dbl>, flu.tot <dbl>, flu.tot.cov <dbl>,
#> # cov.admits <dbl>, cov.admits.cov <dbl>, cov.deaths <dbl>,
#> # cov.deaths.cov <dbl>Time series forecasting
We will first fit a time series model, creating an ensemble model from ARIMA and exponential smoothing models. Time series modeling is based on the tidyverts (https://tidyverts.org/) collection of packages for tidy time series forecasting in R.
Data preparation
We need to initially prepare the data for a time series forecast. The
prep_hdgov_hosp function call below will limit to states
only (no territories), remove any data from an incomplete
epidemiological week (Sunday-Saturday), remove locations with little to
no reported hospitalizations over the last month, and further exclude
Washington DC from downstream analysis. The function aggregates total
number of cases at each epidemiological week at each location. This
function also adds in location FIPS codes, and joins in historical
influenza-like illness (ILI) and hospitalization mean values and ranks
by week. Historical indicators for ILI and hospitalizations are
summarized from CDC
ILINet and CDC
FluSurv-Net respectively.
# Prep data
prepped_hosp <-
hosp %>%
prep_hdgov_hosp(statesonly=TRUE, min_per_week = 0, remove_incomplete = TRUE) %>%
dplyr::filter(abbreviation != "DC")
prepped_hosp
#> # A tibble: 4,284 × 13
#> abbreviation location week_start monday week_end epiyear epiweek
#> <chr> <chr> <date> <date> <date> <dbl> <dbl>
#> 1 US US 2020-10-18 2020-10-19 2020-10-24 2020 43
#> 2 US US 2020-10-25 2020-10-26 2020-10-31 2020 44
#> 3 US US 2020-11-01 2020-11-02 2020-11-07 2020 45
#> 4 US US 2020-11-08 2020-11-09 2020-11-14 2020 46
#> 5 US US 2020-11-15 2020-11-16 2020-11-21 2020 47
#> 6 US US 2020-11-22 2020-11-23 2020-11-28 2020 48
#> 7 US US 2020-11-29 2020-11-30 2020-12-05 2020 49
#> 8 US US 2020-12-06 2020-12-07 2020-12-12 2020 50
#> 9 US US 2020-12-13 2020-12-14 2020-12-19 2020 51
#> 10 US US 2020-12-20 2020-12-21 2020-12-26 2020 52
#> # ℹ 4,274 more rows
#> # ℹ 6 more variables: flu.admits <dbl>, flu.admits.cov <dbl>, ili_mean <dbl>,
#> # ili_rank <int>, hosp_mean <dbl>, hosp_rank <int>Now let’s explore the data.
prepped_hosp %>%
filter(abbreviation %in% c("US", "CA", "TX", "NY")) %>%
ggplot(aes(week_end, flu.admits)) + geom_line() +
facet_wrap(~abbreviation, scale="free_y")
What states had the highest admissions over the 2021-2022 flu season?
prepped_hosp %>%
filter(abbreviation!="US") %>%
filter(week_start>="2021-07-01" & week_end<"2022-06-30") %>%
group_by(abbreviation) %>%
summarize(total.flu.admits=sum(flu.admits)) %>%
arrange(desc(total.flu.admits)) %>%
head(10) %>%
knitr::kable(caption="Top 10 states with highest flu hospitalizations in 2021-2022.")| abbreviation | total.flu.admits |
|---|---|
| TX | 7164 |
| FL | 4550 |
| PA | 3286 |
| CA | 3233 |
| NY | 3073 |
| AZ | 2819 |
| MI | 2511 |
| OK | 2358 |
| OH | 2169 |
| MA | 2141 |
Next let’s turn this into a tsibble. tsibble objects are
tibbles with an index variable describing the inherent ordering
from past to present, and a key that defines observational
units over time. The make_tsibble function provides a
convenience wrapper around tsibble::as_tsibble using the
epidemiological week’s Monday as the weekly index and the location as
the key. Note that the specification of arguments for epidemiological
week / year and location key are passed as “bare” (unquoted) names of
columns storing this information in the original tibble.
prepped_hosp_tsibble <- make_tsibble(prepped_hosp,
epiyear=epiyear,
epiweek=epiweek,
key=location)
prepped_hosp_tsibble
#> # A tsibble: 4,284 x 14 [1W]
#> # Key: location [51]
#> abbreviation location week_start monday yweek week_end epiyear
#> <chr> <chr> <date> <date> <week> <date> <dbl>
#> 1 AL 01 2020-10-18 2020-10-19 2020 W43 2020-10-24 2020
#> 2 AL 01 2020-10-25 2020-10-26 2020 W44 2020-10-31 2020
#> 3 AL 01 2020-11-01 2020-11-02 2020 W45 2020-11-07 2020
#> 4 AL 01 2020-11-08 2020-11-09 2020 W46 2020-11-14 2020
#> 5 AL 01 2020-11-15 2020-11-16 2020 W47 2020-11-21 2020
#> 6 AL 01 2020-11-22 2020-11-23 2020 W48 2020-11-28 2020
#> 7 AL 01 2020-11-29 2020-11-30 2020 W49 2020-12-05 2020
#> 8 AL 01 2020-12-06 2020-12-07 2020 W50 2020-12-12 2020
#> 9 AL 01 2020-12-13 2020-12-14 2020 W51 2020-12-19 2020
#> 10 AL 01 2020-12-20 2020-12-21 2020 W52 2020-12-26 2020
#> # ℹ 4,274 more rows
#> # ℹ 7 more variables: epiweek <dbl>, flu.admits <dbl>, flu.admits.cov <dbl>,
#> # ili_mean <dbl>, ili_rank <int>, hosp_mean <dbl>, hosp_rank <int>Fit a model and forecast
Next, let’s fit a time series model and create forecasts using the
ts_fit_forecast function. This function takes a tsibble
created as above, a forecast horizon in weeks, the name of the outcome
variable to forecast, and optional covariates to use in the ARIMA
model.
Here we fit a non-seasonal ARIMA model with the autoregressive term
(p) restricted to 1:2, order of integration for differencing (d)
restricted to 0:2, and the moving average (q) restricted to 0 (see
?fable::ARIMA for more information). The model also fits a
non-seasonal exponential smoothing model (see ?fable::ETS
for details). In this example we do not fit an autoregressive neural
network model, but we could change nnetar=NULL to
nnetar="AR(P=1)" to do so (see ?fable::nnetar
for details). We can also trim the data used in modeling, and here we
restrict to only include hospitalization data reported after January 1,
2021. By setting ensemble=TRUE we create an ensemble model
that averages the ARIMA and exponential smoothing models.
hosp_fitfor <- ts_fit_forecast(prepped_hosp_tsibble,
horizon=4L,
outcome="flu.admits",
trim_date = "2021-01-01",
covariates=TRUE,
models=list(arima='PDQ(0, 0, 0) + pdq(1:2, 0:2, 0)',
ets='season(method="N")',
nnetar=NULL),
ensemble=TRUE)The function will output messages describing the ARIMA and ETS model
formulas passed to fable::ARIMA and
fable::ETS.
Trimming to 2021-01-01
ARIMA formula: flu.admits ~ PDQ(0, 0, 0) + pdq(1:2, 0:2, 0) + hosp_rank + ili_rank
ETS formula: flu.admits ~ season(method = "N")After fitting the models, the function then forecasts the outcome to
the specified number of weeks. Let’s take a look at the object returned
from this model fit + forecast. We see $tsfit, which gives
us the ARIMA, ETS, and ensemble model fits as list columns, one row per
location; $tsfor gives us the forecast for the next four
weeks, one row per model per location (key in the tsibble);
$formulas gives us the model formulas passed to fable
modeling functions; and if any models failed to converge we would see
those locations in $nullmodels.
hosp_fitfor
#> $tsfit
#> # A mable: 51 x 4
#> # Key: location [51]
#> location arima ets ensemble
#> <chr> <model> <model> <model>
#> 1 01 <LM w/ ARIMA(2,0,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 2 02 <LM w/ ARIMA(1,0,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 3 04 <LM w/ ARIMA(2,1,0) errors> <ETS(M,N,N)> <COMBINATION>
#> 4 05 <LM w/ ARIMA(2,1,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 5 06 <LM w/ ARIMA(2,2,0) errors> <ETS(M,A,N)> <COMBINATION>
#> 6 08 <LM w/ ARIMA(1,1,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 7 09 <LM w/ ARIMA(1,1,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 8 10 <LM w/ ARIMA(2,1,0) errors> <ETS(A,N,N)> <COMBINATION>
#> 9 12 <LM w/ ARIMA(2,2,0) errors> <ETS(M,N,N)> <COMBINATION>
#> 10 13 <LM w/ ARIMA(2,1,0) errors> <ETS(A,N,N)> <COMBINATION>
#> # ℹ 41 more rows
#>
#> $tsfor
#> # A fable: 612 x 10 [1W]
#> # Key: location, .model [153]
#> location .model yweek flu.admits .mean epiweek ili_mean ili_rank hosp_mean
#> <chr> <chr> <week> <dist> <dbl> <dbl> <dbl> <int> <dbl>
#> 1 01 arima 2022 W22 N(16, 56) 16.2 22 1.23 17 0
#> 2 01 arima 2022 W23 N(20, 56) 20.4 23 1.09 14 0
#> 3 01 arima 2022 W24 N(19, 65) 19.2 24 1.02 12 0
#> 4 01 arima 2022 W25 N(21, 65) 21.0 25 0.984 10 0
#> 5 01 ets 2022 W22 N(19, 78) 18.9 22 1.23 17 0
#> 6 01 ets 2022 W23 N(19, 88) 18.9 23 1.09 14 0
#> 7 01 ets 2022 W24 N(19, 99) 18.9 24 1.02 12 0
#> 8 01 ets 2022 W25 N(19, 110) 18.9 25 0.984 10 0
#> 9 01 ensem… 2022 W22 N(18, 63) 17.5 22 1.23 17 0
#> 10 01 ensem… 2022 W23 N(20, 67) 19.7 23 1.09 14 0
#> # ℹ 602 more rows
#> # ℹ 1 more variable: hosp_rank <int>
#>
#> $formulas
#> $formulas$arima
#> flu.admits ~ PDQ(0, 0, 0) + pdq(1:2, 0:2, 0) + hosp_rank + ili_rank
#> <environment: 0x55b086e3bce8>
#>
#> $formulas$ets
#> flu.admits ~ season(method = "N")
#> <environment: 0x55b086e3bce8>
#>
#>
#> $nullmodels
#> # A tibble: 0 × 2
#> # ℹ 2 variables: location <chr>, model <chr>Next, we can format the forecasts for submission
to FluSight using the format_for_submission
function.
formatted_list <- format_for_submission(hosp_fitfor$tsfor, method = "ts", format = "legacy")The list contains separate submission-ready tibbles, one element for each type of model fitted.
formatted_list
#> $arima
#> # A tibble: 4,896 × 7
#> forecast_date target target_end_date location type quantile value
#> <date> <chr> <date> <chr> <chr> <chr> <chr>
#> 1 2022-12-20 1 wk ahead inc f… 2022-06-04 01 point NA 17
#> 2 2022-12-20 2 wk ahead inc f… 2022-06-11 01 point NA 21
#> 3 2022-12-20 3 wk ahead inc f… 2022-06-18 01 point NA 20
#> 4 2022-12-20 4 wk ahead inc f… 2022-06-25 01 point NA 22
#> 5 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.010 0
#> 6 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.010 3
#> 7 2022-12-20 3 wk ahead inc f… 2022-06-18 01 quan… 0.010 1
#> 8 2022-12-20 4 wk ahead inc f… 2022-06-25 01 quan… 0.010 3
#> 9 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.025 2
#> 10 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.025 6
#> # ℹ 4,886 more rows
#>
#> $ensemble
#> # A tibble: 4,896 × 7
#> forecast_date target target_end_date location type quantile value
#> <date> <chr> <date> <chr> <chr> <chr> <chr>
#> 1 2022-12-20 1 wk ahead inc f… 2022-06-04 01 point NA 18
#> 2 2022-12-20 2 wk ahead inc f… 2022-06-11 01 point NA 20
#> 3 2022-12-20 3 wk ahead inc f… 2022-06-18 01 point NA 20
#> 4 2022-12-20 4 wk ahead inc f… 2022-06-25 01 point NA 20
#> 5 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.010 0
#> 6 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.010 1
#> 7 2022-12-20 3 wk ahead inc f… 2022-06-18 01 quan… 0.010 0
#> 8 2022-12-20 4 wk ahead inc f… 2022-06-25 01 quan… 0.010 0
#> 9 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.025 3
#> 10 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.025 4
#> # ℹ 4,886 more rows
#>
#> $ets
#> # A tibble: 4,896 × 7
#> forecast_date target target_end_date location type quantile value
#> <date> <chr> <date> <chr> <chr> <chr> <chr>
#> 1 2022-12-20 1 wk ahead inc f… 2022-06-04 01 point NA 19
#> 2 2022-12-20 2 wk ahead inc f… 2022-06-11 01 point NA 19
#> 3 2022-12-20 3 wk ahead inc f… 2022-06-18 01 point NA 19
#> 4 2022-12-20 4 wk ahead inc f… 2022-06-25 01 point NA 19
#> 5 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.010 0
#> 6 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.010 0
#> 7 2022-12-20 3 wk ahead inc f… 2022-06-18 01 quan… 0.010 0
#> 8 2022-12-20 4 wk ahead inc f… 2022-06-25 01 quan… 0.010 0
#> 9 2022-12-20 1 wk ahead inc f… 2022-06-04 01 quan… 0.025 2
#> 10 2022-12-20 2 wk ahead inc f… 2022-06-11 01 quan… 0.025 1
#> # ℹ 4,886 more rowsWe can check to see if the submission is valid.
Note that this will fail if the expected date of the target dates and
current dates don’t line up. You’ll see a message noting “The
submission target end dates do not line up with expected Saturdays by
horizon. Note if submission forecast date is not Sunday or Monday, then
forecasts are assumed to to start the following week.” If
validation succeeds, the $valid element of the returned
list will be TRUE, and FALSE if any validation
checks fail.
validate_forecast(formatted_list$ensemble)Let’s plot the forecast with the observed data using the
plot_forecast function. Here we plot forecasts with the 50%
prediction interval for US, New York (FIPS 36) and Florida (FIPS
12).
plot_forecast(prepped_hosp, formatted_list$ensemble, loc="US", pi = .5)
plot_forecast(prepped_hosp, formatted_list$ensemble, loc="36", pi = .5)
plot_forecast(prepped_hosp, formatted_list$ensemble, loc="12", pi = .5)
Finally, we can pull out the ARIMA model parameters used for each location to save for posterity or retrospective analysis.
hosp_fitfor$tsfit$arima %>%
map("fit") %>%
map_df("spec") %>%
mutate(location = hosp_fitfor$tsfit$location, .before = "p")
#> # A tibble: 51 × 9
#> location p d q P D Q constant period
#> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <lgl> <dbl>
#> 1 01 2 0 0 0 0 0 TRUE 52
#> 2 02 1 0 0 0 0 0 TRUE 52
#> 3 04 2 1 0 0 0 0 FALSE 52
#> 4 05 2 1 0 0 0 0 FALSE 52
#> 5 06 2 2 0 0 0 0 FALSE 52
#> 6 08 1 1 0 0 0 0 FALSE 52
#> 7 09 1 1 0 0 0 0 FALSE 52
#> 8 10 2 1 0 0 0 0 FALSE 52
#> 9 12 2 2 0 0 0 0 FALSE 52
#> 10 13 2 1 0 0 0 0 FALSE 52
#> # ℹ 41 more rowsCount regression forecasting
The time series methods above assume that the outcome has a continuous distribution. When forecasting counts and especially small counts (e.g., “how many influenza hospitalizations will occur next week in Hawaii?”) alternative methods may have more desirable properties. fiphde provides functionality to leverage count regression models for forecasting influenza hospitalizations. Here we will demonstrate how to implement a count regression modeling approach. The example will step through the process of forecasting hospitalizations in Hawaii by using influenza-like illness (ILI) and indicators of historical hospitalization and ILI severity for the given epidemiological weeks. The usage illustrates an automated tuning procedure that finds the “best” models from possible covariates and model families (e.g., Poisson, Quasipoisson, Negative binomial, etc.).
ILI retrieval and prep
Given that ILI will be a covariate in the count regression model, we
must first retrieve ILI data using the get_cdc_ili()
function, which is a wrapper around ilinet(). Here, we’re
only using recent ILI data (i.e., since 2019) and we further filter the
signal to only include Hawaii. ILI is subject to revision and may be
less reliable when initially reported, so we replace the current week
and one week previous with Nowcast data. Finally, we fit
a time series model to forecast ILI for the next four weeks.
ilidat <-
get_cdc_ili(region=c("state"), years=2019:2022) %>%
filter(region == "Hawaii") %>%
replace_ili_nowcast(., weeks_to_replace=1)
ilifor <- forecast_ili(ilidat, horizon=4L, trim_date="2020-03-01")In the next step we are going to log transform weighted ILI and flu
admissions, so in this step we need to remove all the zeros. The fiphde
package provides the mnz_replace function which replaces
zeros in a variable with the smallest non-zero value of that
variable.
ilidat <- ilidat %>% mutate(weighted_ili=mnz_replace(weighted_ili))
ilifor$ilidat <- ilifor$ilidat %>% mutate(ili=mnz_replace(ili))
ilifor$ili_future <- ilifor$ili_future %>% mutate(ili=mnz_replace(ili))
ilifor$ili_bound <- ilifor$ili_bound %>% mutate(ili=mnz_replace(ili))Finally, we create our modeling dataset by combining the flu admission data prepared above with the forecasted ILI data in future weeks together with the historical severity data by epidemiological week.
dat_hi <-
prepped_hosp %>%
filter(abbreviation=="HI") %>%
dplyr::mutate(date = MMWRweek::MMWRweek2Date(epiyear, epiweek)) %>%
left_join(ilifor$ilidat, by = c("epiyear", "location", "epiweek")) %>%
mutate(ili = log(ili))
dat_hi
#> # A tibble: 84 × 15
#> abbreviation location week_start monday week_end epiyear epiweek
#> <chr> <chr> <date> <date> <date> <dbl> <dbl>
#> 1 HI 15 2020-10-18 2020-10-19 2020-10-24 2020 43
#> 2 HI 15 2020-10-25 2020-10-26 2020-10-31 2020 44
#> 3 HI 15 2020-11-01 2020-11-02 2020-11-07 2020 45
#> 4 HI 15 2020-11-08 2020-11-09 2020-11-14 2020 46
#> 5 HI 15 2020-11-15 2020-11-16 2020-11-21 2020 47
#> 6 HI 15 2020-11-22 2020-11-23 2020-11-28 2020 48
#> 7 HI 15 2020-11-29 2020-11-30 2020-12-05 2020 49
#> 8 HI 15 2020-12-06 2020-12-07 2020-12-12 2020 50
#> 9 HI 15 2020-12-13 2020-12-14 2020-12-19 2020 51
#> 10 HI 15 2020-12-20 2020-12-21 2020-12-26 2020 52
#> # ℹ 74 more rows
#> # ℹ 8 more variables: flu.admits <dbl>, flu.admits.cov <dbl>, ili_mean <dbl>,
#> # ili_rank <int>, hosp_mean <dbl>, hosp_rank <int>, date <date>, ili <dbl>Fit a model and forecast
First we create a list of different count regression models to fit. We can inspect the formulations of each model. Note that these models are specified using the trending package and internally the trendeval package is used to determine the best fitting approach.
models <-
list(
poisson = trending::glm_model(flu.admits ~ ili + hosp_rank + ili_rank, family = "poisson"),
quasipoisson = trending::glm_model(flu.admits ~ ili + hosp_rank + ili_rank, family = "quasipoisson"),
negbin = trending::glm_nb_model(flu.admits ~ ili + hosp_rank + ili_rank)
)
models$poisson
#> Untrained trending model:
#> glm(formula = flu.admits ~ ili + hosp_rank + ili_rank, family = "poisson")
models$quasipoisson
#> Untrained trending model:
#> glm(formula = flu.admits ~ ili + hosp_rank + ili_rank, family = "quasipoisson")
models$negbin
#> Untrained trending model:
#> glm.nb(formula = flu.admits ~ ili + hosp_rank + ili_rank)We must next project covariates four weeks into the future, which we can accomplish by pulling from the forecasted ILI and historical severity data.
new_cov <-
ilifor$ili_future %>%
left_join(fiphde:::historical_severity, by="epiweek") %>%
select(-epiweek,-epiyear) %>%
mutate(ili = log(ili))
new_cov
#> # A tibble: 4 × 6
#> location ili ili_mean ili_rank hosp_mean hosp_rank
#> <chr> <dbl> <dbl> <int> <dbl> <int>
#> 1 15 0.655 1.02 12 0 0
#> 2 15 0.634 0.984 10 0 0
#> 3 15 0.642 0.946 9 0 0
#> 4 15 0.644 0.921 8 0 0Next we use fiphde’s glm_wrap function which attempts to
find the best-fitting model out of all the models supplied, and
specifies the prediction interval quantiles used in FluSight.
res <- glm_wrap(dat_hi,
new_covariates = new_cov,
.models = models,
alpha = c(0.01, 0.025, seq(0.05, 0.5, by = 0.05)) * 2)
res$forecasts$location <- "15"We can take a quick look at the forecast output as well as different aspects of the final fitted model.
head(res$forecasts)
#> epiweek epiyear quantile value location
#> 1 22 2022 0.010 0 15
#> 2 22 2022 0.025 0 15
#> 3 22 2022 0.050 0 15
#> 4 22 2022 0.100 0 15
#> 5 22 2022 0.150 1 15
#> 6 22 2022 0.200 1 15
res$model
#> # A tibble: 2 × 4
#> model_class fit location data$fit_data
#> <chr> <trndng__> <chr> <list>
#> 1 glm <glm> 15 <tibble [84 × 15]>
#> 2 glm <fn> 15 <tibble [84 × 15]>
res$model$fit
#> $fitted_model
#>
#> Call: glm(formula = flu.admits ~ ili + hosp_rank + ili_rank, family = "poisson",
#> data = data)
#>
#> Coefficients:
#> (Intercept) ili hosp_rank ili_rank
#> -0.110528 1.463418 -0.051589 0.003076
#>
#> Degrees of Freedom: 83 Total (i.e. Null); 80 Residual
#> Null Deviance: 229.6
#> Residual Deviance: 175.6 AIC: 234
#>
#> $predict
#> function (newdata, alpha = 0.05, add_pi = TRUE, simulate_pi = TRUE,
#> uncertain = TRUE)
#> {
#> if (missing(newdata)) {
#> newdata <- data[all.vars(formula(model))]
#> }
#> result <- add_confidence_interval(model, newdata, alpha)
#> if (add_pi) {
#> if (simulate_pi) {
#> result <- add_prediction_interval(model, result,
#> alpha, simulate_pi, uncertain)
#> }
#> else {
#> result <- add_prediction_interval(model, result,
#> alpha, simulate_pi, uncertain)
#> }
#> }
#> result
#> }
#> <bytecode: 0x55b093ac6130>
#> <environment: 0x55b089669cb8>
#>
#> attr(,"class")
#> [1] "trending_model_fit" "list"
res$model$fit$fitted_model$family
#>
#> Family: poisson
#> Link function: log
res$model$fit$fitted_model$coefficients
#> (Intercept) ili hosp_rank ili_rank
#> -0.110527528 1.463417923 -0.051589427 0.003075992Next we prepare the data in the quantile format used by FluSight.
hi_glm_prepped <- format_for_submission(res$forecasts, method="CREG", format = "legacy")
hi_glm_prepped
#> $CREG
#> # A tibble: 96 × 7
#> forecast_date target target_end_date location type quantile value
#> <date> <chr> <date> <chr> <chr> <chr> <chr>
#> 1 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.010 0
#> 2 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.025 0
#> 3 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.050 0
#> 4 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.100 0
#> 5 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.150 1
#> 6 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.200 1
#> 7 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.250 1
#> 8 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.300 1
#> 9 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.350 2
#> 10 2024-08-01 1 wk ahead inc f… 2022-06-04 15 quan… 0.400 2
#> # ℹ 86 more rowsFinally, we can visualize the forecasted point estimates with the 90% prediction interval alongside the observed data.
plot_forecast(dat_hi, hi_glm_prepped$CREG, location="15", pi=.9)