# Marketing mix modeling in R, from scratch (Vexpower)
# Guide: https://vexpower.com/mmm/marketing-mix-modeling-r
#
# Base R only: read.csv, stats::filter for adstock, lm for the regression, optim for the
# adstock and saturation parameters. Same model as the Python and Excel tutorials, on the
# same synthetic practice dataset, so you can compare numbers across all three.
#
# Run it from the folder that contains mmm-practice-synthetic.csv:
#   Rscript mmm-r-tutorial.R
# or source it in RStudio. If the CSV isn't there, it is downloaded from vexpower.com.
#
# The practice data is synthetic and exactly linear with no noise: revenue = 20000 + 1.4*search +
# 0.8*social + 0.6*video + 3000*sin(pi*t/26) + 1800*promotion, with no adstock and no saturation.
# Expect decay near zero and half-saturation on its upper bound; that is not a finding about media.

# ---- 1. Load and audit the data ---------------------------------------------------------

csv <- "mmm-practice-synthetic.csv"
if (!file.exists(csv)) {
  download.file("https://vexpower.com/downloads/mmm/mmm-practice-synthetic.csv", csv)
}
d <- read.csv(csv)
d$week <- as.Date(d$week)

channels <- c("search_spend", "social_spend", "video_spend")

print(dim(d))
cat("Missing values:", sum(is.na(d)), "\n")
print(table(diff(d$week)))              # every gap should be 7 days
print(range(d$week))
print(summary(d[, c("revenue", channels, "promotion")]))
print(round(cor(d[, c("revenue", channels, "promotion")])[, "revenue"], 2))

# ---- 2. Adstock and saturation ---------------------------------------------------------

# Geometric adstock: adstock[t] = x[t] + decay * adstock[t-1].
# stats::filter with method = "recursive" computes exactly this recursion.
# (Write stats::filter in full: dplyr::filter has the same name and does something else.)
adstock <- function(x, decay) {
  as.numeric(stats::filter(x, filter = decay, method = "recursive"))
}

# Hill curve with shape 1: share of the maximum effect reached at x; 0.5 when x == half_sat.
hill <- function(x, half_sat) {
  x / (x + half_sat)
}

# ---- 3. Controls and the design matrix -------------------------------------------------

n <- nrow(d)
d$t <- 0:(n - 1)
d$season_sin <- sin(2 * pi * d$t / 52)
d$season_cos <- cos(2 * pi * d$t / 52)

spend <- as.matrix(d[, channels])
mean_spend <- colMeans(spend)
y <- d$revenue

# params = c(decay_search, decay_social, decay_video, halfsat_search, halfsat_social, halfsat_video)
# Half-saturation is a multiple of each channel's average weekly spend.
media_matrix <- function(params) {
  sapply(1:3, function(j) {
    hill(adstock(spend[, j], params[j]), params[3 + j] * mean_spend[j])
  })
}

design_matrix <- function(params) {
  m <- media_matrix(params)
  colnames(m) <- channels
  cbind(intercept = 1, t = d$t, season_sin = d$season_sin, season_cos = d$season_cos,
        promotion = d$promotion, m)
}

sse <- function(params, rows = 1:n) {
  X <- design_matrix(params)[rows, , drop = FALSE]
  sum(lm.fit(X, y[rows])$residuals^2)
}

# ---- 4. Search for decay and half-saturation -------------------------------------------

lower <- c(0, 0, 0, 0.2, 0.2, 0.2)
upper <- c(0.9, 0.9, 0.9, 10, 10, 10)

# Coarse grid for a starting point (same decay and multiple for every channel) ...
grid <- expand.grid(decay = seq(0, 0.9, by = 0.1), mult = c(0.5, 1, 2, 5, 10))
grid_sse <- apply(grid, 1, function(g) sse(c(rep(g["decay"], 3), rep(g["mult"], 3))))
start <- unlist(c(rep(grid$decay[which.min(grid_sse)], 3), rep(grid$mult[which.min(grid_sse)], 3)))

# ... then optim polishes it within the bounds.
fit_params <- optim(start, sse, method = "L-BFGS-B", lower = lower, upper = upper)
best <- fit_params$par

print(data.frame(channel = channels,
                 decay = round(best[1:3], 3),
                 half_sat_x_mean = round(best[4:6], 2),
                 half_sat_dollars = round(best[4:6] * mean_spend)))

# Fit the final model with lm so you get the usual summary output.
model_data <- data.frame(revenue = y, design_matrix(best)[, -1])
model <- lm(revenue ~ t + season_sin + season_cos + promotion +
              search_spend + social_spend + video_spend, data = model_data)
print(summary(model))

fitted_rev <- fitted(model)
mape <- mean(abs((y - fitted_rev) / y))
cat(sprintf("R-squared: %.4f   MAPE: %.2f%%\n", summary(model)$r.squared, 100 * mape))

plot(d$week, y, type = "l", xlab = "Week", ylab = "Revenue", main = "Actual (black) vs model (green)")
lines(d$week, fitted_rev, col = "darkgreen", lty = 2)

# ---- 5. Decomposition ------------------------------------------------------------------

