`forecast.model_combination()` fails for combinations containing a `prophet()` model

[mpjashby]

I'm not sure if this is an issue with fablelite or fable.prophet – I think (from the traceback) it's in fablelite:::forecast.model_combination(), hence opening the issue here.

Combination models including a prophet model can be generated, but calling forecast() on the result leads to an error.

library("fable")
#> Loading required package: fablelite
library("fable.prophet")
#> Loading required package: Rcpp
library("lubridate")
#> 
#> Attaching package: 'lubridate'
#> The following object is masked from 'package:base':
#> 
#>     date
library("tidyverse")

# generate some models based on the first month of `vic_elec` data, using an
# arima and tslm model for comparison
models <- tsibbledata::vic_elec %>% 
    filter(between(Date, ymd("2012-01-01"), ymd("2012-01-31"))) %>% 
    model(
        arima = ARIMA(Demand ~ trend() + season() + Temperature + Holiday),
        prophet = prophet(Demand ~ growth() + season('week') + Temperature + Holiday),
        tslm = TSLM(Demand ~ trend() + season() + Temperature + Holiday),
        combo = combination_model(
            ARIMA(Demand ~ trend() + season() + Temperature + Holiday),
            prophet(Demand ~ growth() + season('week') + Temperature + Holiday),
            TSLM(Demand ~ trend() + season() + Temperature + Holiday)
        )
    )

# check the models are not NULL etc.
report(models$combo[[1]])
#> Series: Demand 
#> Model: COMBINATION 
#> Combination: (Demand + Demand + Demand) * 0.333333333333333
#> 
#> ===========================================================
#> 
#> Series: Demand + Demand + Demand 
#> Model: COMBINATION 
#> Combination: Demand + Demand + Demand
#> 
#> =====================================
#> 
#> Series: Demand + Demand 
#> Model: COMBINATION 
#> Combination: Demand + Demand
#> 
#> ============================
#> 
#> Series: Demand 
#> Model: LM w/ ARIMA(2,0,4) errors 
#> 
#> Coefficients:
#>          ar1      ar2      ma1      ma2      ma3     ma4  trend()
#>       1.9013  -0.9172  -0.2509  -0.3299  -0.1035  0.2574   0.4911
#> s.e.  0.0205   0.0199   0.0335   0.0322   0.0302  0.0260   0.1864
#>       season()  Temperature  Holiday  intercept
#>         9.6897      11.0670   3.9816  4259.0318
#> s.e.    1.5002       6.4805  39.7646   210.9456
#> 
#> sigma^2 estimated as 7479:  log likelihood=-8745.27
#> AIC=17514.55   AICc=17514.76   BIC=17578.21
#> 
#> Series: Demand 
#> Model: prophet 
#> 
#> A model specific report is not available for this model class.
#> 
#> Series: Demand 
#> Model: TSLM 
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1420.1  -506.8    11.1   474.1  1610.1 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 1573.12374   76.13591   20.66  < 2e-16 ***
#> trend()        0.28884    0.04108    7.03 3.13e-12 ***
#> season()       9.86739   33.98112    0.29    0.772    
#> Temperature  143.20583    3.23056   44.33  < 2e-16 ***
#> HolidayTRUE -547.48679   60.83494   -9.00  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 655.4 on 1483 degrees of freedom
#> Multiple R-squared: 0.5918,  Adjusted R-squared: 0.5907
#> F-statistic: 537.6 on 4 and 1483 DF, p-value: < 2.22e-16

# use the second month of data for forecasts
new_data <- tsibbledata::vic_elec %>% 
    filter(between(Date, ymd("2012-02-01"), ymd("2012-02-28")))

# generate combination models and forecast
models %>% 
    select(combo) %>% 
    forecast(new_data = new_data) %>% 
    head()
#> Error in vapply(..., USE.NAMES = FALSE): values must be length 1,
#>  but FUN(X[[1]]) result is length 0

models %>% 
    mutate(combo = (arima + tslm) / 2) %>% 
    forecast(new_data = new_data) %>% 
    head()
#> # A fable: 6 x 7 [30m] <UTC>
#> # Key:     .model [1]
#>   .model Time                Demand .distribution Temperature Date      
#>   <chr>  <dttm>               <dbl> <dist>              <dbl> <date>    
#> 1 arima  2012-02-01 00:00:00  4353. N(4353,   74…        15.5 2012-02-01
#> 2 arima  2012-02-01 00:30:00  4372. N(4372,  278…        15.3 2012-02-01
#> 3 arima  2012-02-01 01:00:00  4367. N(4367,  545…        15.2 2012-02-01
#> 4 arima  2012-02-01 01:30:00  4447. N(4447,  838…        15.1 2012-02-01
#> 5 arima  2012-02-01 02:00:00  4512. N(4512, 1228…        15.0 2012-02-01
#> 6 arima  2012-02-01 02:30:00  4601. N(4601, 1706…        14.8 2012-02-01
#> # … with 1 more variable: Holiday <lgl>

models %>% 
    mutate(combo = (arima + prophet) / 2) %>% 
    forecast(new_data = new_data) %>% 
    head()
#> Error in vapply(..., USE.NAMES = FALSE): values must be length 1,
#>  but FUN(X[[1]]) result is length 0

models %>% 
    mutate(combo = (prophet + tslm) / 2) %>% 
    forecast(new_data = new_data) %>% 
    head()
#> Error in vapply(..., USE.NAMES = FALSE): values must be length 1,
#>  but FUN(X[[1]]) result is length 0

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

[mitchelloharawild]

Thanks for the question.

Currently only combinations of normally distributed forecasts are supported, although I'd like to generalise from this in a future version. The prophet package produces intervals from simulated future paths, which are not necessarily normally distributed.

The error message is confusing though, so I'll make this better. :+1: