Агрегирование прогнозов с помощью Fable

Question

Проблема: с помощью Fable я могу легко создавать прогнозы для временных рядов с сгруппированной структурой и даже использовать Fable.aggregate_key/ reconcileсинтаксис для создания согласованного прогноза верхнего уровня. Однако я не могу легко получить доступ к сводным прогнозам с помощью этого метода, и альтернатива, которую я использую, включает отказ от структуры басни (таблицы прогнозов). Может ли кто-нибудь сказать мне, есть ли более простой / предполагаемый способ сделать это с помощью пакета? Как вы можете видеть в примерах, я могу добраться туда другими методами, но я хотел бы знать, есть ли способ лучше. Любая помощь с благодарностью получена!

Подход 1. Мои попытки обобщить прогноз без использованияaggregate_key/ reconcile в основном использовали dplyr group_by а также summarise, однако интервал прогнозирования для прогноза отформатирован как объект нормального распределения, который, похоже, не поддерживает суммирование с использованием этого метода. Чтобы обойти это, я использовалhilo а также unpack_hiloдля извлечения границ для различных интервалов прогнозирования, которые затем можно суммировать обычным методом. Однако мне бы очень хотелось сохранить структуру басни и объекты раздачи, что невозможно при использовании этого метода.

Подход 2: Альтернатива, использующаяaggregate_key/ reconcile только кажется, поддерживает агрегацию с использованием min_trace. Я понимаю, что этот метод предназначен для оптимального согласования, тогда как мне нужен простой агрегированный прогноз снизу вверх. Кажется, что должен быть простой способ получения прогнозов снизу вверх с использованием этого синтаксиса, но я пока его не нашел. Более того, даже используяmin_trace Я не уверен, как получить доступ к самому агрегированному прогнозу, как вы можете видеть в примере!

Пример использования подхода 1:

library(fable)
#> Loading required package: fabletools
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

lung_deaths_agg <- as_tsibble(cbind(mdeaths, fdeaths))
  
fc_1 <- lung_deaths_agg %>% 
  model(lm = TSLM(value ~ trend() + season())) %>% 
  forecast()

fc_1
#> # A fable: 48 x 5 [1M]
#> # Key:     key, .model [2]
#>    key     .model    index        value .mean
#>    <chr>   <chr>     <mth>       <dist> <dbl>
#>  1 fdeaths lm     1980 Jan N(794, 5940)  794.
#>  2 fdeaths lm     1980 Feb N(778, 5940)  778.
#>  3 fdeaths lm     1980 Mar N(737, 5940)  737.
#>  4 fdeaths lm     1980 Apr N(577, 5940)  577.
#>  5 fdeaths lm     1980 May N(456, 5940)  456.
#>  6 fdeaths lm     1980 Jun N(386, 5940)  386.
#>  7 fdeaths lm     1980 Jul N(379, 5940)  379.
#>  8 fdeaths lm     1980 Aug N(335, 5940)  335.
#>  9 fdeaths lm     1980 Sep N(340, 5940)  340.
#> 10 fdeaths lm     1980 Oct N(413, 5940)  413.
#> # ... with 38 more rows

fc_1 %>%
  hilo() %>% 
  unpack_hilo(c(`80%`, `95%`)) %>% 
  as_tibble() %>% 
  group_by(index) %>% 
  summarise(across(c(.mean, ends_with("upper"), ends_with("lower")), sum))
