From 7c010aac51caf7b5dcb83780ef11ae0bf3896dd3 Mon Sep 17 00:00:00 2001
From: Ben Letham <bletham@gmail.com>
Date: Tue, 15 May 2018 10:07:38 -0700
Subject: [PATCH] Multiplicative seasonality (R)

---
 R/R/diagnostics.R                      |   3 +-
 R/R/plot.R                             |  54 +++--
 R/R/prophet.R                          | 306 ++++++++++++++++++-------
 R/inst/stan/prophet.stan               |   9 +-
 R/man/add_regressor.Rd                 |  10 +-
 R/man/add_seasonality.Rd               |  12 +-
 R/man/make_all_seasonality_features.Rd |   4 +
 R/man/make_holiday_features.Rd         |   1 +
 R/man/plot_forecast_component.Rd       |   4 +-
 R/man/predictive_samples.Rd            |   5 +-
 R/man/prophet.Rd                       |  10 +-
 R/man/prophet_plot_components.Rd       |  10 +-
 R/man/regressor_column_matrix.Rd       |  29 +++
 R/man/sample_model.Rd                  |   8 +-
 R/man/sample_posterior_predictive.Rd   |   3 +-
 R/tests/testthat/test_diagnostics.R    |   5 +-
 R/tests/testthat/test_prophet.R        | 149 +++++++++---
 R/tests/testthat/test_stan_functions.R |  12 +-
 python/fbprophet/forecaster.py         |   6 +-
 python/fbprophet/plot.py               |   2 +-
 20 files changed, 489 insertions(+), 153 deletions(-)
 create mode 100644 R/man/regressor_column_matrix.Rd

diff --git a/R/R/diagnostics.R b/R/R/diagnostics.R
index cc3ea58..5fea36f 100644
--- a/R/R/diagnostics.R
+++ b/R/R/diagnostics.R
@@ -179,6 +179,7 @@ prophet_copy <- function(m, cutoff = NULL) {
     weekly.seasonality = FALSE,
     daily.seasonality = FALSE,
     holidays = m$holidays,
+    seasonality.mode = m$seasonality.mode,
     seasonality.prior.scale = m$seasonality.prior.scale,
     changepoint.prior.scale = m$changepoint.prior.scale,
     holidays.prior.scale = m$holidays.prior.scale,
@@ -253,7 +254,7 @@ performance_metrics <- function(df, metrics = NULL, rolling_window = 0.1) {
     cols <- c(cols, metric)
   }
   df_m <- df_m[cols]
-  return(na.omit(df_m))
+  return(stats::na.omit(df_m))
 }
 
 #' Compute a rolling mean of x
diff --git a/R/R/plot.R b/R/R/plot.R
index af13905..651dfc3 100644
--- a/R/R/plot.R
+++ b/R/R/plot.R
@@ -69,8 +69,9 @@ plot.prophet <- function(x, fcst, uncertainty = TRUE, plot_cap = TRUE,
 }
 
 #' Plot the components of a prophet forecast.
-#' Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if
-#' present, and holidays if present.
+#' Prints a ggplot2 with whichever are available of: trend, holidays, weekly
+#' seasonality, yearly seasonality, and additive and multiplicative extra
+#' regressors.
 #'
 #' @param m Prophet object.
 #' @param fcst Data frame returned by predict(m, df).
@@ -94,11 +95,12 @@ prophet_plot_components <- function(
   yearly_start = 0
 ) {
   # Plot the trend
-  panels <- list(plot_forecast_component(fcst, 'trend', uncertainty, plot_cap))
+  panels <- list(
+    plot_forecast_component(m, fcst, 'trend', uncertainty, plot_cap))
   # Plot holiday components, if present.
   if (!is.null(m$holidays) && ('holidays' %in% colnames(fcst))) {
     panels[[length(panels) + 1]] <- plot_forecast_component(
-      fcst, 'holidays', uncertainty, FALSE)
+      m, fcst, 'holidays', uncertainty, FALSE)
   }
   # Plot weekly seasonality, if present
   if ("weekly" %in% colnames(fcst)) {
@@ -116,10 +118,17 @@ prophet_plot_components <- function(
     }
   }
   # Plot extra regressors
-  if ((length(m$extra_regressors) > 0)
-      & ('extra_regressors' %in% colnames(fcst))) {
-    panels[[length(panels) + 1]] <- plot_forecast_component(
-      fcst, 'extra_regressors', uncertainty, FALSE)
+  regressors <- list(additive = FALSE, multiplicative = FALSE)
+  for (name in names(m$extra_regressors)) {
+    regressors[[m$extra_regressors[[name]]$mode]] <- TRUE
+  }
+  for (mode in c('additive', 'multiplicative')) {
+    if ((regressors[[mode]]) &
+        (paste0('extra_regressors_', mode) %in% colnames(fcst))
+    ) {
+      panels[[length(panels) + 1]] <- plot_forecast_component(
+        m, fcst, paste0('extra_regressors_', mode), uncertainty, FALSE)
+    }
   }
 
   # Make the plot.
@@ -135,6 +144,7 @@ prophet_plot_components <- function(
 
 #' Plot a particular component of the forecast.
 #'
+#' @param m Prophet model
 #' @param fcst Dataframe output of `predict`.
 #' @param name String name of the component to plot (column of fcst).
 #' @param uncertainty Boolean to plot uncertainty intervals.
@@ -145,7 +155,7 @@ prophet_plot_components <- function(
 #'
 #' @export
 plot_forecast_component <- function(
-  fcst, name, uncertainty = TRUE, plot_cap = FALSE
+  m, fcst, name, uncertainty = TRUE, plot_cap = FALSE
 ) {
   gg.comp <- ggplot2::ggplot(
       fcst, ggplot2::aes_string(x = 'ds', y = name, group = 1)) +
@@ -168,6 +178,9 @@ plot_forecast_component <- function(
         fill = "#0072B2",
         na.rm = TRUE)
   }
+  if (name %in% m$component.modes$multiplicative) {
+    gg.comp <- gg.comp + ggplot2::scale_y_continuous(labels = scales::percent)
+  }
   return(gg.comp)
 }
 
@@ -220,6 +233,11 @@ plot_weekly <- function(m, uncertainty = TRUE, weekly_start = 0) {
                            fill = "#0072B2",
                            na.rm = TRUE)
   }
+  if (m$seasonalities$weekly$mode == 'multiplicative') {
+    gg.weekly <- (
+      gg.weekly + ggplot2::scale_y_continuous(labels = scales::percent)
+    )
+  }
   return(gg.weekly)
 }
 
@@ -255,6 +273,11 @@ plot_yearly <- function(m, uncertainty = TRUE, yearly_start = 0) {
                            fill = "#0072B2",
                            na.rm = TRUE)
   }
+  if (m$seasonalities$yearly$mode == 'multiplicative') {
+    gg.yearly <- (
+      gg.yearly + ggplot2::scale_y_continuous(labels = scales::percent)
+    )
+  }
   return(gg.yearly)
 }
 
@@ -299,6 +322,9 @@ plot_seasonality <- function(m, name, uncertainty = TRUE) {
       fill = "#0072B2",
       na.rm = TRUE)
   }
+  if (m$seasonalities[[name]]$mode == 'multiplicative') {
+    gg.s <- gg.s + ggplot2::scale_y_continuous(labels = scales::percent)
+  }
   return(gg.s)
 }
 
@@ -386,8 +412,7 @@ dyplot.prophet <- function(x, fcst, uncertainty=TRUE,
       dygraphs::dyAnnotation(x, text, text, attachAtBottom = TRUE)
   }
   
-  dyBase <- dyBase %>% 
-    dygraphs::dyOptions(colors = RColorBrewer::brewer.pal(3, "Set1")) %>%
+  dyBase <- dyBase %>%
     # plot actual values
     dygraphs::dySeries('y', label=actual.label) %>% 
     # plot forecast and ribbon
