Prepare observed data

To motivate this example we will use incident flu hospitalization reported via HHS Protect. Note that the data here has been aggregated from daily to weekly resolution, and is provided in this format as internal rplanes package data. We select only the columns in which we’re interested (“date”, “location”, and “flu.admits”) and make sure that the date field is formatted as a date:

hosp_all <-
  read.csv(system.file("extdata/observed/hdgov_hosp_weekly.csv", package = "rplanes")) %>%
  select(date, location, flu.admits) %>%
  mutate(date = as.Date(date))

head(hosp_all)

It is important to emphasize that the observed data must at minimum include columns for location (geographic unit such as FIPS code), date (date of reported value; must be date class), and a reported value for each location/date combination.

Convert observed data into a signal with to_signal()

The downstream plausibility analysis requires that observed and/or forecasted data be converted to an S3 “signal” class. The to_signal() function constructs this object and applies an additional “observed” or “forecast” class depending on the type of signal specified.

Here we will convert the observed data from above to a signal object:

observed_signal <- to_signal(input = hosp_all, outcome = "flu.admits", type = "observed", resolution = "weeks", horizon = NULL)

Create a seed from the observed data via plane_seed()

The plane_seed() function provides a means to internally calculate and store baseline characteristics for the observed data in what the package calls a “seed”. These characteristics are summarized independently for each location, and are then used in downstream scoring functions.

In this example, we input the prepped observed signal object and make use of the “cut_date” argument. By restricting to dates on or before the specified cut date, we ensure that the baseline characteristics are not computed based on data that overlaps with the forecast. In practice, it may be possible to use the default (NULL) cut date when analyzing forecasts, so long as the reported data is available up to (and not exceeding) the first horizon of the forecast.

prepped_seed <- plane_seed(observed_signal, cut_date = "2022-10-29")

Prepare forecast data

The forecast data to be assessed must include a point estimate and some representation of uncertainty (i.e., prediction interval) for each location and horizon. For analysis in rplanes, the forecast data must be prepared to include the following:

• location: Geographic unit such as FIPS code
• date: Date corresponding the forecast horizon
• horizon: Forecast horizon
• lower: Lower limit of a given prediction interval for the forecast
• point: Point estimate for the forecast
• upper: Upper limit of a given prediction interval for the forecast

For convenience, the rplanes package includes read_forecast() to convert forecasts from a quantile format used by many forecast hubs^1,2 into the format needed for plausibility analysis. The read_forecast() function reads csv input and returns a tibble with summarized forecast data for each location and horizon in the original file.

For this example, we use a forecast of incident weekly flu hospitalizations. This forecast is originally in the quantile format referenced above and is included in rplanes as example data:

forecast_fp <- system.file("extdata/forecast/2022-10-31-SigSci-TSENS.csv", package = "rplanes")

read.csv(forecast_fp) %>% 
  head(.)

We can use read_forecast() to load the data and convert to the required format:

prepped_forecast <- read_forecast(forecast_fp)

head(prepped_forecast)

Convert forecast data into a signal with to_signal()

With the forecast formatted appropriately, we can convert it to a signal object:

forecast_signal <- 
  prepped_forecast  %>%
  to_signal(., outcome = "flu.admits", type = "forecast", resolution = "weeks", horizon = 4)

Run plane_score() and visualize results

The planes_score() function wraps PLANES scoring for plausibility components across all locations in a single step. The function can accept arguments to specific component functions and allows the user to dictate which components to use (with a default to all). For more about these arguments see ?plane_score.

We pass the input forecast signal to the function along with the prepped seed to run the scoring independently at all locations:

scores <- plane_score(input = forecast_signal, seed = prepped_seed)

Note that if we had other forecasts to evaluate for the same time period, we could prepare each forecast signal appropriately and use the same seed (i.e., we would not have to recreate the seed for each forecast).

The output of plane_score() is a list of scoring results for all locations. We can convert the scores_summary element from that list to a tibble. By default, the score uses equal weights for all components. However, users can modify this behavior using the “weights” argument as needed. The res object below shows the location, the number of flags raised (“n_flags”), the weighted flags (“n_flags_weighted”), the number of components tested (“n_components”), the overall sum of weights for components tested (“weights_denominator”), and the “score” (which is “n_flags_weighted” divided by “weights_denominator”), the names of components tested, and the names of the components that were flagged (if any):

res <-
  scores$scores_summary %>%
  map_df(., as_tibble)
  
head(res)

res %>%
  count(n_flags) %>%
  mutate(n_flags = as.character(n_flags)) %>%
  ggplot(aes(n_flags,n)) +
  geom_col() +
  labs(x = "Number of flags raised", y = "Count")

For the evaluated forecast, most of the locations had one or fewer of the flags raised. However, several locations had two or three of the seven possible components flagged for implausibility.

Prepare observed data

For this example, we will also use the reported HHS Protect flu hospitalization data. However, we will truncate this data to stop at the end of 2022. In this scenario, we will be interested in assessing the most recent reported week (i.e., the week of 2022-12-31) for plausibility.

Note that in practice, it may be of interest to assess multiple time points for recently observed data.

To begin, we will read in the data and filter it to dates prior to 2023:

hosp_pre23 <-
  read.csv(system.file("extdata/observed/hdgov_hosp_weekly.csv", package = "rplanes")) %>%
  select(date, location, flu.admits) %>%
  mutate(date = as.Date(date)) %>%
  filter(date < as.Date("2023-01-01"))

head(hosp_pre23)

Convert observed data into a signal with to_signal()

With the data loaded, we can convert the hospitalization tibble to a signal object:

observed_signal <- to_signal(input = hosp_pre23, outcome = "flu.admits", type = "observed", resolution = "weeks", horizon = NULL)

Create a seed from the observed data via plane_seed()

We can use the observed the signal as the input to plane_score(), but also as the input to plane_seed(). However, to do so we must use an argument that restricts the seeding to only dates prior to the cut date (inclusive).

In this example, because we are assessing the observed data for the week of 2022-12-31 we would want to seed our plausibility analysis with characteristics up through the most recent reported date (i.e., 2022-12-24):

prepped_seed <- plane_seed(observed_signal, cut_date = "2022-12-24")

Run plane_score() and visualize results

As with the forecast assessments, the PLANES scoring is wrapped for convenience in the plane_score() function. Since we are assessing an observed signal in this case, we can explicitly tell the scoring to use certain components. It is important to emphasize that not all components will work for both forecasts and observed signals. Currently, plane_cover(), plane_taper(), plane_trend(), and plane_shape() (all of which are run inside plane_score()) work only for forecasts and not for observed data:

scores <- plane_score(observed_signal, seed = prepped_seed, components = c("repeat","diff"))

res <-
  scores$scores_summary %>%
  map_df(., as_tibble)
  
head(res)

res %>%
  count(n_flags) %>%
  mutate(n_flags = as.character(n_flags)) %>%
  ggplot(aes(n_flags,n)) +
  geom_col() +
  labs(x = "Number of flags raised", y = "Count")

These results indicate that there is just one location for which any flag has been raised. Again, in this case only two components are used for scoring so the results should be interpreted accordingly. Furthermore, we have to consider that for reported signals the date at which the observed data is retrieved is paramount. It will likely be most meaningful to analyze plausibility of observed data in real time (or near real time) as opposed to retrospectively, since backfill and reporting issues may be corrected prior to assessment.