feat: time series validation - v3.9.0 (#545)

- **Feat**: new time series validation via time-series train/val/test dynamic splits and Adjusted R2 and NRMSE metrics reported for each group feature. We are adding an additional `train_size` hyperparameter to pick the size of the training size, which by default will iterate in the range of 0.5-0.8. Given it's a hyperparameter, you can change the range or fix the value manually. Turn on/off this feature using the `ts_validation` new parameter on `robyn_run()`; default is set to `FALSE` for now. This is an important step for the forecasting coming function.
- **Feat**: new `ts_validation()` function to plot time-series validation and convergence results. Generated and exported by default when `ts_validation = TRUE`, and when `export = TRUE`, creating `ts_validation_plot.png` file.
- **Fix**: updated Adjusted R2 calculation for time-series validation using same denominator.
- **Fix**: results are not sorted by lowest errors now to keep iteration results actual order.
- **Feat**: added prophet monthly component to enrich decomposition results.
- **Fix**: correct solID (not 1_1_1) for fixed hyperparameters recreated models.
- **Recode**: reduced the size of `xDecompVec` on `OutputCollect` by keeping pareto-front models only.
- **Docs**: adapted standard inputs for window start and end to include more data (3 years by default)

Co-authored-by: Gufeng Zhou <gufengzhou@fb.com>

Author: laresbernardo

R/DESCRIPTION

@@ -1,7 +1,7 @@
 Package: Robyn
 Type: Package
 Title: Semi-Automated Marketing Mix Modeling (MMM) from Meta Marketing Science 
-Version: 3.8.2
+Version: 3.9.0
 Authors@R: c(
     person("Gufeng", "Zhou", , "gufeng@meta.com", c("aut")),
     person("Leonel", "Sentana", , "leonelsentana@meta.com", c("aut")),

R/NAMESPACE

@@ -38,6 +38,7 @@ export(robyn_train)
 export(robyn_update)
 export(robyn_write)
 export(saturation_hill)
+export(ts_validation)
 import(ggplot2)
 importFrom(doParallel,registerDoParallel)
 importFrom(doParallel,stopImplicitCluster)
@@ -108,6 +109,7 @@ importFrom(patchwork,plot_annotation)
 importFrom(patchwork,plot_layout)
 importFrom(patchwork,wrap_plots)
 importFrom(prophet,add_regressor)
+importFrom(prophet,add_seasonality)
 importFrom(prophet,fit.prophet)
 importFrom(prophet,prophet)
 importFrom(rPref,low)

R/R/auxiliary.R

@@ -4,16 +4,22 @@
 # LICENSE file in the root directory of this source tree.
 
 # Calculate R-squared
-get_rsq <- function(true, predicted, p = NULL, df.int = NULL) {
+get_rsq <- function(true, predicted, p = NULL, df.int = NULL, n_train = NULL) {
   sse <- sum((predicted - true)^2)
   sst <- sum((true - mean(true))^2)
-  rsq <- 1 - sse / sst
+  rsq <- 1 - sse / sst # rsq interpreted as variance explained
+  rsq_out <- rsq
   if (!is.null(p) && !is.null(df.int)) {
-    n <- length(true)
+    if (!is.null(n_train)) {
+      n <- n_train # for oos dataset, use n from train set for adj. rsq
+    } else {
+      n <- length(true)
+    }
     rdf <- n - p - 1
-    rsq <- 1 - (1 - rsq) * ((n - df.int) / rdf)
+    rsq_adj <- 1 - (1 - rsq) * ((n - df.int) / rdf)
+    rsq_out <- rsq_adj
   }
-  return(rsq)
+  return(rsq_out)
 }
 
 # Robyn colors

R/R/checks.R

@@ -6,6 +6,8 @@
 ############# Auxiliary non-exported functions #############
 
 opts_pnd <- c("positive", "negative", "default")
+other_hyps <- c("lambda", "train_size")
+hyps_name <- c("thetas", "shapes", "scales", "alphas", "gammas")
 
 check_nas <- function(df) {
   name <- deparse(substitute(df))
@@ -151,7 +153,7 @@ check_prophet <- function(dt_holidays, prophet_country, prophet_vars, prophet_si
   if (is.null(dt_holidays) || is.null(prophet_vars)) {
     return(invisible(NULL))
   } else {
-    opts <- c("trend", "season", "weekday", "holiday")
+    opts <- c("trend", "season", "monthly", "weekday", "holiday")
     if (!all(prophet_vars %in% opts)) {
       stop("Allowed values for 'prophet_vars' are: ", paste(opts, collapse = ", "))
     }
@@ -426,7 +428,10 @@ check_hyperparameters <- function(hyperparameters = NULL, adstock = NULL,
       "robyn_inputs(InputCollect = InputCollect, hyperparameters = ...)"
     ))
   } else {
-    hyperparameters <- hyperparameters[which(!names(hyperparameters) %in% "lambda")]
+    # Non-adstock hyperparameters check
+    check_train_size(hyperparameters)
+    # Adstock hyperparameters check
+    hyperparameters <- hyperparameters[which(!names(hyperparameters) %in% other_hyps)]
     hyperparameters_ordered <- hyperparameters[order(names(hyperparameters))]
     get_hyp_names <- names(hyperparameters_ordered)
     ref_hyp_name_spend <- hyper_names(adstock, all_media = paid_media_spends)
@@ -470,6 +475,17 @@ check_hyperparameters <- function(hyperparameters = NULL, adstock = NULL,
   }
 }
 
+check_train_size <- function(hyps) {
+  if ("train_size" %in% names(hyps)) {
+    if (!length(hyps$train_size) %in% 1:2) {
+      stop("Hyperparameter 'train_size' must be length 1 (fixed) or 2 (range)")
+    }
+    if (any(hyps$train_size <= 0.1) || any(hyps$train_size > 1)) {
+      stop("Hyperparameter 'train_size' values must be defined between 0.1 and 1")
+    }
+  }
+}
+
 check_hyper_limits <- function(hyperparameters, hyper) {
   hyper_which <- which(endsWith(names(hyperparameters), hyper))
   if (length(hyper_which) == 0) {
@@ -692,8 +708,8 @@ check_hyper_fixed <- function(InputCollect, dt_hyper_fixed, add_penalty_factor)
   hyper_fixed <- !is.null(dt_hyper_fixed)
   # Adstock hyper-parameters
   hypParamSamName <- hyper_names(adstock = InputCollect$adstock, all_media = InputCollect$all_media)
-  # Add lambda hyper-parameter
-  hypParamSamName <- c(hypParamSamName, "lambda")
+  # Add lambda and other hyper-parameters manually
+  hypParamSamName <- c(hypParamSamName, other_hyps)
   # Add penalty factor hyper-parameters names
   if (add_penalty_factor) {
     for_penalty <- names(select(InputCollect$dt_mod, -.data$ds, -.data$dep_var))
@@ -734,10 +750,14 @@ check_init_msg <- function(InputCollect, cores) {
     "Using", InputCollect$adstock, "adstocking with",
     length(InputCollect$hyper_updated), "hyperparameters", det
   )
-  if (check_parallel()) {
-    message(paste(base, "on", cores, "cores"))
+  if (cores == 1) {
+    message(paste(base, "with no parallel computation"))
   } else {
-    message(paste(base, "on 1 core (Windows fallback)"))
+    if (check_parallel()) {
+      message(paste(base, "on", cores, "cores"))
+    } else {
+      message(paste(base, "on 1 core (Windows fallback)"))
+    }
   }
 }
 

R/R/clusters.R

@@ -186,14 +186,6 @@ confidence_calcs <- function(xDecompAgg, cls, all_paid, dep_var_type, k, boot_n
         ci_low <- ifelse(boot_res$ci[1] < 0, 0, boot_res$ci[1])
         ci_up <- boot_res$ci[2]
 
-        ## Experiment with gamma distribution fitting
-        # mod_gamma <- nloptr(x0 = c(1, 1), eval_f = gamma_mle, lb = c(0,0),
-        #                     x = unlist(df_chn$roi_total),
-        #                     opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5))
-        # gamma_params <- mod_gamma$solution
-        # g_low = qgamma(0.025, shape=gamma_params[[1]], scale= gamma_params[[2]])
-        # g_up = qgamma(0.975, shape=gamma_params[[1]], scale= gamma_params[[2]])
-
         # Collect loop results
         chn_collect[[i]] <- df_chn %>%
           mutate(
@@ -283,13 +275,6 @@ errors_scores <- function(df, balance = rep(1, 3)) {
   return(scores)
 }
 
-# gamma_mle <- function(params, x) {
-#   gamma_shape <- params[[1]]
-#   gamma_scale <- params[[2]]
-#   # Negative log-likelihood
-#   return(-sum(dgamma(x, shape = gamma_shape, scale = gamma_scale, log = TRUE)))
-# }
-
 # ROIs data.frame for clustering (from xDecompAgg or pareto_aggregated.csv)
 .prepare_df <- function(x, all_media, dep_var_type) {
   check_opts(all_media, unique(x$rn))

R/R/convergence.R

@@ -42,13 +42,8 @@ robyn_converge <- function(OutputModels, n_cuts = 20, sd_qtref = 3, med_lowb = 2
   stopifnot(n_cuts > min(c(sd_qtref, med_lowb)) + 1)
 
   # Gather all trials
-  get_lists <- as.logical(grepl("trial", names(OutputModels)) * unlist(lapply(OutputModels, is.list)))
-  OutModels <- OutputModels[get_lists]
-  for (i in seq_along(OutModels)) {
-    if (i == 1) df <- data.frame()
-    temp <- OutModels[[i]]$resultCollect$resultHypParam %>% mutate(trial = i)
-    df <- rbind(df, temp)
-  }
+  get_trials <- which(names(OutputModels) %in% paste0("trial", seq(OutputModels$trials)))
+  df <- bind_rows(lapply(OutputModels[get_trials], function(x) x$resultCollect$resultHypParam))
   calibrated <- isTRUE(sum(df$mape) > 0)
 
   # Calculate deciles
@@ -199,3 +194,37 @@ robyn_converge <- function(OutputModels, n_cuts = 20, sd_qtref = 3, med_lowb = 2
 
   return(invisible(cvg_out))
 }
+
+test_cvg <- function() {
+  # Experiment with gamma distribution fitting
+  gamma_mle <- function(params, x) {
+    gamma_shape <- params[[1]]
+    gamma_scale <- params[[2]]
+    # Negative log-likelihood
+    return(-sum(dgamma(x, shape = gamma_shape, scale = gamma_scale, log = TRUE)))
+  }
+  f_geo <- function(a, r, n) {
+    for (i in 2:n) a[i] <- a[i - 1] * r
+    return(a)
+  }
+  seq_nrmse <- f_geo(5, 0.7, 100)
+  df_nrmse <- data.frame(x = 1:100, y = seq_nrmse, type = "true")
+  mod_gamma <- nloptr(
+    x0 = c(1, 1), eval_f = gamma_mle, lb = c(0, 0),
+    x = seq_nrmse,
+    opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5)
+  )
+  gamma_params <- mod_gamma$solution
+  seq_nrmse_gam <- 1 / dgamma(seq_nrmse, shape = gamma_params[[1]], scale = gamma_params[[2]])
+  seq_nrmse_gam <- seq_nrmse_gam / (max(seq_nrmse_gam) - min(seq_nrmse_gam))
+  seq_nrmse_gam <- max(seq_nrmse) * seq_nrmse_gam
+  range(seq_nrmse_gam)
+  range(seq_nrmse)
+  df_nrmse_gam <- data.frame(x = 1:100, y = seq_nrmse_gam, type = "pred")
+  df_nrmse <- bind_rows(df_nrmse, df_nrmse_gam)
+  p <- ggplot(df_nrmse, aes(.data$x, .data$y, color = .data$type)) +
+    geom_line()
+  return(p)
+  # g_low = qgamma(0.025, shape=gamma_params[[1]], scale= gamma_params[[2]])
+  # g_up = qgamma(0.975, shape=gamma_params[[1]], scale= gamma_params[[2]])
+}

R/R/exports.R

@@ -35,7 +35,6 @@ robyn_save <- function(InputCollect,
     )
 
   # Nice and tidy table format for hyper-parameters
-  hyps_name <- c("thetas", "shapes", "scales", "alphas", "gammas")
   regex <- paste(paste0("_", hyps_name), collapse = "|")
   hyps <- filter(OutputCollect$resultHypParam, .data$solID == select_model) %>%
     select(contains(hyps_name)) %>%
@@ -129,7 +128,13 @@ print.robyn_save <- function(x, ...) {
   print(glued(
     "\n\nModel's Performance and Errors:\n    {errors}",
     errors = paste(
-      "R2 (train):", signif(x$errors$rsq_train, 4),
+      sprintf(
+        "R2 (%s): %s)",
+        ifelse(!isTRUE(x$ExportedModel$ts_validation), "train", "test"),
+        ifelse(!isTRUE(x$ExportedModel$ts_validation),
+          signif(x$errors$rsq_train, 4), signif(x$errors$rsq_test, 4)
+        )
+      ),
       "| NRMSE =", signif(x$errors$nrmse, 4),
       "| DECOMP.RSSD =", signif(x$errors$decomp.rssd, 4),
       "| MAPE =", signif(x$errors$mape, 4)
@@ -144,7 +149,7 @@ print.robyn_save <- function(x, ...) {
     replace(., . == "NA", "-") %>% as.data.frame())
 
   print(glued(
-    "\n\nHyper-parameters for channel transformations:\n    Adstock: {x$adstock}"
+    "\n\nHyper-parameters:\n    Adstock: {x$adstock}"
   ))
 
   print(as.data.frame(x$hyper_df))

R/R/imports.R

@@ -35,7 +35,7 @@
 #' @importFrom nloptr nloptr
 #' @importFrom parallel detectCores
 #' @importFrom patchwork guide_area plot_layout plot_annotation wrap_plots
-#' @importFrom prophet add_regressor fit.prophet prophet
+#' @importFrom prophet add_regressor add_seasonality fit.prophet prophet
 #' @importFrom reticulate tuple use_condaenv import conda_create conda_install py_module_available
 #' virtualenv_create py_install use_virtualenv
 #' @importFrom rPref low psel

R/R/inputs.R

@@ -418,7 +418,7 @@ Adstock: {x$adstock}
     },
     hyps = if (!is.null(x$hyperparameters)) {
       glued(
-        "Hyper-parameters for channel transformations:\n{flatten_hyps(x$hyperparameters)}"
+        "Hyper-parameters ranges:\n{flatten_hyps(x$hyperparameters)}"
       )
     } else {
       paste("Hyper-parameters:", "\033[0;31mNot set yet\033[0m")
@@ -514,11 +514,14 @@ Adstock: {x$adstock}
 #' @export
 hyper_names <- function(adstock, all_media) {
   adstock <- check_adstock(adstock)
-  global_name <- c("thetas", "shapes", "scales", "alphas", "gammas", "lambdas")
   if (adstock == "geometric") {
-    local_name <- sort(apply(expand.grid(all_media, global_name[grepl("thetas|alphas|gammas", global_name)]), 1, paste, collapse = "_"))
+    local_name <- sort(apply(expand.grid(all_media, hyps_name[
+      grepl("thetas|alphas|gammas", hyps_name)
+    ]), 1, paste, collapse = "_"))
   } else if (adstock %in% c("weibull_cdf", "weibull_pdf")) {
-    local_name <- sort(apply(expand.grid(all_media, global_name[grepl("shapes|scales|alphas|gammas", global_name)]), 1, paste, collapse = "_"))
+    local_name <- sort(apply(expand.grid(all_media, hyps_name[
+      grepl("shapes|scales|alphas|gammas", hyps_name)
+    ]), 1, paste, collapse = "_"))
   }
   return(local_name)
 }
@@ -769,6 +772,7 @@ prophet_decomp <- function(dt_transform, dt_holidays,
   use_trend <- "trend" %in% prophet_vars
   use_holiday <- "holiday" %in% prophet_vars
   use_season <- "season" %in% prophet_vars | "yearly.seasonality" %in% prophet_vars
+  use_monthly <- "monthly" %in% prophet_vars
   use_weekday <- "weekday" %in% prophet_vars | "weekly.seasonality" %in% prophet_vars
 
   dt_regressors <- bind_cols(recurrence, select(
@@ -791,6 +795,12 @@ prophet_decomp <- function(dt_transform, dt_holidays,
   custom_params$yearly.seasonality <- custom_params$weekly.seasonality <- NULL
   prophet_params <- append(prophet_params, custom_params)
   modelRecurrence <- do.call(prophet, as.list(prophet_params))
+  if (use_monthly) {
+    modelRecurrence <- add_seasonality(
+      modelRecurrence,
+      name = "monthly", period = 30.5, fourier.order = 5
+    )
+  }
 
   # dt_regressors <<- dt_regressors
   # modelRecurrence <<- modelRecurrence
@@ -821,12 +831,13 @@ prophet_decomp <- function(dt_transform, dt_holidays,
       # dt_regressors <<- dt_regressors
     }
     mod <- fit.prophet(modelRecurrence, dt_regressors)
-    forecastRecurrence <- predict(mod, dt_regressors)
+    forecastRecurrence <- predict(mod, dt_regressors) # prophet::prophet_plot_components(modelRecurrence, forecastRecurrence)
   }
 
   these <- seq_along(unlist(recurrence[, 1]))
   if (use_trend) dt_transform$trend <- forecastRecurrence$trend[these]
   if (use_season) dt_transform$season <- forecastRecurrence$yearly[these]
+  if (use_monthly) dt_transform$monthly <- forecastRecurrence$monthly[these]
   if (use_weekday) dt_transform$weekday <- forecastRecurrence$weekly[these]
   if (use_holiday) dt_transform$holiday <- forecastRecurrence$holidays[these]
   return(dt_transform)