Marketing mix modeling in R: base R lm, adstock and optim, then Robyn

Marketing mix modeling in R needs nothing beyond base R: stats::filter computes geometric adstock, a one-line function gives you a Hill saturation curve, lm fits the regression of sales on transformed spend plus controls, and optim searches for the decay and saturation parameters that fit best. This tutorial builds that model step by step on a 104-week practice dataset, then shows where Meta’s Robyn package takes over.

R is where a lot of marketing mix modeling has historically lived, because econometricians and marketing scientists tend to work in it, and Robyn made it the default for many in-house teams. When I built the earlier Vexpower MMM path, Robyn was the first package students met after building a model by hand in a spreadsheet. However, I’d still start with base R, because Robyn makes dozens of decisions for you, and you’ll get much more out of its output once you’ve made those decisions yourself. The model here is the same one built in the Python tutorial and the Excel template, on the same data, so you can compare numbers across all three. For the concepts behind it, start with what marketing mix modeling is.

Marketing mix modeling in R starts with a data audit

The script expects mmm-practice-synthetic.csv in the working directory and downloads it from this site if it isn’t there. It’s an invented dataset: 104 weeks from January 2024 to December 2025 with revenue, spend on search, social and video, and a promotion flag. It was generated from a straight linear equation with no noise, no carryover and no saturation, so the adstock and Hill steps below are there to teach the method rather than to discover anything, and a correct fit should find decay near zero and push saturation to its bound. The first job is an audit, because a missing week or a column in the wrong units will break a model long before the statistics do.

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))

You should see 104 rows and six columns, no missing values, 103 gaps of seven days, and correlations with revenue of 0.05 for search, 0.19 for social, minus 0.10 for video and 0.22 for promotion. Those correlations are a useful warning. Search has the biggest budget and, as we’ll see, the biggest effect, but its raw correlation with revenue is close to zero because a yearly seasonal cycle moves revenue by more than search does. Only a regression with seasonality in it can separate the two, which is the confounding problem in miniature.

Geometric adstock with stats::filter

Adstock is the carryover of advertising from one week into the next. The geometric version says this week’s adstock is this week’s spend plus a fixed share (the decay) of last week’s adstock. That’s a recursive filter, and R has one built in: stats::filter with method = "recursive" computes y[t] = x[t] + decay * y[t-1], starting from zero, which is exactly what we need.

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)
}

Two small details matter. stats::filter returns a time-series object, so as.numeric turns it back into a plain vector. And it’s worth writing stats::filter in full every time, because if you’ve loaded dplyr, plain filter means something completely different and you’ll get a baffling error. The hill function is the Hill saturation curve with its shape fixed at 1, so it rises almost linearly at low spend and flattens towards a maximum, reaching half of it at the half-saturation point.

Build the controls and the design matrix

The regression needs controls for everything other than media that moves revenue: a trend counter, a sine and cosine pair with a 52-week period for seasonality, and the promotion flag. Each channel then gets one column, its spend adstocked and saturated. Half-saturation is expressed as a multiple of each channel’s average weekly spend, which keeps the three channels on the same scale for the optimizer.

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)
}

The sse function is the heart of it. For any six parameters it builds the design matrix, fits ordinary least squares with lm.fit (the fast engine underneath lm, which takes a matrix instead of a formula) and returns the residual sum of squares. Fitting the model then means finding the parameters that make that number as small as possible.

Fit the parameters with optim, then the model with lm

optim with method = "L-BFGS-B" minimizes a function within box bounds, which is what we want: decays between 0 and 0.9, and half-saturation multiples between 0.2 and 10. A coarse grid gives it a sensible place to start, so it doesn’t wander off into a bad local minimum.

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

You could do the same with nls and algorithm = "port", which also accepts bounds and gives you standard errors for the nonlinear parameters. In my experience it’s fussier about starting values and tends to fail exactly when a parameter wants to sit on its bound, which, as you’re about to see, is what happens here. With optim, the linear coefficients are handled by least squares inside the objective, so only the six nonlinear parameters are searched, and that’s much more stable.

Once the parameters are chosen, the script refits the final model with lm so you get the familiar summary table:

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))

The search should leave all three decays at or within a hair of zero and push all three half-saturation multiples to 10, the top of their range. That’s the model telling you there’s no detectable carryover and no detectable diminishing returns in this data, and it’s correct: the practice data was generated with neither. A parameter on its bound is always worth reporting, because it means your assumption, not the data, decided the value. The fit is essentially perfect, with an R-squared that rounds to 1.0000, a MAPE of 0.04 percent and a residual standard error of about 13.6, because the synthetic data has no noise beyond rounding. Real weekly revenue won’t fit like this, and an R-squared this close to 1 on real data should worry you.

Decomposition, ROI and marginal ROI