#> `summarise()` ungrouping output (override with `.groups` argument)
#> # A tibble: 24 x 6
#>       index .mean `80%_upper` `95%_upper` `80%_lower` `95%_lower`
#>       <mth> <dbl>       <dbl>       <dbl>       <dbl>       <dbl>
#>  1 1980 Jan 2751.       3089.       3267.       2414.       2236.
#>  2 1980 Feb 2687.       3024.       3202.       2350.       2171.
#>  3 1980 Mar 2535.       2872.       3051.       2198.       2020.
#>  4 1980 Apr 2062.       2399.       2577.       1725.       1546.
#>  5 1980 May 1597.       1934.       2113.       1260.       1082.
#>  6 1980 Jun 1401.       1738.       1916.       1064.        885.
#>  7 1980 Jul 1343.       1680.       1858.       1006.        827.
#>  8 1980 Aug 1200.       1538.       1716.        863.        685.
#>  9 1980 Sep 1189.       1527.       1705.        852.        674.
#> 10 1980 Oct 1482.       1819.       1998.       1145.        967.
#> # ... with 14 more rows

Пример использования подхода 2:

fc_2 <- lung_deaths_agg %>%
  aggregate_key(key, value = sum(value)) %>% 
  model(lm = TSLM(value ~ trend() + season())) %>%
  reconcile(lm = min_trace(lm)) %>% 
  forecast()

fc_2
#> # A fable: 72 x 5 [1M]
#> # Key:     key, .model [3]
#>    key     .model    index        value .mean
#>    <chr>   <chr>     <mth>       <dist> <dbl>
#>  1 fdeaths lm     1980 Jan N(794, 5606)  794.
#>  2 fdeaths lm     1980 Feb N(778, 5606)  778.
#>  3 fdeaths lm     1980 Mar N(737, 5606)  737.
#>  4 fdeaths lm     1980 Apr N(577, 5606)  577.
#>  5 fdeaths lm     1980 May N(456, 5606)  456.
#>  6 fdeaths lm     1980 Jun N(386, 5606)  386.
#>  7 fdeaths lm     1980 Jul N(379, 5606)  379.
#>  8 fdeaths lm     1980 Aug N(335, 5606)  335.
#>  9 fdeaths lm     1980 Sep N(340, 5606)  340.
#> 10 fdeaths lm     1980 Oct N(413, 5606)  413.
#> # ... with 62 more rows

fc_2 %>% as_tibble() %>% select(key) %>% slice(50:55)
#> # A tibble: 6 x 1
#>   key         
#>   <chr>       
#> 1 <aggregated>
#> 2 <aggregated>
#> 3 <aggregated>
#> 4 <aggregated>
#> 5 <aggregated>
#> 6 <aggregated>

fc_2 %>% as_tibble() %>% select(key) %>% filter(key == "<aggregated>")
#> # A tibble: 0 x 1
#> # ... with 1 variable: key <chr>

Answer 1

Approach 1:

Working with distributions requires more care (than numbers) when adding things together. More specifically, the mean of a Normal distribution can be added without issue:

library(distributional)
mean(dist_normal(2,3) + dist_normal(4,1))
#> [1] 6
mean(dist_normal(2,3)) + mean(dist_normal(4,1))
#> [1] 6

^{Created on 2020-07-03 by the reprex package (v0.3.0)}

However the quantiles (used to produce your 80% and 95% intervals) cannot:

library(distributional)
quantile(dist_normal(2,3) + dist_normal(4,1), 0.9)
#> [1] 10.05262
quantile(dist_normal(2,3), 0.9) + quantile(dist_normal(4,1), 0.9)
#> [1] 11.12621

^{Created on 2020-07-03 by the reprex package (v0.3.0)}

If you want to aggregate distributions, you'll need to compute the sum on the distribution itself:

library(fable)
library(dplyr)
lung_deaths_agg <- as_tsibble(cbind(mdeaths, fdeaths))

fc_1 <- lung_deaths_agg %>% 
  model(lm = fable::TSLM(value ~ trend() + season())) %>% 
  forecast()
fc_1 %>% 
  summarise(value = sum(value), .mean = mean(value))
