healthyverse_tsa

Time Series Analysis and Nested Modeling of the Healthyverse Packages

Steven P. Sanderson II, MPH - Date: 02 December, 2023

This analysis follows a Nested Modeltime Workflow.

Get Data

glimpse(downloads_tbl)
## Rows: 83,636
## Columns: 11
## $ date      <date> 2020-11-23, 2020-11-23, 2020-11-23, 2020-11-23, 2020-11-23,…
## $ time      <Period> 15H 36M 55S, 11H 26M 39S, 23H 34M 44S, 18H 39M 32S, 9H 0M…
## $ date_time <dttm> 2020-11-23 15:36:55, 2020-11-23 11:26:39, 2020-11-23 23:34:…
## $ size      <int> 4858294, 4858294, 4858301, 4858295, 361, 4863722, 4864794, 4…
## $ r_version <chr> NA, "4.0.3", "3.5.3", "3.5.2", NA, NA, NA, NA, NA, NA, NA, N…
## $ r_arch    <chr> NA, "x86_64", "x86_64", "x86_64", NA, NA, NA, NA, NA, NA, NA…
## $ r_os      <chr> NA, "mingw32", "mingw32", "linux-gnu", NA, NA, NA, NA, NA, N…
## $ package   <chr> "healthyR.data", "healthyR.data", "healthyR.data", "healthyR…
## $ version   <chr> "1.0.0", "1.0.0", "1.0.0", "1.0.0", "1.0.0", "1.0.0", "1.0.0…
## $ country   <chr> "US", "US", "US", "GB", "US", "US", "DE", "HK", "JP", "US", …
## $ ip_id     <int> 2069, 2804, 78827, 27595, 90474, 90474, 42435, 74, 7655, 638…

The last day in the data set is 2023-11-30 23:48:34, the file was birthed on: 2022-07-02 23:58:17, and at report knit time is -1.237984^{4} hours old. Happy analyzing!

Now that we have our data lets take a look at it using the skimr package.

skim(downloads_tbl)
   
Name downloads_tbl
Number of rows 83636
Number of columns 11
_______________________  
Column type frequency:  
character 6
Date 1
numeric 2
POSIXct 1
Timespan 1
________________________  
Group variables None

Data summary

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
r_version 57374 0.31 5 5 0 39 0
r_arch 57374 0.31 3 7 0 5 0
r_os 57374 0.31 7 15 0 17 0
package 0 1.00 7 13 0 7 0
version 0 1.00 5 6 0 49 0
country 6951 0.92 2 2 0 147 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
date 0 1 2020-11-23 2023-11-30 2022-07-15 1103

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
size 0 1 1260453.36 1617169.95 357 16847.25 322874 2429696 5677952 ▇▁▂▁▁
ip_id 0 1 10271.52 18218.43 1 184.00 2899 11412 143633 ▇▁▁▁▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
date_time 0 1 2020-11-23 09:00:41 2023-11-30 23:48:34 2022-07-15 19:07:28 50540

Variable type: Timespan

skim_variable n_missing complete_rate min max median n_unique
time 0 1 0 59 51 60

We can see that the following columns are missing a lot of data and for us are most likely not useful anyways, so we will drop them c(r_version, r_arch, r_os)

Plots

Now lets take a look at a time-series plot of the total daily downloads by package. We will use a log scale and place a vertical line at each version release for each package.

Now lets take a look at some time series decomposition graphs.

Feature Engineering

Now that we have our basic data and a shot of what it looks like, let’s add some features to our data which can be very helpful in modeling. Lets start by making a tibble that is aggregated by the day and package, as we are going to be interested in forecasting the next 4 weeks or 28 days for each package. First lets get our base data.

Now we are going to do some basic pre-processing.

data_padded_tbl <- base_data %>%
  pad_by_time(
    .date_var  = date,
    .pad_value = 0
  )

# Get log interval and standardization parameters
log_params  <- liv(data_padded_tbl$value, limit_lower = 0, offset = 1, silent = TRUE)
limit_lower <- log_params$limit_lower
limit_upper <- log_params$limit_upper
offset      <- log_params$offset

data_liv_tbl <- data_padded_tbl %>%
  # Get log interval transform
  mutate(value_trans = liv(value, limit_lower = 0, offset = 1, silent = TRUE)$log_scaled)