X <- design_matrix(best)
beta <- coef(model)
names(beta)[1] <- "intercept"      # lm calls it "(Intercept)"
beta <- beta[colnames(X)]          # same order as the design matrix
contrib <- sweep(X, 2, beta, `*`)
decomp <- data.frame(
  baseline = rowSums(contrib[, c("intercept", "t", "season_sin", "season_cos")]),
  promotion = contrib[, "promotion"],
  search = contrib[, "search_spend"],
  social = contrib[, "social_spend"],
  video = contrib[, "video_spend"]
)
print(round(colSums(decomp) / sum(decomp), 3))

barplot(t(as.matrix(decomp)), border = NA, space = 0, legend.text = colnames(decomp),
        args.legend = list(x = "topleft", ncol = 5, cex = 0.7, bty = "n"),
        main = "Weekly revenue decomposition", ylab = "Revenue")

# ---- 6. ROI and marginal ROI -----------------------------------------------------------

decays <- best[1:3]
half <- best[4:6] * mean_spend
media_beta <- beta[channels]

weekly_response <- function(s, j) {
  a <- s / (1 - decays[j])                    # steady-state adstock
  media_beta[j] * hill(a, half[j])
}
marginal_roi <- function(s, j) {
  a <- s / (1 - decays[j])
  media_beta[j] * half[j] / (a + half[j])^2 / (1 - decays[j])
}

roi <- data.frame(
  channel = channels,
  spend = colSums(spend),
  contribution = colSums(contrib[, channels])
)
roi$roi <- round(roi$contribution / roi$spend, 2)
roi$marginal_roi <- round(sapply(1:3, function(j) marginal_roi(mean_spend[j], j)), 2)
roi$spend <- round(roi$spend)
roi$contribution <- round(roi$contribution)
print(roi, row.names = FALSE)

# ---- 7. Response curves ----------------------------------------------------------------

s_max <- 3 * max(mean_spend)
plot(NULL, xlim = c(0, s_max), ylim = c(0, max(sapply(1:3, function(j) weekly_response(3 * mean_spend[j], j)))),
     xlab = "Weekly spend", ylab = "Weekly revenue contribution", main = "Response curves")
for (j in 1:3) {
  s <- seq(0, 3 * mean_spend[j], length.out = 100)
  lines(s, weekly_response(s, j), col = j, lwd = 2)
  points(mean_spend[j], weekly_response(mean_spend[j], j), col = j, pch = 19)
}
legend("topleft", legend = channels, col = 1:3, lwd = 2, bty = "n")

# ---- 8. Reallocate the same budget -----------------------------------------------------

# Same total weekly budget, each channel within +/- 30% of its average.
# With three channels a fine grid is simple and can't get stuck: choose search and social,
# video takes whatever is left, and keep only plans where video is inside its limits too.
total <- sum(mean_spend)
lo <- 0.7 * mean_spend
hi <- 1.3 * mean_spend
plans <- expand.grid(search = seq(lo[1], hi[1], length.out = 300),
                     social = seq(lo[2], hi[2], length.out = 300))
plans$video <- total - plans$search - plans$social
plans <- plans[plans$video >= lo[3] & plans$video <= hi[3], ]
plans$response <- weekly_response(plans$search, 1) + weekly_response(plans$social, 2) +
  weekly_response(plans$video, 3)
plan <- unlist(plans[which.max(plans$response), c("search", "social", "video")])

current_response <- sum(sapply(1:3, function(j) weekly_response(mean_spend[j], j)))
print(data.frame(channel = channels, current = round(mean_spend), proposed = round(plan),
                 change = round(plan - mean_spend)), row.names = FALSE)
cat(sprintf("Extra weekly revenue for the same budget: %.0f\n", max(plans$response) - current_response))

# ---- 9. Holdout check ------------------------------------------------------------------

train <- 1:78
test <- 79:n
fit_train <- optim(start, sse, rows = train, method = "L-BFGS-B", lower = lower, upper = upper)
X_all <- design_matrix(fit_train$par)
beta_train <- lm.fit(X_all[train, ], y[train])$coefficients
pred <- as.numeric(X_all %*% beta_train)
cat(sprintf("Train MAPE: %.2f%%   Holdout MAPE: %.2f%%\n",
            100 * mean(abs((y[train] - pred[train]) / y[train])),
            100 * mean(abs((y[test] - pred[test]) / y[test]))))

# ---- 10. Check against the answer key --------------------------------------------------

# The practice data was generated with known effects: 1.4, 0.8 and 0.6 of revenue per unit of
# search, social and video spend, with no carryover and no saturation. Half-saturation finished
# on its upper bound of 10, so refit with a much wider bound to see how much the bound shaped ROI.
wide <- optim(start, sse, method = "L-BFGS-B", lower = lower, upper = c(0.9, 0.9, 0.9, 100, 100, 100))
X_wide <- design_matrix(wide$par)
beta_wide <- lm.fit(X_wide, y)$coefficients
roi_wide <- colSums(sweep(X_wide[, channels], 2, beta_wide[channels], `*`)) / colSums(spend)
print(data.frame(channel = channels, roi_bound_10 = roi$roi, roi_bound_100 = round(roi_wide, 2),
                 true_roi = c(1.4, 0.8, 0.6)), row.names = FALSE)

# ---- Next: Robyn ----------------------------------------------------------------------
# install.packages("Robyn") gives you Meta's version of this workflow: ridge regression,
# geometric or Weibull adstock, Hill saturation, an automated hyperparameter search with
# Nevergrad (needs Python) and a budget allocator. See https://vexpower.com/mmm/robyn