@@ -399,8 +424,11 @@ dyplot.prophet <- function(x, fcst, uncertainty=TRUE,
   if (!is.null(x$holidays)) {
     for (i in 1:nrow(x$holidays)) {
       # make a gray line
-      dyBase <- dyBase %>% dygraphs::dyEvent(x$holidays$ds[i],color = "rgb(200,200,200)", strokePattern = "solid")
-      dyBase <- dyBase %>% dygraphs::dyAnnotation(x$holidays$ds[i], x$holidays$holiday[i], x$holidays$holiday[i], attachAtBottom = TRUE)
+      dyBase <- dyBase %>% dygraphs::dyEvent(
+        x$holidays$ds[i],color = "rgb(200,200,200)", strokePattern = "solid")
+      dyBase <- dyBase %>% dygraphs::dyAnnotation(
+        x$holidays$ds[i], x$holidays$holiday[i], x$holidays$holiday[i],
+        attachAtBottom = TRUE)
     }
   }
   return(dyBase)
diff --git a/R/R/prophet.R b/R/R/prophet.R
index e996dd9..86aa9b9 100644
--- a/R/R/prophet.R
+++ b/R/R/prophet.R
@@ -40,6 +40,7 @@ utils::globalVariables(c(
 #'  range of days around the date to be included as holidays. lower_window=-2
 #'  will include 2 days prior to the date as holidays. Also optionally can have
 #'  a column prior_scale specifying the prior scale for each holiday.
+#' @param seasonality.mode 'additive' (default) or 'multiplicative'.
 #' @param seasonality.prior.scale Parameter modulating the strength of the
 #'  seasonality model. Larger values allow the model to fit larger seasonal
 #'  fluctuations, smaller values dampen the seasonality. Can be specified for
@@ -82,6 +83,7 @@ prophet <- function(df = NULL,
                     weekly.seasonality = 'auto',
                     daily.seasonality = 'auto',
                     holidays = NULL,
+                    seasonality.mode = 'additive',
                     seasonality.prior.scale = 10,
                     holidays.prior.scale = 10,
                     changepoint.prior.scale = 0.05,
@@ -105,6 +107,7 @@ prophet <- function(df = NULL,
     weekly.seasonality = weekly.seasonality,
     daily.seasonality = daily.seasonality,
     holidays = holidays,
+    seasonality.mode = seasonality.mode,
     seasonality.prior.scale = seasonality.prior.scale,
     changepoint.prior.scale = changepoint.prior.scale,
     holidays.prior.scale = holidays.prior.scale,
@@ -122,7 +125,9 @@ prophet <- function(df = NULL,
     stan.fit = NULL,
     params = list(),
     history = NULL,
-    history.dates = NULL
+    history.dates = NULL,
+    train.component.cols = NULL,
+    component.modes = NULL
   )
   validate_inputs(m)
   class(m) <- append("prophet", class(m))
@@ -168,6 +173,9 @@ validate_inputs <- function(m) {
       validate_column_name(m, h, check_holidays = FALSE)
     }
   }
+  if (!(m$seasonality.mode %in% c('additive', 'multiplicative'))) {
+    stop("seasonality.mode must be 'additive' or 'multiplicative'")
+  }
 }
 
 #' Validates the name of a seasonality, holiday, or regressor.
@@ -187,8 +195,9 @@ validate_column_name <- function(
     stop('Holiday name cannot contain "_delim_"')
   }
   reserved_names = c(
-    'trend', 'seasonal', 'seasonalities', 'daily', 'weekly', 'yearly',
-    'holidays', 'zeros', 'extra_regressors', 'yhat'
+    'trend', 'additive_terms', 'daily', 'weekly', 'yearly',
+    'holidays', 'zeros', 'extra_regressors_additive', 'yhat',
+    'extra_regressors_multiplicative', 'multiplicative_terms'
   )
   rn_l = paste0(reserved_names,"_lower")
   rn_u = paste0(reserved_names,"_upper")
@@ -511,6 +520,7 @@ make_seasonality_features <- function(dates, period, series.order, prefix) {
 #' @return A list with entries
 #'  holiday.features: dataframe with a column for each holiday.
 #'  prior.scales: array of prior scales for each holiday column.
+#'  holiday.names: array of names of all holidays.
 #'
 #' @importFrom dplyr "%>%"
 #' @keywords internal
@@ -568,7 +578,8 @@ make_holiday_features <- function(m, dates) {
     prior.scales <- c(prior.scales, prior.scales.list[[sn]])
   }
   return(list(holiday.features = holiday.features,
-              prior.scales = prior.scales))
+              prior.scales = prior.scales,
+              holiday.names = names(prior.scales.list)))
 }
 
 #' Add an additional regressor to be used for fitting and predicting.
@@ -579,6 +590,10 @@ make_holiday_features <- function(m, dates) {
 #' coefficient is given a prior with the specified scale parameter.
 #' Decreasing the prior scale will add additional regularization. If no
 #' prior scale is provided, holidays.prior.scale will be used.
+#' Mode can be specified as either 'additive' or 'multiplicative'. If not
+#' specified, m$seasonality.mode will be used. 'additive' means the effect of
+#' the regressor will be added to the trend, 'multiplicative' means it will
+#' multiply the trend.
 #'
 #' @param m Prophet object.
 #' @param name String name of the regressor
@@ -587,11 +602,15 @@ make_holiday_features <- function(m, dates) {
 #' @param standardize Bool, specify whether this regressor will be standardized
 #'  prior to fitting. Can be 'auto' (standardize if not binary), True, or
 #'  False.
+#' @param mode Optional, 'additive' or 'multiplicative'. Defaults to
+#'  m$seasonality.mode.
 #'
 #' @return  The prophet model with the regressor added.
 #'
 #' @export
-add_regressor <- function(m, name, prior.scale = NULL, standardize = 'auto'){
+add_regressor <- function(
+  m, name, prior.scale = NULL, standardize = 'auto', mode = NULL
+){
   if (!is.null(m$history)) {
     stop('Regressors must be added prior to model fitting.')
   }
@@ -599,14 +618,21 @@ add_regressor <- function(m, name, prior.scale = NULL, standardize = 'auto'){
   if (is.null(prior.scale)) {
     prior.scale <- m$holidays.prior.scale
   }
+  if (is.null(mode)) {
+    mode <- m$seasonality.mode
+  }
   if(prior.scale <= 0) {
     stop("Prior scale must be > 0")
   }
+  if (!(mode %in% c('additive', 'multiplicative'))) {
+    stop("mode must be 'additive' or 'multiplicative'")
+  }
   m$extra_regressors[[name]] <- list(
     prior.scale = prior.scale,
     standardize = standardize,
     mu = 0,
-    std = 1.0
+    std = 1.0,
+    mode = mode
   )
   return(m)
 }
@@ -622,17 +648,25 @@ add_regressor <- function(m, name, prior.scale = NULL, standardize = 'auto'){
 #' flexibility, decreasing will dampen it. If not provided, will use the
 #' seasonality.prior.scale provided on Prophet initialization (defaults to 10).
 #'
+#' Mode can be specified as either 'additive' or 'multiplicative'. If not
+#' specified, m$seasonality.mode will be used (defaults to 'additive').
+#' Additive means the seasonality will be added to the trend, multiplicative
+#' means it will multiply the trend.
+#'
 #' @param m Prophet object.
 #' @param name String name of the seasonality component.
 #' @param period Float number of days in one period.
 #' @param fourier.order Int number of Fourier components to use.
-#' @param prior.scale Float prior scale for this component.
+#' @param prior.scale Optional float prior scale for this component.
+#' @param mode Optional 'additive' or 'multiplicative'.
 #'
 #' @return The prophet model with the seasonality added.
 #'
 #' @importFrom dplyr "%>%"
 #' @export
-add_seasonality <- function(m, name, period, fourier.order, prior.scale = NULL) {
+add_seasonality <- function(
+  m, name, period, fourier.order, prior.scale = NULL, mode = NULL
+) {
   if (!is.null(m$history)) {
     stop("Seasonality must be added prior to model fitting.")
   }
@@ -648,10 +682,17 @@ add_seasonality <- function(m, name, period, fourier.order, prior.scale = NULL)
   if (ps <= 0) {
     stop('Prior scale must be > 0')
   }
+  if (is.null(mode)) {
+    mode <- m$seasonality.mode
+  }
+  if (!(mode %in% c('additive', 'multiplicative'))) {
+    stop("mode must be 'additive' or 'multiplicative'")
+  }
   m$seasonalities[[name]] <- list(
     period = period,
     fourier.order = fourier.order,
-    prior.scale = ps
+    prior.scale = ps,
+    mode = mode
   )
   return(m)
 }
@@ -667,11 +708,16 @@ add_seasonality <- function(m, name, period, fourier.order, prior.scale = NULL)
 #'  seasonal.features: Dataframe with regressor features,
 #'  prior.scales: Array of prior scales for each colum of the features
 #'    dataframe.
+#'  component.cols: Dataframe with indicators for which regression components
+#'    correspond to which columns.
+#'  modes: List with keys 'additive' and 'multiplicative' with arrays of
+#'    component names for each mode of seasonality.
 #'
 #' @keywords internal
 make_all_seasonality_features <- function(m, df) {
   seasonal.features <- data.frame(row.names = seq_len(nrow(df)))
   prior.scales <- c()
+  modes <- list(additive = c(), multiplicative = c())
 
   # Seasonality features
   for (name in names(m$seasonalities)) {
@@ -681,6 +727,7 @@ make_all_seasonality_features <- function(m, df) {
     seasonal.features <- cbind(seasonal.features, features)
     prior.scales <- c(prior.scales,
                       props$prior.scale * rep(1, ncol(features)))
+    modes[[props$mode]] <- c(modes[[props$mode]], name)
   }
 
   # Holiday features
@@ -688,20 +735,129 @@ make_all_seasonality_features <- function(m, df) {
     hf <- make_holiday_features(m, df$ds)
     seasonal.features <- cbind(seasonal.features, hf$holiday.features)
     prior.scales <- c(prior.scales, hf$prior.scales)
+    modes[[m$seasonality.mode]] <- c(
+      modes[[m$seasonality.mode]], hf$holiday.names)
   }
 
   # Additional regressors
   for (name in names(m$extra_regressors)) {
+    props <- m$extra_regressors[[name]]
     seasonal.features[[name]] <- df[[name]]
-    prior.scales <- c(prior.scales, m$extra_regressors[[name]]$prior.scale)
+    prior.scales <- c(prior.scales, props$prior.scale)
+    modes[[props$mode]] <- c(modes[[props$mode]], name)
   }
 
+  # Dummy to prevent empty X
   if (ncol(seasonal.features) == 0) {
     seasonal.features <- data.frame(zeros = rep(0, nrow(df)))
     prior.scales <- 1
   }
+
+  components.list <- regressor_column_matrix(m, seasonal.features, modes)
   return(list(seasonal.features = seasonal.features,
-              prior.scales = prior.scales))
+              prior.scales = prior.scales,
+              component.cols = components.list$component.cols,
+              modes = components.list$modes))
+}
+
+#' Dataframe indicating which columns of the feature matrix correspond to
+#' which seasonality/regressor components.
+#'
+#' Includes combination components, like 'additive_terms'. These combination
+#' components will be added to the 'modes' input.
+#'
+#' @param m Prophet object.
+#' @param seasonal.features Constructed seasonal features dataframe.
+#' @param modes List with keys 'additive' and 'multiplicative' with arrays of
+#'  component names for each mode of seasonality.
+#'
+#' @return List with items
+#'  component.cols: A binary indicator dataframe with columns seasonal
+#'    components and rows columns in seasonal.features. Entry is 1 if that
+#'    column is used in that component.
+#'  modes: Updated input with combination components.
+#'
+#' @keywords internal
+regressor_column_matrix <- function(m, seasonal.features, modes) {
+  components <- dplyr::data_frame(component = colnames(seasonal.features)) %>%
+    dplyr::mutate(col = seq_len(n())) %>%
+    tidyr::separate(component, c('component', 'part'), sep = "_delim_",
+                    extra = "merge", fill = "right") %>%
+    dplyr::select(col, component)
+  # Add total for holidays
+  if(!is.null(m$holidays)){
+    components <- add_group_component(
+      components, 'holidays', unique(m$holidays$holiday))
+  }
+  # Add totals for additive and multiplicative components, and regressors
+  for (mode in c('additive', 'multiplicative')) {
+    components <- add_group_component(
+      components, paste0(mode, '_terms'), modes[[mode]])
+    regressors_by_mode <- c()
+    for (name in names(m$extra_regressors)) {
+      if (m$extra_regressors[[name]]$mode == mode) {
+        regressors_by_mode <- c(regressors_by_mode, name)
+      }
+    }
+    components <- add_group_component(
+      components, paste0('extra_regressors_', mode), regressors_by_mode)
+    # Add combination components to modes
+    modes[[mode]] <- c(modes[[mode]], paste0(mode, '_terms'))
+    modes[[mode]] <- c(modes[[mode]], paste0('extra_regressors_', mode))
+  }
+  # After all of the additive/multiplicative groups have been added,
+  modes[[m$seasonality.mode]] <- c(modes[[m$seasonality.mode]], 'holidays')
+  # Convert to a binary matrix
+  component.cols <- as.data.frame.matrix(
+    table(components$col, components$component)
+  )
+  component.cols <- (
+    component.cols[order(as.numeric(row.names(component.cols))), ,
+                   drop = FALSE]
+  )
+  # Add columns for additive and multiplicative terms, if missing
+  for (name in c('additive_terms', 'multiplicative_terms')) {
+    if (!(name %in% colnames(component.cols))) {
+      component.cols[[name]] <- 0
+    }
+  }
+  # Remove the placeholder
+  components <- dplyr::filter(components, component != 'zeros')
+  # Validation
+  if (
+    max(component.cols$additive_terms
+    + component.cols$multiplicative_terms) > 1
+  ) {
+    stop('A bug occurred in seasonal components.')
+  }
+  # Compare to training, if set.
+  if (!is.null(m$train.component.cols)) {
+    component.cols <- component.cols[, colnames(m$train.component.cols)]
+    if (!all(component.cols == m$train.component.cols)) {
+      stop('A bug occurred in constructing regressors.')
+    }
+  }
+  return(list(component.cols = component.cols, modes = modes))
+}
+
+#' Adds a component with given name that contains all of the components
+#' in group.
+#'
+#' @param components Dataframe with components.
+#' @param name Name of new group component.
+#' @param group  List of components that form the group.
+#'
+#' @return Dataframe with components.
+#'
+#' @keywords internal
+add_group_component <- function(components, name, group) {
+  new_comp <- components[(components$component %in% group), ]
+  group_cols <- unique(new_comp$col)
+  if (length(group_cols) > 0) {
+    new_comp <- data.frame(col=group_cols, component=name)
+    components <- rbind(components, new_comp)
+  }
+  return(components)
 }
 
 #' Get number of Fourier components for built-in seasonalities.
@@ -764,7 +920,8 @@ set_auto_seasonalities <- function(m) {
     m$seasonalities[['yearly']] <- list(
       period = 365.25,
       fourier.order = fourier.order,
-      prior.scale = m$seasonality.prior.scale
+      prior.scale = m$seasonality.prior.scale,
+      mode = m$seasonality.mode
     )
   }
 
@@ -775,7 +932,8 @@ set_auto_seasonalities <- function(m) {
     m$seasonalities[['weekly']] <- list(
       period = 7,
       fourier.order = fourier.order,
-      prior.scale = m$seasonality.prior.scale
+      prior.scale = m$seasonality.prior.scale,
+      mode = m$seasonality.mode
     )
   }
 
@@ -786,7 +944,8 @@ set_auto_seasonalities <- function(m) {
     m$seasonalities[['daily']] <- list(
       period = 1,
       fourier.order = fourier.order,
-      prior.scale = m$seasonality.prior.scale
+      prior.scale = m$seasonality.prior.scale,
+      mode = m$seasonality.mode
     )
   }
   return(m)
@@ -891,6 +1050,9 @@ fit.prophet <- function(m, df, ...) {
   out2 <- make_all_seasonality_features(m, history)
   seasonal.features <- out2$seasonal.features
   prior.scales <- out2$prior.scales
+  component.cols <- out2$component.cols
+  m$train.component.cols <- component.cols
+  m$component.modes <- out2$modes
 
   m <- set_changepoints(m)
 
@@ -905,7 +1067,9 @@ fit.prophet <- function(m, df, ...) {
     X = as.matrix(seasonal.features),
     sigmas = array(prior.scales),
     tau = m$changepoint.prior.scale,
-    trend_indicator = as.numeric(m$growth == 'logistic')
+    trend_indicator = as.numeric(m$growth == 'logistic'),
+    s_a = array(component.cols$additive_terms),
+    s_m = array(component.cols$multiplicative_terms)
   )
 
   # Run stan
@@ -1023,7 +1187,7 @@ predict.prophet <- function(object, df = NULL, ...) {
   }
   df <- df[cols]
   df <- dplyr::bind_cols(df, seasonal.components, intervals)
-  df$yhat <- df$trend + df$seasonal
+  df$yhat <- df$trend * (1 + df$multiplicative_terms) + df$additive_terms
   return(df)
 }
 
@@ -1118,76 +1282,37 @@ predict_trend <- function(model, df) {
 #'
 #' @keywords internal
 predict_seasonal_components <- function(m, df) {
-  seasonal.features <- make_all_seasonality_features(m, df)$seasonal.features
+  out <- make_all_seasonality_features(m, df)
+  seasonal.features <- out$seasonal.features
+  component.cols <- out$component.cols
   lower.p <- (1 - m$interval.width)/2
   upper.p <- (1 + m$interval.width)/2
 
-  components <- dplyr::data_frame(component = colnames(seasonal.features)) %>%
-    dplyr::mutate(col = seq_len(n())) %>%
-    tidyr::separate(component, c('component', 'part'), sep = "_delim_",
-                    extra = "merge", fill = "right") %>%
-    dplyr::select(col, component)
-  # Add total for all regression components
-  components <- rbind(
-    components,
-    data.frame(col = seq_len(ncol(seasonal.features)), component = 'seasonal'))
-  # Add totals for seasonality, holiday, and extra regressors
-  components <- add_group_component(
-    components, 'seasonalities', names(m$seasonalities))
-  if(!is.null(m$holidays)){
-    components <- add_group_component(
-      components, 'holidays', unique(m$holidays$holiday))
+  X <- as.matrix(seasonal.features)
+  component.predictions <- data.frame(matrix(ncol = 0, nrow = nrow(X)))
+  for (component in colnames(component.cols)) {
+    beta.c <- m$params$beta * component.cols[[component]]
+
+    comp <- X %*% t(beta.c)
+    if (component %in% m$component.modes$additive) {
+      comp <- comp * m$y.scale
+    }
+    component.predictions[[component]] <- rowMeans(comp, na.rm = TRUE)
+    component.predictions[[paste0(component, '_lower')]] <- apply(
+      comp, 1, stats::quantile, lower.p, na.rm = TRUE)
+    component.predictions[[paste0(component, '_upper')]] <- apply(
+      comp, 1, stats::quantile, upper.p, na.rm = TRUE)
   }
-  components <- add_group_component(
-    components, 'extra_regressors', names(m$extra_regressors))
-  # Remove the placeholder
-  components <- dplyr::filter(components, component != 'zeros')
-
-  component.predictions <- components %>%
-    dplyr::group_by(component) %>% dplyr::do({
-      comp <- (as.matrix(seasonal.features[, .$col])
-               %*% t(m$params$beta[, .$col, drop = FALSE])) * m$y.scale
-      dplyr::data_frame(ix = seq_len(nrow(seasonal.features)),
-                        mean = rowMeans(comp, na.rm = TRUE),
-                        lower = apply(comp, 1, stats::quantile, lower.p,
-                                      na.rm = TRUE),
-                        upper = apply(comp, 1, stats::quantile, upper.p,
-                                      na.rm = TRUE))
-    }) %>%
-    tidyr::gather(stat, value, mean, lower, upper) %>%
-    dplyr::mutate(stat = ifelse(stat == 'mean', '', paste0('_', stat))) %>%
-    tidyr::unite(component, component, stat, sep="") %>%
-    tidyr::spread(component, value) %>%
-    dplyr::select(-ix)
-
   return(component.predictions)
 }
 
-#' Adds a component with given name that contains all of the components
-#' in group.
-#'
-#' @param components Dataframe with components.
-#' @param name Name of new group component.
-#' @param group  List of components that form the group.
-#'
-#' @return Dataframe with components.
-#'
-#' @keywords internal
-add_group_component <- function(components, name, group) {
-  new_comp <- components[(components$component %in% group), ]
-  if (nrow(new_comp) > 0) {
-    new_comp$component <- name
-    components <- rbind(components, new_comp)
-  }
-  return(components)
-}
-
 #' Prophet posterior predictive samples.
 #'
 #' @param m Prophet object.
 #' @param df Prediction dataframe.
 #'
-#' @return List with posterior predictive samples for each component.
+#' @return List with posterior predictive samples for the forecast yhat and
+#'  for the trend component.
 #'
 #' @keywords internal
 sample_posterior_predictive <- function(m, df) {
@@ -1196,16 +1321,24 @@ sample_posterior_predictive <- function(m, df) {
   samp.per.iter <- max(1, ceiling(m$uncertainty.samples / n.iterations))
   nsamp <- n.iterations * samp.per.iter  # The actual number of samples
 
-  seasonal.features <- make_all_seasonality_features(m, df)$seasonal.features
+  out <- make_all_seasonality_features(m, df)
+  seasonal.features <- out$seasonal.features
+  component.cols <- out$component.cols
   sim.values <- list("trend" = matrix(, nrow = nrow(df), ncol = nsamp),
-                     "seasonal" = matrix(, nrow = nrow(df), ncol = nsamp),
                      "yhat" = matrix(, nrow = nrow(df), ncol = nsamp))
 
   for (i in seq_len(n.iterations)) {
     # For each set of parameters from MCMC (or just 1 set for MAP),
     for (j in seq_len(samp.per.iter)) {
       # Do a simulation with this set of parameters,
-      sim <- sample_model(m, df, seasonal.features, i)
+      sim <- sample_model(
+        m = m,
+        df = df,
+        seasonal.features = seasonal.features,
+        iteration = i,
+        s_a = component.cols$additive_terms,
+        s_m = component.cols$multiplicative_terms
+      )
       # Store the results
       for (key in c("trend", "seasonal", "yhat")) {
         sim.values[[key]][,(i - 1) * samp.per.iter + j] <- sim[[key]]
@@ -1221,9 +1354,8 @@ sample_posterior_predictive <- function(m, df) {
 #' @param df Dataframe with dates for predictions (column ds), and capacity
 #'  (column cap) if logistic growth.
 #'
-#' @return A list with items "trend", "seasonal", and "yhat" containing
-#'  posterior predictive samples for that component. "seasonal" is the sum
-#'  of seasonalities, holidays, and added regressors.
+#' @return A list with items "trend" and "yhat" containing
+#'  posterior predictive samples for that component.
 #'
 #' @export
 predictive_samples <- function(m, df) {
@@ -1264,22 +1396,24 @@ predict_uncertainty <- function(m, df) {
 #' @param df Prediction dataframe.
 #' @param seasonal.features Data frame of seasonal features
 #' @param iteration Int sampling iteration to use parameters from.
+#' @param s_a Indicator vector for additive components
+#' @param s_m Indicator vector for multiplicative components
 #'
-#' @return List of trend, seasonality, and yhat, each a vector like df$t.
+#' @return List of trend and yhat, each a vector like df$t.
 #'
 #' @keywords internal
-sample_model <- function(m, df, seasonal.features, iteration) {
+sample_model <- function(m, df, seasonal.features, iteration, s_a, s_m) {
   trend <- sample_predictive_trend(m, df, iteration)
 
   beta <- m$params$beta[iteration,]
-  seasonal <- (as.matrix(seasonal.features) %*% beta) * m$y.scale
+  Xb_a = as.matrix(seasonal.features) %*% (beta * s_a) * m$y.scale
+  Xb_m = as.matrix(seasonal.features) %*% (beta * s_m)
 
   sigma <- m$params$sigma_obs[iteration]
   noise <- stats::rnorm(nrow(df), mean = 0, sd = sigma) * m$y.scale
 
-  return(list("yhat" = trend + seasonal + noise,
-              "trend" = trend,
-              "seasonal" = seasonal))
+  return(list("yhat" = trend * (1 + Xb_m) + Xb_a + noise,
+              "trend" = trend))
 }
 
 #' Simulate the trend using the extrapolated generative model.
diff --git a/R/inst/stan/prophet.stan b/R/inst/stan/prophet.stan
index 746b0cb..4319b6f 100644
--- a/R/inst/stan/prophet.stan
+++ b/R/inst/stan/prophet.stan
@@ -85,6 +85,8 @@ data {
   vector[K] sigmas;     // Scale on seasonality prior
   real<lower=0> tau;    // Scale on changepoints prior
   int trend_indicator;  // 0 for linear, 1 for logistic
+  vector[K] s_a;        // Indicator of additive features
+  vector[K] s_m;        // Indicator of multiplicative features
 }
 
 transformed data {
@@ -102,12 +104,17 @@ parameters {
 
 transformed parameters {
   vector[T] trend;
+  vector[T] Xb_a;
+  vector[T] Xb_m;
 
   if (trend_indicator == 0) {
     trend = linear_trend(k, m, delta, t, A, t_change);
   } else if (trend_indicator == 1) {
     trend = logistic_trend(k, m, delta, t, cap, A, t_change, S);
   }
+
+  Xb_a = X * (beta .* s_a);
+  Xb_m = X * (beta .* s_m);
 }
 
 model {
@@ -119,5 +126,5 @@ model {
   beta ~ normal(0, sigmas);
 
   // Likelihood
-  y ~ normal(trend + X * beta, sigma_obs);
+  y ~ normal(trend .* (1 + Xb_m) + Xb_a, sigma_obs);
 }
diff --git a/R/man/add_regressor.Rd b/R/man/add_regressor.Rd
index f67f0d1..39e16a1 100644
--- a/R/man/add_regressor.Rd
+++ b/R/man/add_regressor.Rd
@@ -4,7 +4,8 @@
 \alias{add_regressor}
 \title{Add an additional regressor to be used for fitting and predicting.}
 \usage{
-add_regressor(m, name, prior.scale = NULL, standardize = "auto")
+add_regressor(m, name, prior.scale = NULL, standardize = "auto",
+  mode = NULL)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -17,6 +18,9 @@ holidays.prior.scale will be used.}
 \item{standardize}{Bool, specify whether this regressor will be standardized
 prior to fitting. Can be 'auto' (standardize if not binary), True, or
 False.}
+
+\item{mode}{Optional, 'additive' or 'multiplicative'. Defaults to
+m$seasonality.mode.}
 }
 \value{
 The prophet model with the regressor added.
@@ -28,4 +32,8 @@ regressor will be standardized unless it is binary. The regression
 coefficient is given a prior with the specified scale parameter.
 Decreasing the prior scale will add additional regularization. If no
 prior scale is provided, holidays.prior.scale will be used.
+Mode can be specified as either 'additive' or 'multiplicative'. If not
+specified, m$seasonality.mode will be used. 'additive' means the effect of
+the regressor will be added to the trend, 'multiplicative' means it will
+multiply the trend.
 }
diff --git a/R/man/add_seasonality.Rd b/R/man/add_seasonality.Rd
index 6df3f19..dc6ef98 100644
--- a/R/man/add_seasonality.Rd
+++ b/R/man/add_seasonality.Rd
@@ -5,7 +5,8 @@
 \title{Add a seasonal component with specified period, number of Fourier
 components, and prior scale.}
 \usage{
-add_seasonality(m, name, period, fourier.order, prior.scale = NULL)
+add_seasonality(m, name, period, fourier.order, prior.scale = NULL,
+  mode = NULL)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -16,7 +17,9 @@ add_seasonality(m, name, period, fourier.order, prior.scale = NULL)
 
 \item{fourier.order}{Int number of Fourier components to use.}
 
-\item{prior.scale}{Float prior scale for this component.}
+\item{prior.scale}{Optional float prior scale for this component.}
+
+\item{mode}{Optional 'additive' or 'multiplicative'.}
 }
 \value{
 The prophet model with the seasonality added.
@@ -30,4 +33,9 @@ seasonalities are 10 and 3 respectively.
 Increasing prior scale will allow this seasonality component more
 flexibility, decreasing will dampen it. If not provided, will use the
 seasonality.prior.scale provided on Prophet initialization (defaults to 10).
+
+Mode can be specified as either 'additive' or 'multiplicative'. If not
+specified, m$seasonality.mode will be used (defaults to 'additive').
+Additive means the seasonality will be added to the trend, multiplicative
+means it will multiply the trend.
 }
diff --git a/R/man/make_all_seasonality_features.Rd b/R/man/make_all_seasonality_features.Rd
index f077ea6..866a867 100644
--- a/R/man/make_all_seasonality_features.Rd
+++ b/R/man/make_all_seasonality_features.Rd
@@ -18,6 +18,10 @@ List with items
  seasonal.features: Dataframe with regressor features,
  prior.scales: Array of prior scales for each colum of the features
    dataframe.
+ component.cols: Dataframe with indicators for which regression components
+   correspond to which columns.
+ modes: List with keys 'additive' and 'multiplicative' with arrays of
+   component names for each mode of seasonality.
 }
 \description{
 Dataframe with seasonality features.
diff --git a/R/man/make_holiday_features.Rd b/R/man/make_holiday_features.Rd
index c502972..0d8d968 100644
--- a/R/man/make_holiday_features.Rd
+++ b/R/man/make_holiday_features.Rd
@@ -15,6 +15,7 @@ make_holiday_features(m, dates)
 A list with entries
  holiday.features: dataframe with a column for each holiday.
  prior.scales: array of prior scales for each holiday column.
+ holiday.names: array of names of all holidays.
 }
 \description{
 Construct a matrix of holiday features.
diff --git a/R/man/plot_forecast_component.Rd b/R/man/plot_forecast_component.Rd
index 06c1fb1..8284d64 100644
--- a/R/man/plot_forecast_component.Rd
+++ b/R/man/plot_forecast_component.Rd
@@ -4,9 +4,11 @@
 \alias{plot_forecast_component}
 \title{Plot a particular component of the forecast.}
 \usage{
-plot_forecast_component(fcst, name, uncertainty = TRUE, plot_cap = FALSE)
+plot_forecast_component(m, fcst, name, uncertainty = TRUE, plot_cap = FALSE)
 }
 \arguments{
+\item{m}{Prophet model}
+
 \item{fcst}{Dataframe output of `predict`.}
 
 \item{name}{String name of the component to plot (column of fcst).}
diff --git a/R/man/predictive_samples.Rd b/R/man/predictive_samples.Rd
index c6ff9de..7bce31b 100644
--- a/R/man/predictive_samples.Rd
+++ b/R/man/predictive_samples.Rd
@@ -13,9 +13,8 @@ predictive_samples(m, df)
 (column cap) if logistic growth.}
 }
 \value{
-A list with items "trend", "seasonal", and "yhat" containing
- posterior predictive samples for that component. "seasonal" is the sum
- of seasonalities, holidays, and added regressors.
+A list with items "trend" and "yhat" containing
+ posterior predictive samples for that component.
 }
 \description{
 Sample from the posterior predictive distribution.
diff --git a/R/man/prophet.Rd b/R/man/prophet.Rd
index dc73bbf..1f085d5 100644
--- a/R/man/prophet.Rd
+++ b/R/man/prophet.Rd
@@ -7,10 +7,10 @@
 prophet(df = NULL, growth = "linear", changepoints = NULL,
   n.changepoints = 25, yearly.seasonality = "auto",
   weekly.seasonality = "auto", daily.seasonality = "auto",
-  holidays = NULL, seasonality.prior.scale = 10,
-  holidays.prior.scale = 10, changepoint.prior.scale = 0.05,
-  mcmc.samples = 0, interval.width = 0.8, uncertainty.samples = 1000,
-  fit = TRUE, ...)
+  holidays = NULL, seasonality.mode = "additive",
+  seasonality.prior.scale = 10, holidays.prior.scale = 10,
+  changepoint.prior.scale = 0.05, mcmc.samples = 0, interval.width = 0.8,
+  uncertainty.samples = 1000, fit = TRUE, ...)
 }
 \arguments{
 \item{df}{(optional) Dataframe containing the history. Must have columns ds
@@ -46,6 +46,8 @@ range of days around the date to be included as holidays. lower_window=-2
 will include 2 days prior to the date as holidays. Also optionally can have
 a column prior_scale specifying the prior scale for each holiday.}
 
+\item{seasonality.mode}{'additive' (default) or 'multiplicative'.}
+
 \item{seasonality.prior.scale}{Parameter modulating the strength of the
 seasonality model. Larger values allow the model to fit larger seasonal
 fluctuations, smaller values dampen the seasonality. Can be specified for
diff --git a/R/man/prophet_plot_components.Rd b/R/man/prophet_plot_components.Rd
index 4bf0f13..932a7de 100644
--- a/R/man/prophet_plot_components.Rd
+++ b/R/man/prophet_plot_components.Rd
@@ -3,8 +3,9 @@
 \name{prophet_plot_components}
 \alias{prophet_plot_components}
 \title{Plot the components of a prophet forecast.
-Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if
-present, and holidays if present.}
+Prints a ggplot2 with whichever are available of: trend, holidays, weekly
+seasonality, yearly seasonality, and additive and multiplicative extra
+regressors.}
 \usage{
 prophet_plot_components(m, fcst, uncertainty = TRUE, plot_cap = TRUE,
   weekly_start = 0, yearly_start = 0)
@@ -33,6 +34,7 @@ Invisibly return a list containing the plotted ggplot objects
 }
 \description{
 Plot the components of a prophet forecast.
-Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if
-present, and holidays if present.
+Prints a ggplot2 with whichever are available of: trend, holidays, weekly
+seasonality, yearly seasonality, and additive and multiplicative extra
+regressors.
 }
diff --git a/R/man/regressor_column_matrix.Rd b/R/man/regressor_column_matrix.Rd
new file mode 100644
index 0000000..2c0c534
--- /dev/null
+++ b/R/man/regressor_column_matrix.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prophet.R
+\name{regressor_column_matrix}
+\alias{regressor_column_matrix}
+\title{Dataframe indicating which columns of the feature matrix correspond to
+which seasonality/regressor components.}
+\usage{
+regressor_column_matrix(m, seasonal.features, modes)
+}
+\arguments{
+\item{m}{Prophet object.}
+
+\item{seasonal.features}{Constructed seasonal features dataframe.}
+
+\item{modes}{List with keys 'additive' and 'multiplicative' with arrays of
+component names for each mode of seasonality.}
+}
+\value{
+List with items
+ component.cols: A binary indicator dataframe with columns seasonal
+   components and rows columns in seasonal.features. Entry is 1 if that
+   column is used in that component.
+ modes: Updated input with combination components.
+}
+\description{
+Includes combination components, like 'additive_terms'. These combination
+components will be added to the 'modes' input.
+}
+\keyword{internal}
diff --git a/R/man/sample_model.Rd b/R/man/sample_model.Rd
index 4757a09..6bb5d9c 100644
--- a/R/man/sample_model.Rd
+++ b/R/man/sample_model.Rd
@@ -4,7 +4,7 @@
 \alias{sample_model}
 \title{Simulate observations from the extrapolated generative model.}
 \usage{
-sample_model(m, df, seasonal.features, iteration)
+sample_model(m, df, seasonal.features, iteration, s_a, s_m)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -14,9 +14,13 @@ sample_model(m, df, seasonal.features, iteration)
 \item{seasonal.features}{Data frame of seasonal features}
 
 \item{iteration}{Int sampling iteration to use parameters from.}
+
+\item{s_a}{Indicator vector for additive components}
+
+\item{s_m}{Indicator vector for multiplicative components}
 }
 \value{
-List of trend, seasonality, and yhat, each a vector like df$t.
+List of trend and yhat, each a vector like df$t.
 }
 \description{
 Simulate observations from the extrapolated generative model.
diff --git a/R/man/sample_posterior_predictive.Rd b/R/man/sample_posterior_predictive.Rd
index cf156b8..a3b4562 100644
--- a/R/man/sample_posterior_predictive.Rd
+++ b/R/man/sample_posterior_predictive.Rd
@@ -12,7 +12,8 @@ sample_posterior_predictive(m, df)
 \item{df}{Prediction dataframe.}
 }
 \value{
-List with posterior predictive samples for each component.
+List with posterior predictive samples for the forecast yhat and
+ for the trend component.
 }
 \description{
 Prophet posterior predictive samples.
diff --git a/R/tests/testthat/test_diagnostics.R b/R/tests/testthat/test_diagnostics.R
index 193232c..7f229f5 100644
--- a/R/tests/testthat/test_diagnostics.R
+++ b/R/tests/testthat/test_diagnostics.R
@@ -144,7 +144,8 @@ test_that("copy", {
     yearly.seasonality = c(TRUE, FALSE),
     weekly.seasonality = c(TRUE, FALSE),
     daily.seasonality = c(TRUE, FALSE),
-    holidays = c('null', 'insert_dataframe')
+    holidays = c('null', 'insert_dataframe'),
+    seasonality.mode = c('additive', 'multiplicative')
   )
   products <- expand.grid(inputs)
   for (i in 1:length(products)) {
@@ -178,7 +179,7 @@ test_that("copy", {
     args <- c('growth', 'changepoints', 'n.changepoints', 'holidays',
               'seasonality.prior.scale', 'holidays.prior.scale',
               'changepoints.prior.scale', 'mcmc.samples', 'interval.width',
-              'uncertainty.samples')
+              'uncertainty.samples', 'seasonality.mode')
     for (arg in args) {
       expect_equal(m1[[arg]], m2[[arg]])
     }
diff --git a/R/tests/testthat/test_prophet.R b/R/tests/testthat/test_prophet.R
index 5607f87..be0c127 100644
--- a/R/tests/testthat/test_prophet.R
+++ b/R/tests/testthat/test_prophet.R
@@ -225,7 +225,7 @@ test_that("piecewise_logistic", {
 })
 
 test_that("holidays", {
-  holidays = data.frame(ds = c('2016-12-25'),
+  holidays <- data.frame(ds = c('2016-12-25'),
                         holiday = c('xmas'),
                         lower_window = c(-1),
                         upper_window = c(0))
@@ -236,12 +236,14 @@ test_that("holidays", {
   out <- prophet:::make_holiday_features(m, df$ds)
   feats <- out$holiday.features
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_equal(nrow(feats), nrow(df))
   expect_equal(ncol(feats), 2)
   expect_equal(sum(colSums(feats) - c(1, 1)), 0)
   expect_true(all(priors == c(10., 10.)))
+  expect_equal(names, c('xmas'))
 
-  holidays = data.frame(ds = c('2016-12-25'),
+  holidays <- data.frame(ds = c('2016-12-25'),
                         holiday = c('xmas'),
                         lower_window = c(-1),
                         upper_window = c(10))
@@ -249,9 +251,11 @@ test_that("holidays", {
   out <- prophet:::make_holiday_features(m, df$ds)
   feats <- out$holiday.features
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_equal(nrow(feats), nrow(df))
   expect_equal(ncol(feats), 12)
   expect_true(all(priors == rep(10, 12)))
+  expect_equal(names, c('xmas'))
   # Check prior specifications
   holidays <- data.frame(
     ds = prophet:::set_date(c('2016-12-25', '2017-12-25')),
@@ -263,7 +267,9 @@ test_that("holidays", {
   m <- prophet(holidays = holidays, fit = FALSE)
   out <- prophet:::make_holiday_features(m, df$ds)
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_true(all(priors == c(5., 5.)))
+  expect_equal(names, c('xmas'))
   # 2 different priors
   holidays2 <- data.frame(
     ds = prophet:::set_date(c('2012-06-06', '2013-06-06')),
@@ -276,7 +282,9 @@ test_that("holidays", {
   m <- prophet(holidays = holidays2, fit = FALSE)
   out <- prophet:::make_holiday_features(m, df$ds)
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_true(all(priors == c(8, 8, 5, 5)))
+  expect_true(all(sort(names) == c('seans-bday', 'xmas')))
   holidays2 <- data.frame(
     ds = prophet:::set_date(c('2012-06-06', '2013-06-06')),
     holiday = c('seans-bday', 'seans-bday'),
@@ -334,7 +342,8 @@ test_that("auto_weekly_seasonality", {
   expect_equal(m$weekly.seasonality, 'auto')
   m <- fit.prophet(m, train.w)
   expect_true('weekly' %in% names(m$seasonalities))
-  true <- list(period = 7, fourier.order = 3, prior.scale = 10)
+  true <- list(
+    period = 7, fourier.order = 3, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$weekly[[name]], true[[name]])
   }
@@ -350,7 +359,8 @@ test_that("auto_weekly_seasonality", {
   m <- prophet(train.w)
   expect_false('weekly' %in% names(m$seasonalities))
   m <- prophet(DATA, weekly.seasonality = 2, seasonality.prior.scale = 3)
-  true <- list(period = 7, fourier.order = 2, prior.scale = 3)
+  true <- list(
+    period = 7, fourier.order = 2, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$weekly[[name]], true[[name]])
   }
@@ -363,7 +373,8 @@ test_that("auto_yearly_seasonality", {
   expect_equal(m$yearly.seasonality, 'auto')
   m <- fit.prophet(m, DATA)
   expect_true('yearly' %in% names(m$seasonalities))
-  true <- list(period = 365.25, fourier.order = 10, prior.scale = 10)
+  true <- list(
+    period = 365.25, fourier.order = 10, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$yearly[[name]], true[[name]])
   }
@@ -375,7 +386,8 @@ test_that("auto_yearly_seasonality", {
   m <- prophet(train.y, yearly.seasonality = TRUE)
   expect_true('yearly' %in% names(m$seasonalities))
   m <- prophet(DATA, yearly.seasonality = 7, seasonality.prior.scale = 3)
-  true <- list(period = 365.25, fourier.order = 7, prior.scale = 3)
+  true <- list(
+    period = 365.25, fourier.order = 7, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$yearly[[name]], true[[name]])
   }
@@ -388,7 +400,8 @@ test_that("auto_daily_seasonality", {
   expect_equal(m$daily.seasonality, 'auto')
   m <- fit.prophet(m, DATA2)
   expect_true('daily' %in% names(m$seasonalities))
-  true <- list(period = 1, fourier.order = 4, prior.scale = 10)
+  true <- list(
+    period = 1, fourier.order = 4, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$daily[[name]], true[[name]])
   }
@@ -400,7 +413,8 @@ test_that("auto_daily_seasonality", {
   m <- prophet(train.y, daily.seasonality = TRUE)
   expect_true('daily' %in% names(m$seasonalities))
   m <- prophet(DATA2, daily.seasonality = 7, seasonality.prior.scale = 3)
-  true <- list(period = 1, fourier.order = 7, prior.scale = 3)
+  true <- list(
+    period = 1, fourier.order = 7, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$daily[[name]], true[[name]])
   }
@@ -424,7 +438,8 @@ test_that("custom_seasonality", {
                          prior_scale = c(4))
   m <- prophet(holidays=holidays)
   m <- add_seasonality(m, name='monthly', period=30, fourier.order=5)
-  true <- list(period = 30, fourier.order = 5, prior.scale = 10)
+  true <- list(
+    period = 30, fourier.order = 5, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$monthly[[name]], true[[name]])
   }
@@ -436,12 +451,25 @@ test_that("custom_seasonality", {
   )
   m <- add_seasonality(m, name='weekly', period=30, fourier.order=5)
   # Test priors
-  m <- prophet(holidays = holidays, yearly.seasonality = FALSE)
+  m <- prophet(
+    holidays = holidays, yearly.seasonality = FALSE,
+    seasonality.mode = 'multiplicative')
   m <- add_seasonality(
-    m, name='monthly', period=30, fourier.order=5, prior.scale = 2)
+    m, name='monthly', period=30, fourier.order=5, prior.scale = 2,
+    mode = 'additive')
   m <- fit.prophet(m, DATA)
-  prior.scales <- prophet:::make_all_seasonality_features(
-    m, m$history)$prior.scales
+  expect_equal(m$seasonalities$monthly$mode, 'additive')
+  expect_equal(m$seasonalities$weekly$mode, 'multiplicative')
+  out <- prophet:::make_all_seasonality_features(m, m$history)
+  prior.scales <- out$prior.scales
+  component.cols <- out$component.cols
+  expect_equal(sum(component.cols$monthly), 10)
+  expect_equal(sum(component.cols$special_day), 1)
+  expect_equal(sum(component.cols$weekly), 6)
+  expect_equal(sum(component.cols$additive_terms), 10)
+  expect_equal(sum(component.cols$multiplicative_terms), 7)
+  expect_equal(sum(component.cols$monthly[1:11]), 10)
+  expect_equal(sum(component.cols$weekly[11:17]), 6)
   expect_true(all(prior.scales == c(rep(2, 10), rep(10, 6), 4)))
 })
 
@@ -450,10 +478,13 @@ test_that("added_regressors", {
   m <- prophet()
   m <- add_regressor(m, 'binary_feature', prior.scale=0.2)
   m <- add_regressor(m, 'numeric_feature', prior.scale=0.5)
+  m <- add_regressor(
+    m, 'numeric_feature2', prior.scale=0.5, mode = 'multiplicative')
   m <- add_regressor(m, 'binary_feature2', standardize=TRUE)
   df <- DATA
   df$binary_feature <- c(rep(0, 255), rep(1, 255))
   df$numeric_feature <- 0:509
+  df$numeric_feature2 <- 0:509
   # Require all regressors in df
   expect_error(
     fit.prophet(m, df)
@@ -461,7 +492,9 @@ test_that("added_regressors", {
   df$binary_feature2 <- c(rep(1, 100), rep(0, 410))
   m <- fit.prophet(m, df)
   # Check that standardizations are correctly set
-  true <- list(prior.scale = 0.2, mu = 0, std = 1, standardize = 'auto')
+  true <- list(
+    prior.scale = 0.2, mu = 0, std = 1, standardize = 'auto', mode = 'additive'
+  )
   for (name in names(true)) {
     expect_equal(true[[name]], m$extra_regressors$binary_feature[[name]])
   }
@@ -470,6 +503,7 @@ test_that("added_regressors", {
     expect_equal(true[[name]], m$extra_regressors$numeric_feature[[name]],
                  tolerance = 1e-5)
   }
+  expect_equal(m$extra_regressors$numeric_feature2$mode, 'multiplicative')
   true <- list(prior.scale = 10., mu = 0.1960784, std = 0.3974183)
   for (name in names(true)) {
     expect_equal(true[[name]], m$extra_regressors$binary_feature2[[name]],
@@ -484,28 +518,91 @@ test_that("added_regressors", {
   out <- prophet:::make_all_seasonality_features(m, df2)
   seasonal.features <- out$seasonal.features
   prior.scales <- out$prior.scales
-  expect_true('binary_feature' %in% colnames(seasonal.features))
-  expect_true('numeric_feature' %in% colnames(seasonal.features))
-  expect_true('binary_feature2' %in% colnames(seasonal.features))
-  expect_equal(ncol(seasonal.features), 29)
-  expect_true(all(sort(prior.scales[27:29]) == c(0.2, 0.5, 10.)))
+  component.cols <- out$component.cols
+  modes <- out$modes
+  expect_equal(ncol(seasonal.features), 30)
+  r_names <- c('binary_feature', 'numeric_feature', 'binary_feature2')
+  true.priors <- c(0.2, 0.5, 10.)
+  for (i in seq_along(r_names)) {
+    name <- r_names[i]
+    expect_true(name %in% colnames(seasonal.features))
+    expect_equal(sum(component.cols[[name]]), 1)
+    expect_equal(sum(prior.scales * component.cols[[name]]), true.priors[i])
+  }
   # Check that forecast components are reasonable
   future <- data.frame(
-    ds = c('2014-06-01'), binary_feature = c(0), numeric_feature = c(10))
+    ds = c('2014-06-01'),
+    binary_feature = c(0),
+    numeric_feature = c(10),
+    numeric_feature2 = c(10)
+  )
   expect_error(predict(m, future))
   future$binary_feature2 <- 0.
   fcst <- predict(m, future)
-  expect_equal(ncol(fcst), 31)
+  expect_equal(ncol(fcst), 37)
   expect_equal(fcst$binary_feature[1], 0)
-  expect_equal(fcst$extra_regressors[1],
+  expect_equal(fcst$extra_regressors_additive[1],
                fcst$numeric_feature[1] + fcst$binary_feature2[1])
-  expect_equal(fcst$seasonalities[1], fcst$yearly[1] + fcst$weekly[1])
-  expect_equal(fcst$seasonal[1],
-               fcst$seasonalities[1] + fcst$extra_regressors[1])
-  expect_equal(fcst$yhat[1], fcst$trend[1] + fcst$seasonal[1])
+  expect_equal(fcst$extra_regressors_multiplicative[1],
+               fcst$numeric_feature2[1])
+  expect_equal(fcst$additive_terms[1],
+               fcst$yearly[1] + fcst$weekly[1]
+               + fcst$extra_regressors_additive[1])
+  expect_equal(fcst$multiplicative_terms[1],
+               fcst$extra_regressors_multiplicative[1])
+  expect_equal(
+    fcst$yhat[1],
+    fcst$trend[1] * (1 + fcst$multiplicative_terms[1]) + fcst$additive_terms[1]
+  )
   # Check fails if constant extra regressor
   df$constant_feature <- 5
   m <- prophet()
   m <- add_regressor(m, 'constant_feature')
   expect_error(fit.prophet(m, df))
 })
+
+test_that("set_seasonality_mode", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  m <- prophet()
+  expect_equal(m$seasonality.mode, 'additive')
+  m <- prophet(seasonality.mode = 'multiplicative')
+  expect_equal(m$seasonality.mode, 'multiplicative')
+  expect_error(prophet(seasonality.mode = 'batman'))
+})
+
+test_that("seasonality_modes", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  holidays <- data.frame(ds = c('2016-12-25'),
+                        holiday = c('xmas'),
+                        lower_window = c(-1),
+                        upper_window = c(0))
+  m <- prophet(seasonality.mode = 'multiplicative', holidays = holidays)
+  m <- add_seasonality(
+    m, name = 'monthly', period = 30, fourier.order = 3, mode = 'additive')
+  m <- add_regressor(m, name = 'binary_feature', mode = 'additive')
+  m <- add_regressor(m, name = 'numeric_feature')
+  # Construct seasonal features
+  df <- DATA
+  df$binary_feature <- c(rep(0, 255), rep(1, 255))
+  df$numeric_feature <- 0:509
+  out <- prophet:::setup_dataframe(m, df, initialize_scales = TRUE)
+  df <- out$df
+  m <- out$m
+  m$history <- df
+  m <- prophet:::set_auto_seasonalities(m)
+  out <- prophet:::make_all_seasonality_features(m, df)
+  component.cols <- out$component.cols
+  modes <- out$modes
+  expect_equal(sum(component.cols$additive_terms), 7)
+  expect_equal(sum(component.cols$multiplicative_terms), 29)
+  expect_equal(
+    sort(modes$additive),
+    c('additive_terms', 'binary_feature', 'extra_regressors_additive',
+      'monthly')
+  )
+  expect_equal(
+    sort(modes$multiplicative),
+    c('extra_regressors_multiplicative', 'holidays', 'multiplicative_terms',
+      'numeric_feature', 'weekly', 'xmas', 'yearly')
+  )
+})
diff --git a/R/tests/testthat/test_stan_functions.R b/R/tests/testthat/test_stan_functions.R
index c3d8a13..9687096 100644
--- a/R/tests/testthat/test_stan_functions.R
+++ b/R/tests/testthat/test_stan_functions.R
@@ -1,9 +1,15 @@
 library(prophet)
 context("Prophet stan model tests")
 
-rstan::expose_stan_functions(
-  rstan::stanc(file="../..//inst/stan/prophet.stan")
-)
+fn <- tryCatch({
+  rstan::expose_stan_functions(
+    rstan::stanc(file="../../inst/stan/prophet.stan")
+  )
+}, error = function(e) {
+  rstan::expose_stan_functions(
+    rstan::stanc(file=system.file("stan/prophet.stan", package="prophet"))
+  )
+})
 
 DATA <- read.csv('data.csv')
 N <- nrow(DATA)
diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py
index aecf016..2852171 100644
--- a/python/fbprophet/forecaster.py
+++ b/python/fbprophet/forecaster.py
@@ -692,7 +692,7 @@ class Prophet(object):
         # Convert to a binary matrix
         component_cols = pd.crosstab(
             components['col'], components['component'],
-        )
+        ).sort_index(level='col')
         # Add columns for additive and multiplicative terms, if missing
         for name in ['additive_terms', 'multiplicative_terms']:
             if name not in component_cols:
@@ -729,7 +729,7 @@ class Prophet(object):
         new_comp = components[components['component'].isin(set(group))].copy()
         group_cols = new_comp['col'].unique()
         if len(group_cols) > 0:
-            new_comp = pd.DataFrame({'component': name, 'col': group_cols})
+            new_comp = pd.DataFrame({'col': group_cols, 'component': name})
             components = components.append(new_comp)
         return components
 
@@ -1248,6 +1248,8 @@ class Prophet(object):
         df: Prediction dataframe.
         seasonal_features: pd.DataFrame of seasonal features.
         iteration: Int sampling iteration to use parameters from.
+        s_a: Indicator vector for additive components
+        s_m: Indicator vector for multiplicative components
 
         Returns
         -------
diff --git a/python/fbprophet/plot.py b/python/fbprophet/plot.py
index 2b55723..535a128 100644
--- a/python/fbprophet/plot.py
+++ b/python/fbprophet/plot.py
@@ -344,7 +344,7 @@ def plot_seasonality(m, name, ax=None, uncertainty=True):
     ax.xaxis.set_major_formatter(FuncFormatter(
         lambda x, pos=None: fmt_str.format(dt=num2date(x))))
     ax.set_xlabel('ds')
-    ax.set_ylabel('{}'.format(name))
+    ax.set_ylabel(name)
     if m.seasonalities[name]['mode'] == 'multiplicative':
         ax = set_y_as_percent(ax)
     return artists