TutorialMarkdown.Rmd

---
title: "Uncertainty Tutorial"
output: html_document
date: "2023-06-22"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#library(remotes)
#remotes::install_github("https://github.com/eco4cast/neon4cast")
library(tidyverse)
library(neon4cast)
library(lubridate)
library(rMR)
library(glue)
library(tsibble)
library(fable)
library(arrow)
library(forecast)
library(here)
library(sensitivity)
library(mvtnorm)

library(ranger) #needed for random forest implementation
library(doParallel)
library(tidymodels)
library(vip)
library(butcher)
library(bundle)
library(neon4cast)
library(lubridate)
library(decor)


Nproc <- 8
```

# Forecast Uncertainty Tutorial: Using a forecast of Phenology as an example

This tutorial will cover information on uncertainty analyses for forecasting models, as well as an application of methods to two models of phenology, one parametric and the other a machine learning model.

## The Models

The models we used are adapted from the EFI Theory Groups automated forecasts, which can be found here: LINK. We used an ARIMA model as our parametric model, and a Random Forest model as our machine learning model.

For both models, we will need to load in the data from NEON to make the forecasts.

### ARIMA Model (parametric)

First, we load in the required data, including driver data. Then, we will format the data.

```{r}
download_target <- function(theme = c("aquatics", "beetles",
                                      "phenology", "terrestrial_30min",
                                      "terrestrial_daily","ticks")){  
  theme <- match.arg(theme)
  
  target_file <- switch(theme,
                        aquatics = "aquatics-targets.csv.gz",
                        beetles = "beetles-targets.csv.gz",
                        phenology = "phenology-targets.csv.gz",
                        terrestrial_daily = "terrestrial_daily-targets.csv.gz",
                        terrestrial_30min = "terrestrial_30min-targets.csv.gz",
                        ticks = "ticks-targets.csv.gz"
  )
  download_url <- paste0("https://data.ecoforecast.org/neon4cast-targets/",
                         theme, "/", target_file)
  
  readr::read_csv(download_url, show_col_types = FALSE,
                  lazy = FALSE, progress = FALSE)#%>% 
  #as_tibble(index=time, key=siteID)
}
target = download_target(theme="phenology")
target = target |> filter(site_id == "HARV", variable == "gcc_90")

#### Step 1: Define team name, team members, and theme
model_id = "tg_arima"
model_themes = c("phenology") #This model is only relevant for three themes. I am registered for all three
model_types = c("phenology")

#### Step 2: Get NOAA driver data
forecast_date <- as.Date("2023-05-04")
noaa_date <- forecast_date - lubridate::days(1)  #Need to use yesterday's NOAA forecast because today's is not available yet

#We're going to get data for all sites relevant to this model, so as to not have to re-load data for the same sites
site_data <- readr::read_csv("https://raw.githubusercontent.com/eco4cast/neon4cast-targets/main/NEON_Field_Site_Metadata_20220412.csv") %>%
  filter(if_any(matches(model_types),~.==1))
all_sites = site_data$field_site_id

# specify meteorological variables needed to make predictions
variables <- c('air_temperature',
               "surface_downwelling_shortwave_flux_in_air",
               "precipitation_flux",
               "relative_humidity")

# Load stage 2 data
endpoint = "data.ecoforecast.org"
use_bucket <- paste0("neon4cast-drivers/noaa/gefs-v12/stage2/parquet/0/", noaa_date)
use_s3 <- arrow::s3_bucket(use_bucket, endpoint_override = endpoint, anonymous = TRUE)
noaa_future <- arrow::open_dataset(use_s3) |>
  dplyr::filter(site_id %in% 'HARV',
                datetime >= forecast_date,
                #                reference_datetime == lubridate::as_datetime(forecast_date),
                variable %in% variables) |>
  dplyr::collect()

# Format met forecasts
noaa_future_daily <- noaa_future |> 
  mutate(datetime = lubridate::as_date(datetime)) |> 
  # mean daily forecasts at each site per ensemble
  group_by(datetime, parameter, variable) |> 
  summarize(prediction = mean(prediction)) |>
  pivot_wider(names_from = variable, values_from = prediction) |>
  # convert to Celsius
  mutate(air_temperature = air_temperature - 273.15) |> 
  select(datetime, all_of(variables), parameter)