# Get Standardization Params
std_params <- standard_vec(data_liv_tbl$value_trans, silent = TRUE)
std_mean   <- std_params$mean
std_sd     <- std_params$sd

data_transformed_tbl <- data_liv_tbl %>%
  # get standardization
  mutate(value_trans = standard_vec(value_trans, silent = TRUE)$standard_scaled) %>%
  select(-value)

Since this is panel data we can follow one of two different modeling strategies. We can search for a global model in the panel data or we can use nested forecasting finding the best model for each of the time series. Since we only have 5 panels, we will use nested forecasting.

To do this we will use the nest_timeseries and split_nested_timeseries functions to create a nested tibble.

horizon <- 4*7

nested_data_tbl <- data_transformed_tbl %>%
    
    # 1. Extending: We'll predict n days into the future.
    extend_timeseries(
        .id_var        = package,
        .date_var      = date,
        .length_future = horizon
    ) %>%
    
    # 2. Nesting: We'll group by id, and create a future dataset
    #    that forecasts n days of extended data and
    #    an actual dataset that contains n*2 days
    nest_timeseries(
        .id_var        = package,
        .length_future = horizon
        #.length_actual = horizon*2
    ) %>%
    
   # 3. Splitting: We'll take the actual data and create splits
   #    for accuracy and confidence interval estimation of n das (test)
   #    and the rest is training data
    split_nested_timeseries(
        .length_test = horizon
    )

nested_data_tbl
## # A tibble: 7 × 4
##   package       .actual_data       .future_data      .splits         
##   <fct>         <list>             <list>            <list>          
## 1 TidyDensity   <tibble [550 × 2]> <tibble [28 × 2]> <split [522|28]>
## 2 healthyR      <tibble [547 × 2]> <tibble [28 × 2]> <split [519|28]>
## 3 healthyR.ai   <tibble [547 × 2]> <tibble [28 × 2]> <split [519|28]>
## 4 healthyR.data <tibble [547 × 2]> <tibble [28 × 2]> <split [519|28]>
## 5 healthyR.ts   <tibble [542 × 2]> <tibble [28 × 2]> <split [514|28]>
## 6 healthyverse  <tibble [537 × 2]> <tibble [28 × 2]> <split [509|28]>
## 7 tidyAML       <tibble [288 × 2]> <tibble [28 × 2]> <split [260|28]>

Now it is time to make some recipes and models using the modeltime workflow.

Modeltime Workflow

Recipe Object

recipe_base <- recipe(
  value_trans ~ date
  , data = extract_nested_test_split(nested_data_tbl)
  )

recipe_base

recipe_date <- recipe_base %>%
    step_mutate(date = as.numeric(date))

Models

# Models ------------------------------------------------------------------

# Auto ARIMA --------------------------------------------------------------

model_spec_arima_no_boost <- arima_reg() %>%
  set_engine(engine = "auto_arima")

wflw_auto_arima <- workflow() %>%
  add_recipe(recipe = recipe_base) %>%
  add_model(model_spec_arima_no_boost)

# NNETAR ------------------------------------------------------------------

model_spec_nnetar <- nnetar_reg(
  mode              = "regression"
  , seasonal_period = "auto"
) %>%
  set_engine("nnetar")

wflw_nnetar <- workflow() %>%
  add_recipe(recipe = recipe_base) %>%
  add_model(model_spec_nnetar)

# TSLM --------------------------------------------------------------------

model_spec_lm <- linear_reg() %>%
  set_engine("lm")

wflw_lm <- workflow() %>%
  add_recipe(recipe = recipe_base) %>%
  add_model(model_spec_lm)

# MARS --------------------------------------------------------------------

model_spec_mars <- mars(mode = "regression") %>%
  set_engine("earth")

wflw_mars <- workflow() %>%
  add_recipe(recipe = recipe_base) %>%
  add_model(model_spec_mars)

Nested Modeltime Tables