Multiplying each column of the design matrix by its coefficient splits every week’s revenue into baseline (intercept, trend and seasonality), promotion and each channel. sweep does that in one line. Note that lm names the intercept (Intercept), so the script renames it to match the design matrix before lining the coefficients up.

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))

Over the two years the model attributes 76.0 percent of revenue to the baseline, 14.5 percent to search, 5.8 percent to social, 3.1 percent to video and 0.5 percent to promotions. ROI is each channel’s total contribution divided by its total spend, and marginal ROI, the return on the next dollar at average spend, is the slope of the response curve at that point:

marginal_roi <- function(s, j) {
  a <- s / (1 - decays[j])
  media_beta[j] * half[j] / (a + half[j])^2 / (1 - decays[j])
}

You should get an average ROI of 1.54 for search, 0.87 for social and 0.65 for video, and marginal ROIs of 1.40, 0.80 and 0.60. These are revenue per dollar, so convert them to profit with your margin before calling anything profitable; the breakeven ROAS page does that. The gap between average and marginal is the reason budget decisions should use marginal ROAS.

Response curves and a budget reallocation

The response curves plot each channel’s steady-state weekly contribution against weekly spend, with a dot at today’s average. On this data they’re nearly straight across the observed range, which matches what the parameter search found. For the reallocation, the script keeps the total weekly budget fixed and lets each channel move up to 30 percent either way. With three channels the simplest reliable approach in base R is a fine grid: choose search and social spend on a 300 by 300 grid, give video whatever is left, and keep only the plans where video is inside its limits too.

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], ]

The best plan puts search at its ceiling of about 3,153 a week (up 728), social at about 1,353 (down 363) and video at about 855 (down 364), for roughly 474 a week of extra revenue at the same total spend. The Python version, using a continuous optimizer instead of a grid, lands within two dollars of the same split. On real data with genuine saturation the answer would stop short of the limits, at the point where the channels’ marginal returns are equal, and the budget allocation page explains why the limits matter so much.

Holdout check and the answer key

The last two steps refit the whole model (parameters and coefficients) on the first 78 weeks and predict the last 26, which gives a training MAPE of 0.03 percent and a holdout MAPE of 0.09 percent, and then compare the ROI estimates with the effects the data was generated from: 1.4, 0.8 and 0.6 of revenue per dollar of search, social and video. The marginal ROIs match exactly. The average ROIs are about 10 percent too high, because the half-saturation bound of 10 forces the Hill curve to bend slightly; widening the bound to 100 brings them to 1.41, 0.81 and 0.61. The model validation page covers what else to check before you trust a model, and the general lesson is that the bounds you set on saturation are assumptions that change the ROI you report.

The viewer below shows each step side by side in R, Python and Excel, using the code from the three downloads, so you can translate between them or explain the model to a colleague who works in a different tool.

d <- read.csv(csv)
d$week <- as.Date(d$week)

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

df = pd.read_csv("mmm-practice-synthetic.csv", parse_dates=["week"])

print(df.shape)
print("Missing values:", int(df.isna().sum().sum()))
print("Week gaps:", df["week"].diff().dropna().value_counts().to_dict())
print(df[["revenue"] + CHANNELS + ["promotion"]].corr().round(2)["revenue"])

Data sheet: one row per week, columns A:F = week, revenue, search_spend,
social_spend, video_spend, promotion (rows 2 to 105).

Missing values:   =COUNTBLANK(A2:F105)
Week gaps:        =A3-A2   (fill down; every value should be 7)
Correlation:      =CORREL($B$2:$B$105,C2:C105)

adstock <- function(x, decay) {
  as.numeric(stats::filter(x, filter = decay, method = "recursive"))
}

def geometric_adstock(x, decay):
    """adstock[t] = x[t] + decay * adstock[t-1]"""
    out = np.zeros(len(x))
    carry = 0.0
    for t, value in enumerate(x):
        carry = value + decay * carry
        out[t] = carry
    return out

Transform!B2:  =Data!C2
Transform!B3:  =Data!C3+Inputs!$C$6*B2      (fill down to B105)

Excel 365 / Google Sheets, one cell:
=SCAN(0,Data!C2:C105,LAMBDA(prev,x,x+Inputs!$C$6*prev))

hill <- function(x, half_sat) {
  x / (x + half_sat)
}

def hill(x, half_sat):
    """Share of the maximum effect reached at x. Equals 0.5 when x == half_sat."""
    return x / (x + half_sat)

Inputs!F6:     =D6*E6                     (multiple x average spend)
Transform!I2:  =B2/(B2+Inputs!$F$6)       (hill_search, fill down)

model <- lm(revenue ~ t + season_sin + season_cos + promotion +
              search_spend + social_spend + video_spend, data = model_data)
print(summary(model))

def fit_ols(params, rows=slice(None)):
    X = np.column_stack([controls, media_columns(params)])
    beta, *_ = np.linalg.lstsq(X[rows], y[rows], rcond=None)
    return beta, X