## grab past met data
met = neon4cast::noaa_stage3() |>
  filter(site_id == "HARV",variable %in% variables,datetime < forecast_date) |>
  collect()
met_daily = met |>
  mutate(datetime = lubridate::as_date(datetime)) |> 
  # mean daily forecasts at each site per ensemble
  group_by(datetime, variable) |> 
  summarize(prediction = mean(prediction)) |>
  pivot_wider(names_from = variable, values_from = prediction) |>
  # convert to Celsius
  mutate(air_temperature = air_temperature - 273.15) |> 
  select(datetime, all_of(variables))
```

Now, we can use the data to fit our model.

```{r}
site = "HARV"
target_variable = "gcc_90"
horiz = 35
#forecast_site <- function(site, target_variable, horiz,step) {

message(paste0("Running site: ", site))

# Get site information for elevation
#site_info <- site_data |> dplyr::filter(field_site_id == site)

# Format site data for arima model
site_target_raw <- target |>
  dplyr::select(datetime, site_id, variable, observation) |>
  dplyr::filter(variable == target_variable, 
                site_id == site) |> 
  tidyr::pivot_wider(names_from = "variable", values_from = "observation")

#  if(!target_variable%in%names(site_target_raw)||sum(!is.na(site_target_raw[target_variable]))==0){
#    message(paste0("No target observations at site ",site,". Skipping forecasts at this site."))
#    return()

# } else {

# if(theme %in% c("ticks","beetles")){
#   site_target = site_target_raw %>%
#     filter(wday(datetime,label = T)=="Mon")|>
#     complete(datetime = full_seq(datetime,step),site_id)
#   #Find the most recent Monday
#   mon = Sys.Date()-abs(1-as.numeric(strftime(Sys.Date(), "%u")))
#   h = as.numeric(floor((mon-max(site_target$datetime))/step)+horiz)
# } else {
site_target = site_target_raw |>
  complete(datetime = full_seq(datetime,1),site_id)
h = as.numeric(forecast_date-max(site_target$datetime)+horiz)
#}

site_target_past = site_target |> filter(datetime < forecast_date) |>
  right_join(met_daily,"datetime") ## merge in covariate data

# Fit arima model
if(sum(site_target[target_variable]<0,na.rm=T)>0){#If there are any negative values, don't consider transformation
  fit = forecast::auto.arima(site_target_past[target_variable])
} else {
  fit = forecast::auto.arima(site_target_past[target_variable], 
                             xreg = as.matrix(site_target_past[,c("air_temperature","relative_humidity")]),
                             #                       lambda = "auto",
                             max.d = 0, max.D = 0,max.q=0)
}

lambda = fit$lambda
boxcox = function(x,lambda){(x^lambda-1)/lambda}
inv.boxcox = function(x,lambda){
  ((x * lambda) + 1)^(1 / lambda) - 1
}
```

Now, since we want to create an ensemble forecast, we create a function that will do that

```{r}
arima.fx <- function(inputs,drivers,epsilon,
                     horiz=35,lag=1,param.names = names(fit$coef)){
  IC = inputs[,"dayof"]
  met.ens = inputs[,"noaa_ensemble_member"]
  sig.ens = inputs[,"sigma.ens"]
  param   = inputs[,param.names]
  betas = param[,which(colnames(param) %in% variables)]
  if(is.null(dim(IC))) IC = as.matrix(IC,ncol=1)
  IC.bc = IC - param[,"intercept"] #boxcox(IC,lambda)
  X = matrix(NA,nrow=nrow(IC),ncol=horiz+lag+1)
  X[,seq_len(lag)] = IC.bc
  dates = sort(unique(drivers$datetime))
  for(t in lag + (0:horiz)){
    met = drivers |>
      filter(datetime == dates[t]) |>
      ungroup() |>
      select(air_temperature,relative_humidity)
    met = met[met.ens,]
    ## gap filling
    met$air_temperature[is.na(met$air_temperature)]     = mean(met$air_temperature,na.rm = TRUE)
    met$relative_humidity[is.na(met$relative_humidity)] = mean(met$relative_humidity,na.rm = TRUE)
    #met[,nrow(met)+1] = apply(met,2,mean,na.rm=TRUE)
    
#    X[,t+1] = unlist(X[,t] + met[,1] * betas[,1] + met[,2]*betas[,2] + epsilon[sig.ens,t] )
    X[,t+1] = unlist(met[,1] * betas[,1] + met[,2]*betas[,2] + epsilon[sig.ens,t] )
    

  }
  y = X + param[,"intercept"]#inv.boxcox(X[,-lag],lambda)
  return(y)
}
```

### Machine Learning Model (Random Forest)

```{r}