nested_modeltime_tbl <- modeltime_nested_fit(
  # Nested Data
  nested_data = nested_data_tbl,
   control = control_nested_fit(
     verbose = TRUE,
     allow_par = FALSE
   ),
  # Add workflows
  wflw_auto_arima,
  wflw_lm,
  wflw_mars,
  wflw_nnetar
)
nested_modeltime_tbl
## # Nested Modeltime Table
##   # A tibble: 7 × 5
##   package       .actual_data .future_data .splits          .modeltime_tables 
##   <fct>         <list>       <list>       <list>           <list>            
## 1 TidyDensity   <tibble>     <tibble>     <split [522|28]> <mdl_tm_t [4 × 5]>
## 2 healthyR      <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [4 × 5]>
## 3 healthyR.ai   <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [4 × 5]>
## 4 healthyR.data <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [4 × 5]>
## 5 healthyR.ts   <tibble>     <tibble>     <split [514|28]> <mdl_tm_t [4 × 5]>
## 6 healthyverse  <tibble>     <tibble>     <split [509|28]> <mdl_tm_t [4 × 5]>
## 7 tidyAML       <tibble>     <tibble>     <split [260|28]> <mdl_tm_t [4 × 5]>

Model Accuracy

nested_modeltime_tbl %>%
  extract_nested_test_accuracy() %>%
  knitr::kable()
package .model_id .model_desc .type mae mape mase smape rmse rsq
TidyDensity 1 ARIMA Test 0.6796987 505.96108 0.9235829 100.99325 0.8754368 0.1107497
TidyDensity 2 LM Test 0.7348597 470.37200 0.9985364 108.77928 0.9411671 0.2481394
TidyDensity 3 EARTH Test 5.3709826 6145.23181 7.2981573 146.71059 5.9581383 0.2481394
TidyDensity 4 NNAR Test 0.7118461 242.61432 0.9672653 130.28580 0.9413293 0.1724163
healthyR 1 ARIMA Test 0.6091090 119.66201 0.5527547 143.96117 0.7379648 0.5146458
healthyR 2 LM Test 0.7395173 105.02823 0.6710976 198.13336 0.9049509 0.0000225
healthyR 3 EARTH Test 0.7192852 190.77120 0.6527374 146.70618 0.9059985 0.0000225
healthyR 4 NNAR Test 0.6933730 128.00874 0.6292225 160.82319 0.8374558 0.1708856
healthyR.ai 1 ARIMA Test 0.7270730 191.23895 0.6274416 143.88073 0.9151967 0.3653932
healthyR.ai 2 LM Test 0.7710211 191.38627 0.6653675 151.77680 0.9745299 0.0064395
healthyR.ai 3 EARTH Test 0.7608267 248.77298 0.6565701 140.17368 0.9679174 0.0064395
healthyR.ai 4 NNAR Test 0.7459107 156.35336 0.6436980 153.50026 0.9194790 0.1928014
healthyR.data 1 ARIMA Test 0.6421820 245.51445 0.8805174 133.82596 0.7541704 0.0645056
healthyR.data 2 LM Test 0.6695574 96.44137 0.9180527 187.75415 0.7718094 0.0049291
healthyR.data 3 EARTH Test 0.6883625 433.43752 0.9438371 125.32968 0.8189929 0.0049291
healthyR.data 4 NNAR Test 0.6393800 122.64982 0.8766756 168.93479 0.7266487 0.1905260
healthyR.ts 1 ARIMA Test 0.8383359 257.79310 0.7292450 109.65921 1.1484092 0.0539207
healthyR.ts 2 LM Test 0.8368915 338.28930 0.7279886 102.29293 1.1286492 0.0459793
healthyR.ts 3 EARTH Test 1.7787674 699.20727 1.5473000 173.91985 2.1433727 0.0459793
healthyR.ts 4 NNAR Test 0.8719026 146.75143 0.7584438 131.64770 1.2005270 0.1626047
healthyverse 1 ARIMA Test 0.4866996 398.26274 0.6832245 86.09955 0.6311693 0.3488767
healthyverse 2 LM Test 0.5831935 538.47619 0.8186817 92.66371 0.7215907 0.0095300
healthyverse 3 EARTH Test 0.6089132 641.33869 0.8547867 92.24968 0.7423738 0.0095300
healthyverse 4 NNAR Test 0.5627061 392.87718 0.7899216 94.57654 0.7043305 0.0409782
tidyAML 1 ARIMA Test 0.6886549 296.24727 0.8684913 104.13976 0.9314680 0.1884987
tidyAML 2 LM Test 0.8209775 515.22613 1.0353689 100.05065 1.1122287 0.0573557
tidyAML 3 EARTH Test 0.7540684 366.32729 0.9509870 104.58815 0.9970733 0.0573557
tidyAML 4 NNAR Test 0.7391225 253.62481 0.9321381 118.48682 0.9558382 0.0373486

