R/beetle_null_forecast.R

In eco4cast/NEONeco4cast:

##' @title Generate null forecast from the Beetle group
##' @return None
##' @param targets_dir, directory where target csv is located
##' @param forecast_dir, directory where forecast will be saved
##' @export
##'
##' @author Carl Boettiger
##' @author Anna Spiers
##'
##'
beetle_null_forecast <- function(targets_dir, forecast_dir){

  ## For illustrative purposes, restrict forecast to 2019
  forecast_year <- 2019

  ## Read in the target data.
  ## NOTE: in general a forecast may instead be made directly from the
  ## the raw data, and may include other drivers.  Using only the target
  ## variables for prediction is merely a minimal model.
  richness <- readr::read_csv(file.path(targets_dir,"richness.csv"))
  abund <- readr::read_csv(file.path(targets_dir,"abund.csv"))

  ## Forecast is just based on historic mean/sd by siteID & month
  richness_model <- richness %>%
    dplyr::filter(year < forecast_year) %>%
    dplyr::group_by(month, siteID) %>%
    dplyr::summarize(mean = mean(n, na.rm = TRUE),
              sd = sd(n, na.rm = TRUE)) %>%
    dplyr::mutate(sd = replace_na(sd, mean(sd, na.rm=TRUE))) %>%
    dplyr::mutate(year = forecast_year)


  abund_model <- abund %>%
    dplyr::filter(year < forecast_year) %>%
    dplyr::group_by(month, siteID) %>%
    dplyr::summarize(mean = mean(abund, na.rm=TRUE),
              sd = sd(abund, na.rm=TRUE))  %>%
    dplyr::mutate(sd = replace_na(sd, mean(sd, na.rm=TRUE))) %>%
    dplyr::mutate(year = forecast_year)

  ### Express forecasts in terms of replicates instead of analytic mean, sd.
  ### This allows for scoring using CRPS, and generalizes to MCMC-based forecasts

  mcmc_samples <- function(df, n_reps = 500){
    ids <- df %>%
      dplyr::mutate(id = paste(siteID, year, month, sep="-")) %>%
      dplyr::pull(id)
    purrr::map_dfr(seq_along(ids),
            function(i) data.frame(id = ids[i],
                                   rep = 1:n_reps,
                                   y = rnorm(n_reps, df$mean[[i]], df$sd[[i]])))
  }

  richness_forecast <- mcmc_samples(richness_model, n_reps)
  abund_forecast <- mcmc_samples(abund_model, n_reps)


  ## Publish the forecast products

  readr::write_csv(richness_forecast, file.path(forecast_dir,"richness_forecast.csv"))
  readr::write_csv(abund_forecast,  file.path(forecast_dir,"abund_forecast.csv"))

  #publish(c(file.path(forecast_dir,"richness_forecast.csv"), file.path(forecast_dir,"abund_forecast.csv")))




}