site_target_past_rf <- site_target_past[!is.na(site_target_past$gcc_90), ]
site_target_past_rf <- site_target_past_rf |> select(-datetime, -site_id)

train_site <- function(sites, site_target_past, target_variable) {
  message(paste0("Running ",target_variable," at all sites"))
  
  if(!target_variable%in%names(site_target_past)){
    message(paste0("No target observations at site ",site,". Skipping forecasts at this site."))
    return()
    
  } else if(sum(!is.na(site_target_past$air_temperature)&!is.na(site_target_past[target_variable]))==0){
    message(paste0("No historical air temp data that corresponds with target observations at site ",site,". Skipping forecasts at this site."))
    return()
    
  } else {
    # Tune and fit lasso model - making use of tidymodels
    
    #Recipe for training models
    rec_base <- recipe(site_target_past)|>
      #step_rm(c("datetime", "site_id", "parameter"))|>
      update_role(everything(), new_role = "predictor")|>
      update_role({{target_variable}}, new_role = "outcome")|>
      #step_dummy(site_id)|> #Random forest handles categorical predictor without need to convert to dummy 
      step_normalize(all_numeric(), -all_outcomes())
    
    ## Set up tuning and fitting engine
    tune_randfor <- rand_forest(
      mtry = tune(),
      trees = 500,
      min_n = tune()) |>
      set_mode("regression") %>%
      set_engine("ranger", importance = "impurity", keep.inbag = TRUE, write.forest = TRUE) 
    
    #k-fold cross-validation
    randfor_resamp <- vfold_cv(site_target_past, v = 10, repeats = 5)# define k-fold cross validation procedure 
    ## Assemble workflow and tune
    wf <- workflow() %>%
      add_recipe(rec_base)
    
    #Tune models
    #If running in parallel  
    library(doParallel)
    cl <- makePSOCKcluster(Nproc) #SET 
    registerDoParallel(cl) 
    randfor_grid <- 
      tune_grid(
        wf %>% add_model(tune_randfor),
        resamples = randfor_resamp,
        grid = 20 
      )
    
    ## Select best model via RMSE
    best_mod<-randfor_grid|>
      select_best("rmse")
    
    #select model with best tuning parameter by RMSE, cross-validation approach
    final_mod <- finalize_workflow(
      wf %>% add_model(tune_randfor),
      best_mod
    )
    
    final_fit <- fit(final_mod, site_target_past)
    
    vip <- final_fit|>extract_fit_parsnip()|>vip()|>pluck("data")|>
      pivot_wider(names_from = "Variable", values_from = "Importance", names_prefix = "importance_")
    
    final_preds <- predict(final_fit, site_target_past)|>
      bind_cols(site_target_past)
    
    final_rmse<-rmse(final_preds, estimate = .pred, truth = {{target_variable}})
    #try to extract fit and write to tibble variable importance as columns bind_cols
    
    #save model fit in minimal form
    res_bundle <-
      final_fit %>%            
      bundle()
    
    return(res_bundle)
    
  }
}

rf_gcc_model <- train_site(site = site, target_variable = target_variable, site_target_past = site_target_past_rf)

```

## Methods of Uncertianty Analysis

### Sources of Uncertainty

The inputs to our models are not perfect values. They contain uncertainty, which can be binned to a number of sources. The main sources we consider are uncertainty from parameters, drivers, initial conditions, random effects, observation error, and process error.

To propagate and analyze these uncertainties, we first need to create design matrices that sample from each of these inputs. In the code below we create two such samples, x1 and x2. We will use the first immediately to generate an ensemble forecast, and then will use both together later in a Sobol analysis, which will help us partition out the different sources of uncertainty.


```{r}
library('sensitivity')
library('mvtnorm')

#Number of design samples
n <- 100000
ne <- 100000
lag = 1