Plot Models

nested_modeltime_tbl %>%
  extract_nested_test_forecast() %>%
  group_by(package) %>%
  plot_modeltime_forecast(
    .interactive = FALSE,
    .conf_interval_show  = FALSE,
    .facet_scales = "free"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Best Model

best_nested_modeltime_tbl <- nested_modeltime_tbl %>%
  modeltime_nested_select_best(
    metric = "rmse",
    minimize = TRUE,
    filter_test_forecasts = TRUE
  )

best_nested_modeltime_tbl %>%
  extract_nested_best_model_report()
## # Nested Modeltime Table
##   # A tibble: 7 × 10
##   package       .model_id .model_desc .type   mae  mape  mase smape  rmse    rsq
##   <fct>             <int> <chr>       <chr> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 TidyDensity           1 ARIMA       Test  0.680  506. 0.924 101.  0.875 0.111 
## 2 healthyR              1 ARIMA       Test  0.609  120. 0.553 144.  0.738 0.515 
## 3 healthyR.ai           1 ARIMA       Test  0.727  191. 0.627 144.  0.915 0.365 
## 4 healthyR.data         4 NNAR        Test  0.639  123. 0.877 169.  0.727 0.191 
## 5 healthyR.ts           2 LM          Test  0.837  338. 0.728 102.  1.13  0.0460
## 6 healthyverse          1 ARIMA       Test  0.487  398. 0.683  86.1 0.631 0.349 
## 7 tidyAML               1 ARIMA       Test  0.689  296. 0.868 104.  0.931 0.188
best_nested_modeltime_tbl %>%
  extract_nested_test_forecast() %>%
  #filter(!is.na(.model_id)) %>%
  group_by(package) %>%
  plot_modeltime_forecast(
    .interactive = FALSE,
    .conf_interval_alpha = 0.2,
    .facet_scales = "free"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Refitting and Future Forecast

Now that we have the best models, we can make our future forecasts.

nested_modeltime_refit_tbl <- best_nested_modeltime_tbl %>%
    modeltime_nested_refit(
        control = control_nested_refit(verbose = TRUE)
    )
nested_modeltime_refit_tbl
## # Nested Modeltime Table
##   # A tibble: 7 × 5
##   package       .actual_data .future_data .splits          .modeltime_tables 
##   <fct>         <list>       <list>       <list>           <list>            
## 1 TidyDensity   <tibble>     <tibble>     <split [522|28]> <mdl_tm_t [1 × 5]>
## 2 healthyR      <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [1 × 5]>
## 3 healthyR.ai   <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [1 × 5]>
## 4 healthyR.data <tibble>     <tibble>     <split [519|28]> <mdl_tm_t [1 × 5]>
## 5 healthyR.ts   <tibble>     <tibble>     <split [514|28]> <mdl_tm_t [1 × 5]>
## 6 healthyverse  <tibble>     <tibble>     <split [509|28]> <mdl_tm_t [1 × 5]>
## 7 tidyAML       <tibble>     <tibble>     <split [260|28]> <mdl_tm_t [1 × 5]>
nested_modeltime_refit_tbl %>%
  extract_nested_future_forecast() %>%
  mutate(across(.value:.conf_hi, .fns = ~ standard_inv_vec(
    x    = .,
    mean = std_mean,
    sd   = std_sd
  )$standard_inverse_value)) %>%
  mutate(across(.value:.conf_hi, .fns = ~ liiv(
    x = .,
    limit_lower = limit_lower,
    limit_upper = limit_upper,
    offset      = offset
  )$rescaled_v)) %>%
  group_by(package) %>%
  plot_modeltime_forecast(
    .interactive = FALSE,
    .conf_interval_alpha = 0.2,
    .facet_scales = "free"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")