healthyverse_tsa
Time Series Analysis, Modeling and Forecasting of the Healthyverse Packages
Steven P. Sanderson II, MPH - Date: 2026-11-04
Introduction
This analysis follows a Nested Modeltime Workflow from modeltime
along with using the NNS package. I use this to monitor the
downloads of all of my packages:
Get Data
glimpse(downloads_tbl)
Rows: 195,087
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 2026-11-02 23:49:59, the file was birthed on: 2025-10-31 10:47:59.603742, and at report knit time is 8817.03 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 | 195087 |
| 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 | 147575 | 0.24 | 1 | 17 | 0 | 56 | 0 |
| r_arch | 147575 | 0.24 | 1 | 7 | 0 | 8 | 0 |
| r_os | 147575 | 0.24 | 6 | 33 | 0 | 43 | 0 |
| package | 0 | 1.00 | 7 | 13 | 0 | 8 | 0 |
| version | 0 | 1.00 | 5 | 17 | 0 | 63 | 0 |
| country | 18740 | 0.90 | 2 | 2 | 0 | 173 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| date | 0 | 1 | 2020-11-23 | 2026-11-02 | 2024-04-10 | 2164 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| size | 0 | 1 | 1165916.34 | 1479414.12 | 355 | 48913.5 | 363140 | 2367799 | 5677952 | ▇▁▂▁▁ |
| ip_id | 0 | 1 | 14000.57 | 30208.35 | 1 | 193.0 | 2860 | 12949 | 429286 | ▇▁▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| date_time | 0 | 1 | 2020-11-23 09:00:41 | 2026-11-02 23:49:59 | 2024-04-10 21:37:24 | 127948 |
Variable type: Timespan
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| time | 0 | 1 | 0 | 59 | 12H 15M 58S | 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.
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
Now lets take a look at some time series decomposition graphs.
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
Seasonal Diagnostics:
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
ACF and PACF Diagnostics:
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
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.
Call:
stats::lm(formula = .formula, data = df)
Residuals:
Min 1Q Median 3Q Max
-148.85 -38.25 -11.67 28.26 828.18
Coefficients:
Estimate Std. Error
(Intercept) -1.369e+02 4.520e+01
date 8.946e-03 2.378e-03
lag(value, 1) 9.954e-02 2.144e-02
lag(value, 7) 7.117e-02 2.203e-02
lag(value, 14) 7.387e-02 2.186e-02
lag(value, 21) 8.861e-02 2.191e-02
lag(value, 28) 7.913e-02 2.185e-02
lag(value, 35) 3.530e-02 2.187e-02
lag(value, 42) 6.137e-02 2.195e-02
lag(value, 49) 6.699e-02 2.188e-02
month(date, label = TRUE).L -7.671e+00 4.655e+00
month(date, label = TRUE).Q 2.268e+00 4.588e+00
month(date, label = TRUE).C -1.526e+01 4.594e+00
month(date, label = TRUE)^4 -1.065e+01 4.619e+00
month(date, label = TRUE)^5 -5.958e+00 4.604e+00
month(date, label = TRUE)^6 2.953e-01 4.604e+00
month(date, label = TRUE)^7 -2.057e+00 4.521e+00
month(date, label = TRUE)^8 -4.297e+00 4.462e+00
month(date, label = TRUE)^9 -5.238e-01 4.451e+00
month(date, label = TRUE)^10 -4.639e-01 4.447e+00
month(date, label = TRUE)^11 -5.717e-01 4.443e+00
fourier_vec(date, type = "sin", K = 1, period = 7) -1.056e+01 2.014e+00
fourier_vec(date, type = "cos", K = 1, period = 7) 7.140e+00 2.068e+00
t value Pr(>|t|)
(Intercept) -3.028 0.002490 **
date 3.762 0.000174 ***
lag(value, 1) 4.642 3.67e-06 ***
lag(value, 7) 3.231 0.001254 **
lag(value, 14) 3.379 0.000740 ***
lag(value, 21) 4.043 5.46e-05 ***
lag(value, 28) 3.621 0.000300 ***
lag(value, 35) 1.614 0.106706
lag(value, 42) 2.796 0.005226 **
lag(value, 49) 3.062 0.002229 **
month(date, label = TRUE).L -1.648 0.099554 .
month(date, label = TRUE).Q 0.494 0.621084
month(date, label = TRUE).C -3.321 0.000913 ***
month(date, label = TRUE)^4 -2.306 0.021207 *
month(date, label = TRUE)^5 -1.294 0.195760
month(date, label = TRUE)^6 0.064 0.948861
month(date, label = TRUE)^7 -0.455 0.649141
month(date, label = TRUE)^8 -0.963 0.335648
month(date, label = TRUE)^9 -0.118 0.906332
month(date, label = TRUE)^10 -0.104 0.916931
month(date, label = TRUE)^11 -0.129 0.897635
fourier_vec(date, type = "sin", K = 1, period = 7) -5.244 1.73e-07 ***
fourier_vec(date, type = "cos", K = 1, period = 7) 3.452 0.000567 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 59.86 on 2092 degrees of freedom
(49 observations deleted due to missingness)
Multiple R-squared: 0.1959, Adjusted R-squared: 0.1875
F-statistic: 23.17 on 22 and 2092 DF, p-value: < 2.2e-16
NNS Forecasting
This is something I have been wanting to try for a while. The NNS
package is a great package for forecasting time series data.
library(NNS)
data_list <- base_data |>
select(package, value) |>
group_split(package)
data_list |>
imap(
\(x, idx) {
obj <- x
x <- obj |> pull(value) |> tail(7*52)
train_set_size <- length(x) - 56
pkg <- obj |> pluck(1) |> unique()
# sf <- NNS.seas(x, modulo = 7, plot = FALSE)$periods
seas <- t(
sapply(
1:25,
function(i) c(
i,
sqrt(
mean((
NNS.ARMA(x,
h = 28,
training.set = train_set_size,
method = "lin",
seasonal.factor = i,
plot=FALSE
) - tail(x, 28)) ^ 2)))
)
)
colnames(seas) <- c("Period", "RMSE")
sf <- seas[which.min(seas[, 2]), 1]
cat(paste0("Package: ", pkg, "\n"))
NNS.ARMA.optim(
variable = x,
h = 28,
training.set = train_set_size,
#seasonal.factor = seq(12, 60, 7),
seasonal.factor = sf,
pred.int = 0.95,
plot = TRUE
)
title(
sub = paste0("\n",
"Package: ", pkg, " - NNS Optimization")
)
}
)
Package: healthyR
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 18 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 12.6844310262358"
[1] "BEST method = 'lin' PATH MEMBER = c( 18 )"
[1] "BEST lin OBJECTIVE FUNCTION = 12.6844310262358"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 18 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 3.28840603457942"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 18 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 3.28840603457942"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 18 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 4.29900717145417"
[1] "BEST method = 'both' PATH MEMBER = c( 18 )"
[1] "BEST both OBJECTIVE FUNCTION = 4.29900717145417"
Package: healthyR.ai
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 73.3393187220088"
[1] "BEST method = 'lin' PATH MEMBER = c( 2 )"
[1] "BEST lin OBJECTIVE FUNCTION = 73.3393187220088"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 4.73196669246817"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 2 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 4.73196669246817"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 7.98217618960141"
[1] "BEST method = 'both' PATH MEMBER = c( 2 )"
[1] "BEST both OBJECTIVE FUNCTION = 7.98217618960141"
Package: healthyR.data
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 4.18687951464587"
[1] "BEST method = 'lin' PATH MEMBER = c( 24 )"
[1] "BEST lin OBJECTIVE FUNCTION = 4.18687951464587"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 2.71506789212433"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 24 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 2.71506789212433"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 2.56340333683306"
[1] "BEST method = 'both' PATH MEMBER = c( 24 )"
[1] "BEST both OBJECTIVE FUNCTION = 2.56340333683306"
Package: healthyR.ts
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 16 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 30.5183213997252"
[1] "BEST method = 'lin' PATH MEMBER = c( 16 )"
[1] "BEST lin OBJECTIVE FUNCTION = 30.5183213997252"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 16 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 11.8090474821482"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 16 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 11.8090474821482"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 16 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 19.5898995848496"
[1] "BEST method = 'both' PATH MEMBER = c( 16 )"
[1] "BEST both OBJECTIVE FUNCTION = 19.5898995848496"
Package: healthyverse
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 31.8115996188069"
[1] "BEST method = 'lin' PATH MEMBER = c( 24 )"
[1] "BEST lin OBJECTIVE FUNCTION = 31.8115996188069"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 23.2245557981719"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 24 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 23.2245557981719"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 24 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 36.1694605526754"
[1] "BEST method = 'both' PATH MEMBER = c( 24 )"
[1] "BEST both OBJECTIVE FUNCTION = 36.1694605526754"
Package: RandomWalker
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 6.32437296292934"
[1] "BEST method = 'lin' PATH MEMBER = c( 2 )"
[1] "BEST lin OBJECTIVE FUNCTION = 6.32437296292934"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 4.58180947642768"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 2 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 4.58180947642768"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 2 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 4.99621406397537"
[1] "BEST method = 'both' PATH MEMBER = c( 2 )"
[1] "BEST both OBJECTIVE FUNCTION = 4.99621406397537"
Package: tidyAML
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 5 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 18.7953407949233"
[1] "BEST method = 'lin' PATH MEMBER = c( 5 )"
[1] "BEST lin OBJECTIVE FUNCTION = 18.7953407949233"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 5 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 7.0124900020093"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 5 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 7.0124900020093"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 5 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 10.1009816823077"
[1] "BEST method = 'both' PATH MEMBER = c( 5 )"
[1] "BEST both OBJECTIVE FUNCTION = 10.1009816823077"
Package: TidyDensity
[1] "CURRNET METHOD: lin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 15 ) ...)"
[1] "CURRENT lin OBJECTIVE FUNCTION = 13.0701171916155"
[1] "BEST method = 'lin' PATH MEMBER = c( 15 )"
[1] "BEST lin OBJECTIVE FUNCTION = 13.0701171916155"
[1] "CURRNET METHOD: nonlin"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 15 ) ...)"
[1] "CURRENT nonlin OBJECTIVE FUNCTION = 13.8826845223752"
[1] "BEST method = 'nonlin' PATH MEMBER = c( 15 )"
[1] "BEST nonlin OBJECTIVE FUNCTION = 13.8826845223752"
[1] "CURRNET METHOD: both"
[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 15 ) ...)"
[1] "CURRENT both OBJECTIVE FUNCTION = 13.9582021727639"
[1] "BEST method = 'both' PATH MEMBER = c( 15 )"
[1] "BEST both OBJECTIVE FUNCTION = 13.9582021727639"
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULL
[[5]]
NULL
[[6]]
NULL
[[7]]
NULL
[[8]]
NULL
Pre-Processing
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 %>%
group_by(package) %>%
# get standardization
mutate(value_trans = standard_vec(value_trans, silent = TRUE)$standard_scaled) %>%
tk_augment_fourier(
.date_var = date,
.periods = c(7, 14, 30, 90, 180),
.K = 2
) %>%
tk_augment_timeseries_signature(
.date_var = date
) %>%
ungroup() %>%
select(-c(value, -year.iso))
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 %>%
# 0. Filter out column where package is NA
filter(!is.na(package)) %>%
# 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: 8 × 4
package .actual_data .future_data .splits
<fct> <list> <list> <list>
1 healthyR.data <tibble [2,152 × 50]> <tibble [28 × 50]> <split [2124|28]>
2 healthyR <tibble [2,146 × 50]> <tibble [28 × 50]> <split [2118|28]>
3 healthyR.ts <tibble [2,082 × 50]> <tibble [28 × 50]> <split [2054|28]>
4 healthyverse <tibble [1,951 × 50]> <tibble [28 × 50]> <split [1923|28]>
5 healthyR.ai <tibble [1,887 × 50]> <tibble [28 × 50]> <split [1859|28]>
6 TidyDensity <tibble [1,740 × 50]> <tibble [28 × 50]> <split [1712|28]>
7 tidyAML <tibble [1,343 × 50]> <tibble [28 × 50]> <split [1315|28]>
8 RandomWalker <tibble [767 × 50]> <tibble [28 × 50]> <split [739|28]>
Now it is time to make some recipes and models using the modeltime workflow.
Modeltime Workflow
Recipe Object
recipe_base <- recipe(
value_trans ~ .
, data = extract_nested_test_split(nested_data_tbl)
)
recipe_base
recipe_date <- recipe(
value_trans ~ date
, data = extract_nested_test_split(nested_data_tbl)
)
Models
# Models ------------------------------------------------------------------
# Auto ARIMA --------------------------------------------------------------
model_spec_arima_no_boost <- arima_reg() %>%
set_engine(engine = "auto_arima")
wflw_auto_arima <- workflow() %>%
add_recipe(recipe = recipe_date) %>%
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_date) %>%
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_tbl[!is.na(nested_modeltime_tbl$package),]
Model Accuracy
nested_modeltime_tbl %>%
extract_nested_test_accuracy() %>%
filter(!is.na(package)) %>%
knitr::kable()
| package | .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
|---|---|---|---|---|---|---|---|---|---|
| healthyR.data | 1 | ARIMA | Test | 0.6390046 | 145.41843 | 0.6619628 | 114.96847 | 0.8000284 | 0.0003527 |
| healthyR.data | 2 | LM | Test | 0.9987961 | 196.63237 | 1.0346809 | 170.43429 | 1.1923765 | 0.0716657 |
| healthyR.data | 3 | EARTH | Test | 0.8831485 | 424.70999 | 0.9148784 | 101.14699 | 1.0575506 | 0.0171136 |
| healthyR.data | 4 | NNAR | Test | 0.9270897 | 131.87687 | 0.9603983 | 155.22161 | 1.1399700 | 0.0085174 |
| healthyR | 1 | ARIMA | Test | 0.6236056 | 188.82459 | 0.8547425 | 151.78291 | 0.8261691 | 0.0531245 |
| healthyR | 2 | LM | Test | 0.8355000 | 314.99298 | 1.1451747 | 153.43705 | 1.0286232 | 0.0046704 |
| healthyR | 3 | EARTH | Test | 1.0227010 | 974.46513 | 1.4017611 | 130.99040 | 1.1649737 | 0.0000241 |
| healthyR | 4 | NNAR | Test | 0.6419218 | 223.94837 | 0.8798476 | 145.55025 | 0.8381883 | 0.0018528 |
| healthyR.ts | 1 | ARIMA | Test | 0.7864207 | 260.36248 | 0.9324022 | 167.68897 | 1.0225103 | 0.0742807 |
| healthyR.ts | 2 | LM | Test | 0.7535258 | 210.08026 | 0.8934011 | 148.78125 | 0.9668758 | 0.0170318 |
| healthyR.ts | 3 | EARTH | Test | 0.7613605 | 626.67628 | 0.9026901 | 106.74086 | 0.9088390 | 0.0001175 |
| healthyR.ts | 4 | NNAR | Test | 0.6709391 | 257.48147 | 0.7954840 | 146.22297 | 0.8668881 | 0.0296466 |
| healthyverse | 1 | ARIMA | Test | 0.7498031 | 152.81587 | 0.8511218 | 75.29417 | 0.9266849 | 0.0218056 |
| healthyverse | 2 | LM | Test | 0.7239930 | 100.70878 | 0.8218240 | 88.36626 | 0.9379998 | 0.0012007 |
| healthyverse | 3 | EARTH | Test | 0.8701898 | 82.25926 | 0.9877760 | 118.79077 | 1.0522450 | 0.0024569 |
| healthyverse | 4 | NNAR | Test | 0.9744019 | 110.62901 | 1.1060699 | 135.35760 | 1.1509197 | 0.0812302 |
| healthyR.ai | 1 | ARIMA | Test | 0.6420321 | 125.11823 | 0.7311706 | 113.14446 | 0.7750812 | 0.0777528 |
| healthyR.ai | 2 | LM | Test | 0.9568014 | 186.34010 | 1.0896419 | 164.06869 | 1.1179020 | 0.0333757 |
| healthyR.ai | 3 | EARTH | Test | 0.8461858 | 341.11263 | 0.9636685 | 92.03720 | 1.0295793 | 0.0007256 |
| healthyR.ai | 4 | NNAR | Test | 0.8330621 | 141.95997 | 0.9487228 | 162.86880 | 0.9954640 | 0.0055734 |
| TidyDensity | 1 | ARIMA | Test | 0.6715099 | 131.01341 | 0.8298905 | 185.49061 | 0.7932926 | 0.0830199 |
| TidyDensity | 2 | LM | Test | 0.5798722 | 211.77380 | 0.7166394 | 122.81124 | 0.7335301 | 0.1103472 |
| TidyDensity | 3 | EARTH | Test | 0.6355769 | 129.08218 | 0.7854824 | 156.91961 | 0.7559654 | 0.1120048 |
| TidyDensity | 4 | NNAR | Test | 0.7428078 | 337.81615 | 0.9180045 | 134.20181 | 0.8829575 | 0.0244800 |
| tidyAML | 1 | ARIMA | Test | 0.7456235 | 151.60214 | 1.0308694 | 106.66637 | 0.9148184 | 0.0387333 |
| tidyAML | 2 | LM | Test | 0.9328374 | 153.51249 | 1.2897039 | 164.29231 | 1.0919806 | 0.0023960 |
| tidyAML | 3 | EARTH | Test | 0.7516625 | 209.58986 | 1.0392187 | 95.87034 | 0.9548689 | 0.3651713 |
| tidyAML | 4 | NNAR | Test | 0.8142920 | 125.26801 | 1.1258078 | 141.89858 | 1.0135785 | 0.0044530 |
| RandomWalker | 1 | ARIMA | Test | 0.5447159 | 131.84374 | 0.7573790 | 172.67069 | 0.7025051 | 0.1567579 |
| RandomWalker | 2 | LM | Test | 0.4626852 | 153.96640 | 0.6433226 | 128.54087 | 0.6043155 | 0.0132084 |
| RandomWalker | 3 | EARTH | Test | 0.4826263 | 99.70126 | 0.6710490 | 170.78898 | 0.6261649 | 0.0338361 |
| RandomWalker | 4 | NNAR | Test | 0.5388196 | 201.56788 | 0.7491808 | 141.96181 | 0.6785070 | 0.0025156 |
Plot Models
nested_modeltime_tbl %>%
extract_nested_test_forecast() %>%
group_by(package) %>%
filter_by_time(.date_var = .index, .start_date = max(.index) - 60) %>%
ungroup() %>%
plot_modeltime_forecast(
.interactive = FALSE,
.conf_interval_show = FALSE,
.facet_scales = "free"
) +
theme_minimal() +
facet_wrap(~ package, nrow = 3) +
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: 8 × 10
package .model_id .model_desc .type mae mape mase smape rmse rsq
<fct> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 healthyR.da… 1 ARIMA Test 0.639 145. 0.662 115. 0.800 3.53e-4
2 healthyR 1 ARIMA Test 0.624 189. 0.855 152. 0.826 5.31e-2
3 healthyR.ts 4 NNAR Test 0.671 257. 0.795 146. 0.867 2.96e-2
4 healthyverse 1 ARIMA Test 0.750 153. 0.851 75.3 0.927 2.18e-2
5 healthyR.ai 1 ARIMA Test 0.642 125. 0.731 113. 0.775 7.78e-2
6 TidyDensity 2 LM Test 0.580 212. 0.717 123. 0.734 1.10e-1
7 tidyAML 1 ARIMA Test 0.746 152. 1.03 107. 0.915 3.87e-2
8 RandomWalker 2 LM Test 0.463 154. 0.643 129. 0.604 1.32e-2
best_nested_modeltime_tbl %>%
extract_nested_test_forecast() %>%
#filter(!is.na(.model_id)) %>%
group_by(package) %>%
filter_by_time(.date_var = .index, .start_date = max(.index) - 60) %>%
ungroup() %>%
plot_modeltime_forecast(
.interactive = FALSE,
.conf_interval_alpha = 0.2,
.facet_scales = "free"
) +
facet_wrap(~ package, nrow = 3) +
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: 8 × 5
package .actual_data .future_data .splits .modeltime_tables
<fct> <list> <list> <list> <list>
1 healthyR.data <tibble> <tibble> <split [2124|28]> <mdl_tm_t [1 × 5]>
2 healthyR <tibble> <tibble> <split [2118|28]> <mdl_tm_t [1 × 5]>
3 healthyR.ts <tibble> <tibble> <split [2054|28]> <mdl_tm_t [1 × 5]>
4 healthyverse <tibble> <tibble> <split [1923|28]> <mdl_tm_t [1 × 5]>
5 healthyR.ai <tibble> <tibble> <split [1859|28]> <mdl_tm_t [1 × 5]>
6 TidyDensity <tibble> <tibble> <split [1712|28]> <mdl_tm_t [1 × 5]>
7 tidyAML <tibble> <tibble> <split [1315|28]> <mdl_tm_t [1 × 5]>
8 RandomWalker <tibble> <tibble> <split [739|28]> <mdl_tm_t [1 × 5]>
nested_modeltime_refit_tbl %>%
extract_nested_future_forecast() %>%
group_by(package) %>%
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)) %>%
filter_by_time(.date_var = .index, .start_date = max(.index) - 60) %>%
ungroup() %>%
plot_modeltime_forecast(
.interactive = FALSE,
.conf_interval_alpha = 0.2,
.facet_scales = "free"
) +
facet_wrap(~ package, nrow = 3) +
theme_minimal() +
theme(legend.position = "bottom")