#Design X1 and X2
par1 <- rmvnorm(n, fit$coef, fit$var.coef)
par2 <- rmvnorm(n, fit$coef, fit$var.coef)

## sample drivers by sampling their ensemble member numbers
d1map <- as.matrix(sample(31, n, replace = TRUE))
colnames(d1map) <- c('noaa_ensemble_member')
d2map <- as.matrix(sample(31, n, replace = TRUE))
colnames(d2map) <- c('noaa_ensemble_member')

## sample initial conditions
ICuncert <- forecast(fit, h=0, level=0.68, xreg = as.matrix(site_target_past[, c('air_temperature', 'relative_humidity')]))
IC.mean = head(ICuncert$mean,n=lag)
IC.sd   = (head(ICuncert$upper,n=lag) - head(ICuncert$lower,n=lag))/2 
ic1 <- rmvnorm(n,as.vector(IC.mean),diag(as.vector(IC.sd)^2,nrow = lag))
colnames(ic1) <- tail(c('lag2', 'lag1', 'dayof'),n = lag)
ic2 <- rmvnorm(n, as.vector(IC.mean), diag(as.vector(IC.sd)^2, nrow = lag))
colnames(ic2) <- tail(c('lag2', 'lag1', 'dayof'), n =lag)

## sample process error by sampling row numbers from a matrix of Normal random numbers. This approach is specifically to support the Sobol analyses, otherwise one could simply have used `epsilon` directly
epsilon = rmvnorm(ne,rep(0,horiz+lag),diag(rep(fit$sigma2,horiz+lag)))  ## here our SIGMA matrix is diagonal, but using rmvnorm to show how this could be generalized to autocorrelated errors
sigma.ens1 = data.frame(sigma.ens = 1:ne)
sigma.ens2 = data.frame(sigma.ens = sample(1:ne,ne,replace=TRUE))

## combine samples together into overall input matrix
x1 <- data.frame(cbind(par1, d1map, ic1, sigma.ens1))
x2 <- data.frame(cbind(par2, d2map, ic2, sigma.ens2))
```



### Validation Analysis

When looking at model uncertainty there is a conceptual distinction between the _predictive uncertainty_ and the _validation uncertainty_. Predictive uncertainty is the result of uncertainty propagation through a model and integrates some or all of the sources of uncertainty discussed above. Validation uncertainty looks at the realized performance of a forecast through some choice of skill score (e.g. RMSE, predictive quantiles, CPRS). Ideally, in a well calibrated forecast the predictive and valdation uncertainties should converge, but this doesn't always happen in practice (e.g. predictive forecast uncertainties may be over- or under-confident) so it is always good to check model performance. In this section we look at the validation uncertainty of the ARIMA forecast.

```{r}
y1 = arima.fx(x1,drivers,epsilon,horiz=35)

#### Validation uncertainty #### 
yobs = site_target |>
  filter(between(datetime,forecast_date,forecast_date+lubridate::days(horiz+1)))

# model time series with CI and data
ybar = apply(y1,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)

plot(yobs$datetime,ybar[2,],lwd=3,
     ylim=range(ybar,na.rm = TRUE),type='l',
     ylab="gcc90",xlab="date")
lines(yobs$datetime,ybar[1,],lty=2)
lines(yobs$datetime,ybar[3,],lty=2)
points(yobs$datetime,yobs$gcc_90,pch="+",cex=2,col=2)
```
Figure V1: Forecast median and 95% CI (black) compared to observations.

```{r}
# raw error vs lead time
resid = ybar[2,] - yobs$gcc_90
plot(yobs$datetime,resid,lwd=3,type='l',
     ylab="model - data",xlab="date")