Model!B5:I9 (array):
=LINEST(Data!$B$2:$B$105,Transform!$E$2:$K$105,TRUE,TRUE)

One labelled coefficient (trend = 7th LINEST column):
=INDEX(LINEST(Data!$B$2:$B$105,Transform!$E$2:$K$105,TRUE,TRUE),1,7)

fit_params <- optim(start, sse, method = "L-BFGS-B", lower = lower, upper = upper)
best <- fit_params$par

result = minimize(sse, x0=start, bounds=bounds, method="L-BFGS-B")
best = result.x

Model!C22:  =SUMSQ(F32:F135)      (sum of squared residuals)

Solver: Set Objective $C$22, To: Min,
By Changing Inputs!$C$6:$C$8, subject to 0 <= decay <= 0.9,
Method: GRG Nonlinear

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"]
)

contrib = pd.DataFrame(X * beta, columns=control_names + CHANNELS, index=df["week"])
decomp = pd.DataFrame({
    "baseline": contrib[["intercept", "trend", "season_sin", "season_cos"]].sum(axis=1),
    "promotion": contrib["promotion"],
    "search": contrib["search_spend"],
    "social": contrib["social_spend"],
    "video": contrib["video_spend"],
})

Decomposition!B2 (baseline):
=Model!$B$13+Model!$C$13*Transform!E2+Model!$D$13*Transform!F2+Model!$E$13*Transform!G2
Decomposition!D2 (search):
=Model!$G$13*Transform!I2

marginal_roi <- function(s, j) {
  a <- s / (1 - decays[j])
  media_beta[j] * half[j] / (a + half[j])^2 / (1 - decays[j])
}
roi$roi <- round(roi$contribution / roi$spend, 2)

def marginal_roi(s, j):
    a = s / (1 - decays[j])
    return media_beta[j] * half[j] / (a + half[j]) ** 2 / (1 - decays[j])

roi["roi"] = (roi["contribution"] / roi["spend"]).round(2)

ROI!D4:  =C4/B4
ROI!E4:  =F4*H4/((Inputs!E6/(1-G4))+H4)^2/(1-G4)

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], ]

plan = minimize(
    neg_total_response,
    x0=mean_spend,
    bounds=[(0.7 * m, 1.3 * m) for m in mean_spend],
    constraints=[{"type": "eq", "fun": lambda s: s.sum() - total}],
    method="SLSQP",
).x

Solver on 'Response curves': Set Objective $C$12, To: Max,
By Changing $C$7:$C$9, subject to $C$10 = $B$10,
$C$7:$C$9 >= $E$7:$E$9 and $C$7:$C$9 <= $F$7:$F$9

train <- 1:78
test <- 79:n
fit_train <- optim(start, sse, rows = train, method = "L-BFGS-B", lower = lower, upper = upper)

train = slice(0, 78)
test = slice(78, n)
holdout_fit = minimize(sse, x0=start, args=(train,), bounds=bounds, method="L-BFGS-B")

=LINEST(Data!$B$2:$B$79,Transform!$E$2:$K$79,TRUE,TRUE)
then compare fitted and actual revenue for rows 80 to 105.

• You wrote stats::filter, not filter, so dplyr can't hijack it
• The audit found 103 seven-day gaps and no missing values (or you know why not)
• Seasonality and promotions are in the model as controls
• Any parameter sitting on its bound is reported, not hidden
• Holdout MAPE is close to training MAPE
• ROI has been converted to profit before it reaches a budget decision

Where Robyn takes over

Once you’ve built the model by hand, Robyn is much easier to understand, because it’s doing the same job with more machinery. Instead of OLS it fits a ridge regression, which shrinks coefficients to cope with correlated channels. Instead of a geometric decay and a fixed Hill shape it offers geometric or Weibull adstock and a Hill curve with a free shape. Instead of one optim run minimizing squared error, it uses Nevergrad (a Python optimization library, called from R) to search thousands of hyperparameter combinations against two objectives at once, prediction error (NRMSE) and how far each channel’s share of effect strays from its share of spend (DECOMP.RSSD), and hands you a Pareto front of candidate models to choose from. Experiment results can enter as a third objective, and the budget allocator does what our grid did, with more options.

To get started, install.packages("Robyn") installs the CRAN version (3.12.1, July 2025), and the repository’s demo.R walks through a full run on simulated data; you’ll also need Python with Nevergrad, which the README explains how to set up. Robyn models a single national series, so if you have regional data, look at the open-source MMM comparison for the Bayesian options, which are mostly in Python. For experiments, GeoLift and CausalImpact both run in R, and the geo lift testing page explains when you’d use them. If you’d like to practise the judgment calls rather than the code, the Vexpower MMM course works through a messier invented brand from brief to budget decision.