#> # A fable: 24 x 3 [1M]
#>       index          value .mean
#>       <mth>         <dist> <dbl>
#>  1 1980 Jan N(2751, 40520) 2751.
#>  2 1980 Feb N(2687, 40520) 2687.
#>  3 1980 Mar N(2535, 40520) 2535.
#>  4 1980 Apr N(2062, 40520) 2062.
#>  5 1980 May N(1597, 40520) 1597.
#>  6 1980 Jun N(1401, 40520) 1401.
#>  7 1980 Jul N(1343, 40520) 1343.
#>  8 1980 Aug N(1200, 40520) 1200.
#>  9 1980 Sep N(1189, 40520) 1189.
#> 10 1980 Oct N(1482, 40520) 1482.
#> # … with 14 more rows

^{Created on 2020-07-03 by the reprex package (v0.3.0)}

Note that this will require the development versions of fabletools (>=0.2.0.9000) and distributional (>=0.1.0.9000) as I have added new features to make this example work.

Approach 2:

Experimental support for bottom up reconciliation is available using fabletools:::bottom_up(). This is currently an internal function as I'm still working on some details of how reconciliation can be done more generally in fabletools.

Matching aggregated values should be done with is_aggregated().

fc_2 <- lung_deaths_agg %>%
  aggregate_key(key, value = sum(value)) %>% 
  model(lm = TSLM(value ~ trend() + season())) %>%
  reconcile(lm = min_trace(lm)) %>% 
  forecast()

fc_2 %>% 
  filter(is_aggregated(key))
#> # A fable: 24 x 5 [1M]
#> # Key:     key, .model [1]
#>    key          .model    index          value .mean
#>    <chr>        <chr>     <mth>         <dist> <dbl>
#>  1 <aggregated> lm     1980 Jan N(2751, 24989) 2751.
#>  2 <aggregated> lm     1980 Feb N(2687, 24989) 2687.
#>  3 <aggregated> lm     1980 Mar N(2535, 24989) 2535.
#>  4 <aggregated> lm     1980 Apr N(2062, 24989) 2062.
#>  5 <aggregated> lm     1980 May N(1597, 24989) 1597.
#>  6 <aggregated> lm     1980 Jun N(1401, 24989) 1401.
#>  7 <aggregated> lm     1980 Jul N(1343, 24989) 1343.
#>  8 <aggregated> lm     1980 Aug N(1200, 24989) 1200.
#>  9 <aggregated> lm     1980 Sep N(1189, 24989) 1189.
#> 10 <aggregated> lm     1980 Oct N(1482, 24989) 1482.
#> # … with 14 more rows

^{Created on 2020-07-03 by the reprex package (v0.3.0)}

Comparing an aggregated vector with "<aggregated>" is ambiguous, as your key's character value may be "<aggregated>" without the value being <aggregated>. I've now updated fabletools to match "<aggregated>" with aggregated values with a warning and hint, so this code now gives:

fc_2 %>% 
  filter(key == "<aggregated>")
#> Warning: <aggregated> character values have been converted to aggregated values.
#> Hint: If you're trying to compare aggregated values, use `is_aggregated()`.
#> # A fable: 24 x 5 [1M]
#> # Key:     key, .model [1]
#>    key          .model    index          value .mean
#>    <chr>        <chr>     <mth>         <dist> <dbl>
#>  1 <aggregated> lm     1980 Jan N(2751, 24989) 2751.
#>  2 <aggregated> lm     1980 Feb N(2687, 24989) 2687.
#>  3 <aggregated> lm     1980 Mar N(2535, 24989) 2535.
#>  4 <aggregated> lm     1980 Apr N(2062, 24989) 2062.
#>  5 <aggregated> lm     1980 May N(1597, 24989) 1597.
#>  6 <aggregated> lm     1980 Jun N(1401, 24989) 1401.
#>  7 <aggregated> lm     1980 Jul N(1343, 24989) 1343.
#>  8 <aggregated> lm     1980 Aug N(1200, 24989) 1200.
#>  9 <aggregated> lm     1980 Sep N(1189, 24989) 1189.
#> 10 <aggregated> lm     1980 Oct N(1482, 24989) 1482.
#> # … with 14 more rows

^{Created on 2020-07-03 by the reprex package (v0.3.0)}