```
Figure V2: Model residual error versus lead time. If we'd run multiple forecasts we would want to average these errors (e.g., RMSE), with lead time rather than date on the x-axis. Unlike traditional models, where we hope that errors are independent and there are no patterns in the residuals, with a forecast we expect the mean validation error to increase with distance / lead time and it is very common for forecast errors to be autocorrelated.


```{r}
# quantile error vs lead time
yquant = rep(NA,ncol(y1))
for(t in seq_len(ncol(y1))){
  yquant[t] = findInterval(yobs$gcc_90[t],vec = sort(y1[,t]))/nrow(y1)
}
plot(yobs$datetime,yquant,type = 'l',ylab="Predictive Quantile",xlab="date")
hist(yquant) ## should be uniform!
```
Figure V3: Model error expressed in terms of predictive quantiles. In a well-calibrated model the histogram of model quantiles should be flat (e.g., when we predict something will occur 10% of the time it actually occurs 10% of the time). While a model that had the predictive quantiles unimodally clumped around 0.5 (50%) would clearly be accurate, this would suggest that the predictive uncertainties were under-confident (i.e., that the predictive distribution is wider [more pessimistic] than it needs to be). At the other extreme, a U-shaped histogram would suggest that the model is over-confident -- the predictive distribution is too narrow relative to the realized model performance.


```{r}
# crps vs lead time
crps = rep(NA,ncol(y1))
for(t in seq_len(ncol(y1)-1)){
  crps[t] = scoringRules::crps_sample(yobs$gcc_90[t],dat = y1[,t])
}
plot(yobs$datetime,crps,type = 'l',xlab="date")
```
Figure V4: Continous Ranked Probability Score (CRPS) versus lead time. Most traditional model performance metrics (RMSE, R^2, correlation, MAE, etc.) are calculated around the _mean_ of the forecast and ignore the distribution of predictions. By contrast, CRPS is a composite score that uses _all_ of the ensemble members in our forecast, accounting for both model accuracy / bias (via a calculation of mean absolute error across ensemble members) and model precision (via a term that accounts for ensemble spread). For a given model error, the CRPS score penalizes for forecasts that are either over- or under-confident.


### More complex approaches to assessing model validation uncertainty:

When assessing the performance of a forecast that is running iteratively over time and potentially over many locations, it can be difficult to understand the patterns of forecast error and uncertainty. We recommend using flexible models (e.g. ML, GAMs) to analyze the patterns in forecast uncertainty as a function of forecast distance (e.g., lead time, spatial distance, environmental distance, phylogentic distance) and other relevant factors that may explain the heterogenity in forecast performance (e.g., environmental drivers). For example, the first paper on the EFI NEON Phenology forecast challenge used GAMs to look at forecast error as a smooth function of lead time and day of year, with site and model as additional explainatory variables

```{r}
horizon = 0:36 ## when averaging over multiple forecasts, will come from dataframe
crps.fit = data.frame(crps = crps, horizon = horizon) 
fit2 <- mgcv::gam(crps ~  s(horizon),
              data = crps.fit,
            method="REML")
summary(fit2)

hnew <- 0:35
crps_smooth <- predict(fit2)
plot(hnew,crps_smooth,xlab="Lead Time",ylab="predicted CRPS",
     type='l',lwd=1,lty=2)
```

### One-at-a-time Analysis

The one-at-a-time analysis works by incorporating the uncertainty of  

### Sobol Analysis

The Sobol analysis uses the Sobol Method, where two design matrices are created, and the program "swaps" values from design 1 with design 2 and calculates how these changes affect changes in the model ensemble output. It uses this to calculate how much the variance in the input contributes to the variance, or uncertainty, in the output. Sobol analyses result in set of Sobol Indices. The important Sobol Indices are the Sobol First Order indices and the Sobol Total Order indices. First order indices tell us about the direct contribution of the uncertainty in each input, while total order indices tell us about the total contribution with interactions that each input has.

This analysis makes use of the two different ensemble of inputs, x1 and x2, we defined earlier. With those in place we can run our Sobol analysis, which gives us estimates of our Sobol first and total order indices.

```{r}
#Sobol function call
sobolOuts <- soboljansen(model = arima.fx, X1 = x1, X2 = x2, drivers = noaa_future_daily, horiz = 35, epsilon = epsilon, conf = 0.95)
```

Now, we clean up our data, since Sobol indices cannot be below 0, but our estimates can end up sightly below zero.

```{r}
#Clean up negatives (set to 0)
sobolOuts$S[sobolOuts$S < 0] <- 0
sobolOuts$T[sobolOuts$T < 0] <- 0
```

Note that with this function, we get a first and total order index at every time point we predict for, in this case 36 days of predictions. Below is a plot of the first order indices.

```{r}
plot(sobolOuts$S[1,], type = 'l', ylim = c(0,1.15), lwd = 3, main = 'Sobol First Order Indices',
     ylab = 'Index', xlab = 'Day')
for(i in 2:nrow(sobolOuts$S)){lines(sobolOuts$S[i,], ylim = c(0,1.15), col = i, lwd = 3)}
legend('topright', legend = rownames(sobolOuts$S), col = 1:7, lty=1, lwd=3, cex = 0.4)
```

As we can see, for most of the prediction, the Sobol first order indices show that sigma.ens, representing process error is the dominant source of uncertainty. The second largest source for much of the time series is the uncertainty in our driver, NOAA weather forecasts. Parameter uncertainty contributes only a small amount to overall uncertainty. Finally, initial condition uncertainty is dominant during our first day of prediction, but quickly decays to being a very minor contributor to overall uncertainty after that.

Now, let's look at the plot of total order indices, which will include the effect of interactions.

```{r}
plot(sobolOuts$T[1,], type = 'l', ylim = c(0,1.15), lwd = 3, main = 'Sobol Total Order Indices',
     ylab = 'Index', xlab = 'Day')
for(i in 2:nrow(sobolOuts$T)){lines(sobolOuts$S[i,], ylim = c(0,1.15), col = i, lwd = 3)}
legend('topright', legend = rownames(sobolOuts$T), col = 1:7, lty=1, lwd=3, cex = 0.4)
```

Again, for most of the prediction, the Sobol total order indices show that sigma.ens, representing process error is the dominant source of uncertainty. The second largest source for much of the time series is the uncertainty in our driver, NOAA weather forecasts. Parameter uncertainty in the parameters accounting for the effect of air temperature and relative humidity are very small, while the intercept parameter contributes a small but not negligible amount. Since this behavior appears in the total order index but the first order index is much smaller, we can conclude that most of the effect of the intercept parameter uncertainty on overall uncertainty is due to interactions, not a direct effect. Finally, initial condition uncertainty is dominant during our first day of prediction, but quickly decays to being a very minor contributor to overall uncertainty after that.



### Sobol Analysis Machine Learning Model

```{r}

rec <- prep(extract_preprocessor(rf_gcc_model |> unbundle()))
N_ens <- 1000
noaa_futures <- noaa_future_daily |>
  filter(datetime >= noaa_date & datetime <= as.Date("2023-06-07")) |> 
  group_by(datetime) |> 
  group_split()

sobol_time <- tibble(
  forecast_date = seq.Date(forecast_date, as.Date("2023-06-07"), by = "day"),
  ensemble = (noaa_futures))

library(future)
library(furrr)
library(future.callr)
plan(callr, workers = Nproc)


source("R/predict_rf_ensemble.R")
model_test <- function(X, mod = rf_gcc_model$object$fit$fit$object$fit) {
  y <- predict_tidy_ensemble(mod, X)
  #  y <- y$.pred
  y
}

sobol_time_fit <- sobol_time |> 
  mutate(sobol = future_map(noaa_futures, function(x) {
    noaa_future2 <- x |> 
      select(-datetime) |> 
      sample_n(N_ens, replace = TRUE)
    
    noaa_future2 <- bake(rec, noaa_future2) |> 
      mutate(parameter_seed = sample.int(n()))
    
    noaa_future_a <- noaa_future2 |>
      filter(parameter_seed %in% 1:round(N_ens/2))
    noaa_future_b <- noaa_future2 |>
      filter(parameter_seed %in% round(N_ens/2 + 1):N_ens) 
    
    for(i in seq_len(ncol(noaa_future_b))) {
      noaa_future_b[[i]] <- sample(noaa_future_b[[i]], length(noaa_future_b[[i]]), replace = FALSE)
    }
    
    s2 <- soboljansen(model = model_test, X1 = noaa_future_a, X2 = noaa_future_b, nboot = 0, conf = 0.95)
    
    return(s2)
    
  }))


sobol_time_fit_data <- sobol_time_fit |> 
  mutate(T = map(sobol, function(x) {
    x$T |> rownames_to_column(var = "variable")
  })) |> 
  select(forecast_date, T) |> 
  unnest(T) |> 
  mutate(variable = fct_recode(variable, tree_uncertainty = "parameter_seed")) |> 
  mutate(variable = fct_relevel(variable, "tree_uncertainty", after = Inf))

ggplot(sobol_time_fit_data, aes(x = forecast_date, fill = variable, y = original)) +
  geom_area()

```