Executive Summary
Where Net Revenue came from, and how much of the split can be believed.
The short answer
Across 156 weeks, marketing drove 36.2 percent of Net Revenue, with TV Spend the largest single contributor at 14.6 percent. Display Spend delivers the best return per dollar at 0.611, and all 4 channels are statistically distinguishable from zero. Channel budgets vary independently (highest VIF 1.03), so the per-channel split is stable.
The detail
The model tracks Net Revenue closely (R-squared 0.971) and splits it into 63.8 percent baseline, trend and seasonality and 36.2 percent marketing. TV Spend is the largest single contributor at 14.6 percent of total Net Revenue on 1,399,198 of spend. The best return per unit of spend is Display Spend at 0.611 (95 percent interval 0.516 to 0.706), and it is running at 81.7 percent of its own modelled ceiling. 4 of 4 channel coefficients are distinguishable from zero at the 5 percent level. Channel spends vary independently enough (highest variance inflation factor 1.03) that the per-channel split is stable. No channel's spend tracks the fitted baseline closely enough to suggest it is simply following demand.
What this can't tell you
These are associations fitted to spend already chosen, not experimental effects. The honest test of a causal claim is a holdout or geo experiment, which this analysis is not.
Analysis Overview
A media mix model of Net Revenue across 156 weekly periods and 4 channel(s).
The short answer
Marketing accounts for 36.2 percent of Net Revenue across 156 weeks, with the remaining 63.8 percent driven by baseline, trend, and seasonality. The model evaluated 376 candidate fits and tracks actual revenue movement with 97.1 percent accuracy, meaning the channel split is based on real variation in spend, not just statistical noise.
The detail
The model splits Net Revenue into a non-marketing baseline that absorbs 63.8 percent and 4 modelled channels that account for 36.2 percent. The baseline includes an intercept, a linear trend term, and annual seasonality (harmonic 1). Every channel's spend was first carried forward using a fitted geometric decay rate, then passed through a negative-exponential saturation curve, so the model credits spend beyond the week it was bought and stops crediting it at a constant rate once the channel is heavily loaded. Both shape parameters were fitted by grid search on the Akaike information criterion, evaluated across 376 candidate model fits over 156 weekly periods, then validated against the last 32 weeks held out of the fit.
What this can't tell you
This is an association fitted to spend already chosen for business reasons. A channel funded because demand was already rising will be credited with that demand. Only an experiment that varies spend on purpose can establish whether changing spend would move revenue by the amounts shown.
Data Preparation
How the timeline, the outcome and the spend columns were cleaned.
The short answer
All 156 weeks of data loaded cleanly with no rows dropped and no missing spend values. The weekly cadence was inferred from date gaps and has no internal gaps. Every spend column varied over time and was usable.
The detail
156 rows loaded; 156 weekly periods of Net Revenue were fitted. The cadence was inferred as weekly from the typical gap between dates in 'Week Starting'. The weekly timeline has no gaps. Every mapped spend column was usable. Spend columns were required to vary over time: a column with the same value in every period carries no information about what moves Net Revenue. 0 rows removed, 0 channels dropped.
What this can't tell you
Data quality was not a limiting factor. The analysis proceeds with the full observed period.
Where the Outcome Came From, Period by Period
The fitted split of Net Revenue into baseline and each channel.
The short answer
TV Spend is the largest marketing contributor at 14.6 percent of total Net Revenue, followed by Paid Search at 9.68 percent. The baseline—the revenue that would arrive with no marketing—carries 63.8 percent of the fitted total, which is why the channel blocks appear smaller than a raw spend-versus-sales chart would suggest.
The detail
Each bar is one week of Net Revenue split into the non-marketing baseline and each channel's fitted contribution. The baseline block carries 63.8 percent of the total and is the part that would have arrived with no marketing at all under this model—it is the reason the channel blocks are smaller than a raw spend-versus-sales chart would suggest. The largest marketing block belongs to TV Spend at 14.6 percent of the total. The blocks sum to the model's fitted value, not to the actual outcome; the model reproduces 97.1 percent of the movement, and the rest is residual.
What this can't tell you
The decomposition reflects fitted contributions, not causation. Residual error (2.9 percent) remains unexplained by the model.
Contribution, Return and Fitted Shape by Channel
Per-channel contribution, share, return per unit of spend, and the fitted carryover and saturation.
| Channel | Contribution | Share PCT | Spend | Return Per Spend | Marginal Return | Decay Rate | Saturation Point | PCT Of Ceiling | Significance |
|---|---|---|---|---|---|---|---|---|---|
| TV Spend | 4.664e+05 | 14.56 | 1.399e+06 | 0.3333 | 0.2118 | 0.6 | 2.808e+04 | 52.1 | p below 0.001 |
| Paid Search Spend | 3.101e+05 | 9.68 | 6.318e+05 | 0.4909 | 0.2069 | 0.1 | 8074 | 68.5 | p below 0.001 |
| Social Spend | 2.396e+05 | 7.48 | 4.324e+05 | 0.554 | 0.1857 | 0.29 | 3660 | 82.5 | p below 0.001 |
| Display Spend | 142120 | 4.44 | 2.328e+05 | 0.6106 | 0.2057 | 0.56 | 2021 | 81.7 | p below 0.001 |
The short answer
TV Spend dominates in absolute terms (14.6 percent of revenue on 1,399,197.9 spent), but Display Spend delivers the highest return per dollar at 0.6106. Social Spend is the most saturated at 82.5 percent of its ceiling, meaning additional budget there yields the least. All four channels show statistically significant effects.
The detail
Contribution is the coefficient multiplied by the channel's total transformed spend, net of carryover and saturation. TV Spend: 466,351.9 contribution, 14.56 percent share, 0.3333 return per spend, 0.60 decay rate (60.0 percent carry-forward), 52.1 percent of ceiling. Paid Search Spend: 310,132.4 contribution, 9.68 percent share, 0.4909 return per spend, 0.1 decay rate, 68.5 percent of ceiling. Social Spend: 239,553.1 contribution, 7.48 percent share, 0.554 return per spend, 0.29 decay rate, 82.5 percent of ceiling. Display Spend: 142,120 contribution, 4.44 percent share, 0.6106 return per spend, 0.56 decay rate, 81.7 percent of ceiling. All coefficients are distinguishable from zero (p below 0.001). TV Spend carries the most spend forward (decay rate 0.60). Social Spend is the most saturated, running at 82.5 percent of its modelled ceiling, so the model implies extra budget there buys the least.
What this can't tell you
These are associations with spend already chosen. A coefficient not distinguishable from zero would not be a number to spend against, but all four pass that test here.
Return per Unit of Spend, with 95 Percent Intervals
Modelled Net Revenue returned per unit of spend, by channel.
The short answer
All four channels show positive returns with tight confidence intervals that do not span zero. Display Spend's interval is widest (0.516 to 0.706), making its point estimate the least defensible; TV Spend's interval is tightest. The intervals exclude uncertainty in carryover and saturation, so true intervals are wider.
The detail
Display Spend: 0.6106 return per unit spend, 95 percent interval 0.516 to 0.706. Social Spend: 0.554 return per unit spend, 95 percent interval 0.4988 to 0.6093. Paid Search Spend: 0.4909 return per unit spend, 95 percent interval 0.4673 to 0.5145. TV Spend: 0.3333 return per unit spend, 95 percent interval 0.3169 to 0.3497. 0 of 4 channel intervals span zero, which means the data cannot rule out that any channel returns nothing. The widest interval belongs to Display Spend (0.516 to 0.706), so its point estimate is the least worth defending. These intervals come from the regression coefficients only. They do not include the uncertainty in the fitted carryover and saturation parameters, so the true intervals are wider than the ones drawn here.
What this can't tell you
Intervals omit shape parameter uncertainty, understating true confidence bounds. These are regression-based ranges, not experimental confidence.
Response Curves — Where Each Channel Flattens Out
Modelled outcome at each sustained spend level, per channel.
The short answer
Social Spend is the most saturated, at 82.5 percent of its ceiling at average spend levels, making it the flattest place to add budget. TV Spend is least saturated at 52.1 percent of its ceiling and has the most curve left. Curves beyond the highest observed spend are the model's shape assumption, not data evidence.
The detail
Each curve shows the modelled Net Revenue a channel would deliver per week if that spend level were sustained, once carryover has settled. The curves bend because the model uses a saturating response: the first units of spend buy more than the last. Social Spend is furthest along its curve, at 82.5 percent of its ceiling at an average of 2,772 per week—the flattest place to add budget. TV Spend is least saturated at 52.1 percent of its ceiling, so it has the most curve left. The curves are extrapolated to 1.5 times the highest spend ever observed for each channel; beyond the range you have actually spent, they are the model's shape assumption rather than anything the data has seen.
What this can't tell you
Extrapolation beyond observed spend range rests on the fitted saturation shape, not on empirical evidence of how channels behave at higher levels.
Can the Per-Channel Split Be Believed?
Variance inflation, pairwise spend correlation, and how closely each channel tracks the baseline.
| Channel | Vif | Max Pair Correlation | Baseline Correlation | Curve Span | Verdict |
|---|---|---|---|---|---|
| TV Spend | 1.03 | 0.048 | 0.023 | 0.648 | separately identified |
| Display Spend | 1.03 | 0.04 | -0.079 | 0.867 | separately identified |
| Paid Search Spend | 1.01 | 0.048 | -0.034 | 0.915 | separately identified |
| Social Spend | 1.01 | 0.045 | 0.014 | 0.84 | separately identified |
The short answer
Channel budgets vary independently (highest VIF 1.03), so the per-channel split is stable and trustworthy. Every channel traverses at least 0.50 of its saturation curve, confirming that contribution levels are identified by the data. No channel's spend tracks the baseline closely, ruling out simple demand-following.
The detail
The variance inflation factor measures how much a channel's coefficient uncertainty is multiplied by because other channels explain its movement. The highest here is TV Spend at 1.03. Channel spends vary independently enough (highest variance inflation factor 1.03) that the per-channel split is stable. The curve span column is a second identification check: it is how much of its own saturation curve a channel actually traverses across the observed spend range, on a 0 to 1 scale. TV Spend: 0.648 curve span, separately identified. Display Spend: 0.867 curve span, separately identified. Paid Search Spend: 0.915 curve span, separately identified. Social Spend: 0.84 curve span, separately identified. Every channel here spans at least 0.50 of its curve, so the levels are identified by the data rather than by the curve shape. The baseline correlation column is a third warning light: a channel whose spend tracks the fitted baseline is a channel whose budget follows demand. TV Spend baseline correlation 0.023, Display Spend -0.079, Paid Search Spend -0.034, Social Spend 0.014. No channel's spend tracks the fitted baseline closely enough to suggest it is simply following demand.
What this can't tell you
These checks confirm the split is stable within this data. They do not establish causation or rule out unmeasured confounding.
Model Fit Diagnostics
How well the model reproduces the outcome, and where its fit statistics flatter it.
| Metric | Value | Interpretation |
|---|---|---|
| R-squared | 0.971 | share of the period-to-period movement in Net Revenue the model reproduces |
| Adjusted R-squared | 0.970 | the same, penalised for the 8 fitted terms |
| Periods used | 156 | observed weekly periods after cleaning |
| Parameters fitted | 8 | intercept, trend, 2 seasonal term(s) and 4 channel coefficient(s) |
| Durbin-Watson | 1.96 | residuals show no strong autocorrelation |
| Highest channel VIF | 1.03 | collinearity between channels — the per-channel split is stable |
| Holdout RMSE (last 32 weeks) | 347.1 | typical error on week periods never used to fit the model |
| Holdout MAPE | 1.4 percent | the same error as a share of the actual outcome |
The short answer
The model reproduces 97.1 percent of Net Revenue movement with no strong residual autocorrelation (Durbin-Watson 1.96), and holdout error on unseen weeks is 1.4 percent of actual revenue. High R-squared on an MMM reflects working controls (trend, seasonality) rather than proof of the channel split; collinearity diagnostics are the decisive test.
The detail
The model reproduces 97.1 percent of the period-to-period movement in Net Revenue (97.0 percent after penalising the 8 fitted terms). 156 periods used, 8 parameters fitted (intercept, trend, 2 seasonal terms, 4 channel coefficients). The Durbin-Watson statistic of 1.96 shows no strong residual autocorrelation, so the intervals are not obviously understated. Holdout RMSE on the last 32 weeks never used to fit the model is 347.1, equivalent to 1.4 percent mean absolute percentage error. A high R-squared on a media mix model is easy to reach—trend and seasonality alone usually explain most of a business series—so it is evidence that the controls are working, not evidence that the channel split is right.
What this can't tell you
In-sample fit always flatters a model with 8 fitted terms. The collinearity card is the one that speaks to whether the per-channel split is believable.
Methods and Disclosure
The exact model form, the search, and what the numbers do and do not license.
| Item | Detail |
|---|---|
| Model form | Net Revenue in period t is modelled as a baseline plus a linear trend plus seasonality plus, for each channel, a coefficient times a saturated, carried-over version of that channel's spend. |
| Carryover (adstock) | Geometric: carried spend in period t equals this period's spend plus a decay rate times the previous period's carried spend. The decay rate was searched over 0.0 to 0.8 and fitted per channel, not assumed. |
| Diminishing returns (saturation) | Negative exponential: effect equals 1 minus the exponential of minus carried spend divided by a saturation scale. Concave everywhere, so the model can never imply that spend keeps paying back at a constant rate. |
| Baseline controls | Intercept, a linear trend term, and Fourier seasonality (annual cycle, harmonic 1). Without these every channel absorbs whatever the business was going to do anyway. |
| Parameter search | Coordinate-wise grid search: 9 decay values times up to 7 saturation scales per channel, repeated over 2 pass(es), 376 model fits evaluated. A further 128 candidate scales were rejected before fitting because they compressed a channel into less than 0.50 of its response curve, where the coefficient stops being separable from the intercept. |
| Selection criterion | Chosen on the Akaike information criterion over all 156 fitted weekly periods, with the shape parameters counted as parameters, then checked against the last 32 weeks held out of the fit. |
| Estimation | Ordinary least squares on 156 observed weekly periods with 8 fitted terms. No regularisation is applied: shrinking the coefficients would make an entangled per-channel split look calmer than it is, and this analysis reports the entanglement instead. |
| Uncertainty | 95 percent intervals come from the coefficient standard errors and scale straight through to contribution, return per unit of spend and marginal return, because each is that coefficient multiplied by a fixed quantity. They do NOT include the uncertainty in the fitted decay and saturation parameters, so they are narrower than the truth. |
| Marginal return | Computed by lifting each channel's entire spend path by 1 percent, re-running that channel's own fitted carryover and saturation, and dividing the change in modelled outcome by the change in spend. |
| Collinearity diagnostic | Variance inflation factor per channel, computed by regressing each channel's transformed spend on the other model terms. Above 10, the channels move together too closely for their individual coefficients to be separated. |
| Neighbouring questions | This answers how much of Net Revenue each channel's SPEND is associated with over time. It is not an incrementality test, which measures the lift of one specific intervention against markets deliberately held back; and it is not multi-touch attribution, which splits credit for conversions that already happened across the touchpoints inside each individual journey. The three answer different questions and will not agree. |
| What this is not | This is a model fitted to observed spend, not an experiment. Channels whose budgets follow demand will be credited with demand. The honest test of a causal claim about Net Revenue is a holdout or geo experiment in which spend is deliberately varied; this analysis is not one. |
The short answer
The model is ordinary least squares on 156 weeks, using geometric adstock (decay rate fitted per channel from 0.0 to 0.8) and negative-exponential saturation (concave, so spend never pays back at a constant rate). Grid search evaluated 376 candidate fits and chose the best by Akaike information criterion, then validated on 32 held-out weeks. Intervals omit shape parameter uncertainty, understating true confidence bounds.
The detail
Model form: Net Revenue in period t equals a baseline plus a linear trend plus seasonality plus, for each channel, a coefficient times a saturated, carried-over version of that channel's spend. Carryover (adstock) is geometric: carried spend in period t equals this period's spend plus a decay rate times the previous period's carried spend. The decay rate was searched over 0.0 to 0.8 and fitted per channel, not assumed. Diminishing returns (saturation) uses negative exponential: effect equals 1 minus the exponential of minus carried spend divided by a saturation scale. Concave everywhere, so the model can never imply that spend keeps paying back at a constant rate. Baseline controls: intercept, a linear trend term, and Fourier seasonality (annual cycle, harmonic 1). Parameter search: coordinate-wise grid search, 9 decay values times up to 7 saturation scales per channel, repeated over 2 passes, 376 model fits evaluated. Selection criterion: Akaike information criterion over all 156 fitted weekly periods, with shape parameters counted as parameters, then checked against the last 32 weeks held out of the fit. Estimation: ordinary least squares on 156 observed weekly periods with 8 fitted terms. No regularisation applied. 95 percent intervals come from coefficient standard errors and do NOT include the uncertainty in the fitted decay and saturation parameters, so they are narrower than the truth. A channel funded because demand was already rising will be credited with that demand.
What this can't tell you
Only an experiment that varies spend on purpose can establish whether changing spend would change Net Revenue by the amounts shown. This analysis fitted to spend already chosen for business reasons.
Methodology
Statistical methodology and diagnostics for Media Mix Model — Lite
Statistical Method
Which marketing channels actually drive sales? Regression-based media mix attribution with adstock and saturation. Map a dated outcome (sales, revenue, conversions) and the spend columns for each channel, and get a per-channel contribution and share of the outcome, a return and marginal return per unit of spend with 95 percent intervals, response curves showing where each channel flattens out, the whole outcome decomposed over time into a non-marketing baseline plus each channel, and a collinearity diagnostic that says out loud when the per-channel split cannot be trusted. Carryover decay and diminishing returns are FITTED per channel by grid search plus refinement, never assumed.
- Each row is one time period, and the outcome and spend columns cover the same period
- The cadence is regular enough to infer (daily, weekly or monthly); duplicate periods are summed
- Advertising effect carries over geometrically and saturates as a negative exponential — the shapes are fitted, but the FAMILY of shapes is an assumption
- The non-marketing baseline is well described by an intercept, a linear trend and Fourier seasonality
- Nothing outside the mapped columns changed at the same time as the spend (no concurrent price change, stockout, competitor exit or PR shock)
- Spend was not set as a direct response to the outcome within the same period
- This is a model fitted to spend the business already chose, not an experiment: a channel funded because demand was rising will be credited with that demand
- When channels' budgets move together, individual coefficients are unstable even though the total fit is good — the collinearity diagnostic reports this and the per-channel split should then not be used
- A channel that never moves far along its own saturation curve has a contribution LEVEL that leans on the assumed curve shape rather than on observed variation
- The 95 percent intervals come from the regression coefficients only; they treat the fitted carryover and saturation parameters as known exactly, so they are narrower than the truth
Analysis Code
Complete R source code for this analysis
Media Mix Model — Lite
Estimates how much of a dated outcome (sales, revenue, conversions) each marketing channel's spend is associated with, using the three things that separate a media mix model from a naive regression on raw spend: fitted carryover (adstock), fitted diminishing returns (saturation), and an explicit non-marketing baseline with trend and seasonality controls.
Why This Method?
Regressing an outcome on raw weekly spend gets two things wrong at once. It assumes advertising works only in the period it is bought (it does not — it decays), and it assumes the tenth unit of spend works as hard as the first (it does not — it saturates). Without a baseline, trend and seasonality, every channel also absorbs whatever the business was going to do anyway, and every channel looks effective.
What This Analysis Covers
- Per-channel carryover rate and saturation point, fitted by grid search
- Contribution and share of the outcome, with the non-marketing baseline
- Return per unit of spend and marginal return, each with a 95% interval
- Response curves showing where each channel flattens out
- The decomposition over time as a stacked view
- Model fit diagnostics and a collinearity diagnostic that says out loud
when the per-channel split cannot be trusted
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {date, outcome, spend_1..spend_N}. All narrative is derived from the user's own column names and computed values.
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))Helpers
Core Analysis Pipeline
compute_shared <- function(df, params, col_map = list()) {
# === SHARED EXPORTS ===
# initial_rows/final_rows/rows_removed $ row accounting
# date_name / outcome_name $ humanized single-word keys
# channel_names $ named chr, semantic -> user name
# used_channels / dropped_df $ kept vs excluded spend columns
# cadence / period_word / n_periods $ time grid
# n_missing_periods / n_dup_rows / n_blank_spend / n_missing_outcome
# par $ fitted theta/kappa per channel
# selection_rule / selection_label $ holdout vs information criterion
# channel_df $ channel, contribution, share_pct, spend,
# return_per_spend, marginal_return, decay_rate,
# saturation_point, p_value (card dataset)
# roi_df $ channel, return_per_spend, ci_low, ci_high
# decomp_df $ period(ISO), component, contribution
# curves_df $ channel, spend_level, response
# collinearity_df $ channel, vif, max_pair_correlation,
# baseline_correlation, verdict
# fit_df_out $ metric, value, interpretation
# methods_df $ item, detail
# r2 / adj_r2 / dw / holdout_rmse / holdout_mape
# vif_max / collin_verdict / collin_phrase
# neg_channels / demand_followers / aliased_channels
# metrics / json_output
# === /SHARED EXPORTS ===
initial_rows <- nrow(df)
date_name <- humanize_semantic("date", col_map)[1]
outcome_name <- humanize_semantic("outcome", col_map)[1]Step 1: Find the mapped columns
if (!("date" %in% names(df))) {
stop(sprintf("A dated period column must be mapped(expected '%s') — a media mix model needs one row per time period.",
date_name))
}
if (!("outcome" %in% names(df))) {
stop(sprintf("An outcome column must be mapped(expected '%s') — sales, revenue or conversions per period.",
outcome_name))
}
sp_cols <- grep("^spend_[0-9]+$", names(df), value = TRUE)
sp_cols <- sp_cols[order(as.integer(sub("^spend_", "", sp_cols)))]
if (length(sp_cols) < 1) {
stop("At least one marketing-spend column must be mapped — there is nothing to attribute the outcome to.")
}
channel_names <- setNames(humanize_semantic(sp_cols, col_map), sp_cols)Step 2: Coerce the outcome (95% rule) and the spend columns
yv <- df$outcome
if (!is.numeric(yv)) {
conv <- suppressWarnings(as.numeric(as.character(yv)))
n_orig <- sum(!is.na(yv) & nzchar(trimws(as.character(yv))))
if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
yv <- conv
} else {
stop(sprintf("The outcome column '%s' is not numeric — a media mix model needs a numeric outcome such as sales, revenue or conversions.",
outcome_name))
}
}
yv <- as.numeric(yv)
dropped_names <- character(0)
dropped_reason <- character(0)
spend_raw <- list()
n_blank_spend <- 0L
for (sc in sp_cols) {
v <- df[[sc]]
if (!is.numeric(v)) {
conv <- suppressWarnings(as.numeric(as.character(v)))
n_orig <- sum(!is.na(v) & nzchar(trimws(as.character(v))))
if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
v <- conv
} else {
dropped_names <- c(dropped_names, channel_names[[sc]])
dropped_reason <- c(dropped_reason, "not numeric")
next
}
}
v <- as.numeric(v)
n_blank_spend <- n_blank_spend + sum(is.na(v))
v[is.na(v)] <- 0 # a blank spend cell is read as no spend that period
spend_raw[[sc]] <- v
}
if (length(spend_raw) < 1) {
stop(sprintf("None of the mapped spend columns(%s) could be read as numbers.",
paste(channel_names[sp_cols], collapse = ", ")))
}Step 3: Parse the dates
raw_date <- df$date
non_blank <- !(is.na(raw_date) | !nzchar(trimws(as.character(raw_date))))
d <- parse_dates_robust(raw_date)
n_nonblank <- sum(non_blank)
if (n_nonblank == 0) {
stop(sprintf("The period column '%s' is empty — there are no dates to model over.", date_name))
}
n_unparsed <- sum(non_blank & is.na(d))
if (n_unparsed / n_nonblank > 0.05) {
stop(sprintf("%d of %d values in '%s' could not be read as dates. Please use a recognizable date format, for example 2024-01-31 or 01/31/2024.",
n_unparsed, n_nonblank, date_name))
}Step 4: Keep complete rows, aggregate duplicate periods by sum
keep <- !is.na(d) & !is.na(yv)
n_missing_outcome <- sum(!is.na(d) & is.na(yv))
if (sum(keep) < 2) {
stop(sprintf("Fewer than two periods have both a readable '%s' and a '%s' value.",
date_name, outcome_name))
}
work <- data.frame(date = d[keep], outcome = yv[keep], stringsAsFactors = FALSE)
for (sc in names(spend_raw)) work[[sc]] <- spend_raw[[sc]][keep]
n_dup_rows <- nrow(work) - length(unique(work$date))
val_cols <- setdiff(names(work), "date")
agg <- aggregate(work[, val_cols, drop = FALSE],
by = list(date = work$date), FUN = sum)
agg <- agg[order(agg$date), , drop = FALSE]
rownames(agg) <- NULL
if (nrow(agg) < 2) {
stop(sprintf("Only one distinct period in '%s' after cleaning — a media mix model needs a time series.",
date_name))
}Step 5: Infer the cadence and lay a regular time grid
med_gap <- median(as.numeric(diff(agg$date)))
if (med_gap <= 1.5) {
cadence <- "daily"; period_word <- "day"; step_days <- 1
season_periods <- c(7, 365.25); season_labels <- c("weekly", "annual")
theta_grid <- seq(0, 0.9, by = 0.1)
} else if (med_gap >= 5.5 && med_gap <= 8.5) {
cadence <- "weekly"; period_word <- "week"; step_days <- 7
season_periods <- c(365.25 / 7); season_labels <- "annual"
theta_grid <- seq(0, 0.8, by = 0.1)
} else if (med_gap >= 26 && med_gap <= 35) {
cadence <- "monthly"; period_word <- "month"; step_days <- NA
season_periods <- c(12); season_labels <- "annual"
theta_grid <- seq(0, 0.8, by = 0.1)
} else {
cadence <- "irregular"; period_word <- sprintf("%d-day period", max(1, round(med_gap)))
step_days <- max(1, round(med_gap))
season_periods <- numeric(0); season_labels <- character(0)
theta_grid <- seq(0, 0.8, by = 0.1)
}
if (cadence == "monthly") {
snapped <- as.Date(format(agg$date, "%Y-%m-01"))
agg <- aggregate(agg[, val_cols, drop = FALSE], by = list(date = snapped), FUN = sum)
agg <- agg[order(agg$date), , drop = FALSE]
grid_dates <- seq(min(agg$date), max(agg$date), by = "month")
} else {
grid_dates <- seq(min(agg$date), max(agg$date), by = step_days)
}
n_grid <- length(grid_dates)
pos <- match(agg$date, grid_dates)
# Periods whose date does not land on the grid (irregular reporting) are
# snapped to the nearest grid period rather than discarded.
if (any(is.na(pos))) {
for (i in which(is.na(pos))) {
pos[i] <- safe_which_max(-abs(as.numeric(grid_dates - agg$date[i])))
}
if (anyNA(pos)) {
agg <- agg[!is.na(pos), , drop = FALSE]
pos <- pos[!is.na(pos)]
if (!nrow(agg)) stop(sprintf("No period in '%s' could be placed on the inferred timeline.", date_name))
}
}
y_grid <- rep(NA_real_, n_grid)
spend_grid <- lapply(names(spend_raw), function(sc) rep(0, n_grid))
names(spend_grid) <- names(spend_raw)
for (i in seq_along(pos)) {
k <- pos[i]
y_grid[k] <- if (is.na(y_grid[k])) agg$outcome[i] else y_grid[k] + agg$outcome[i]
for (sc in names(spend_grid)) spend_grid[[sc]][k] <- spend_grid[[sc]][k] + agg[[sc]][i]
}
obs_idx <- which(!is.na(y_grid))
n_obs <- length(obs_idx)
n_missing_periods <- n_grid - n_obs
if (n_grid > 0 && n_missing_periods / n_grid > 0.3) {
stop(sprintf("%d of the %d %s periods between the first and last '%s' have no data (%s of the timeline). The series is too broken to model.",
n_missing_periods, n_grid, cadence, date_name,
fmt_pct(100 * n_missing_periods / n_grid)))
}Step 6: Drop spend columns with no variation over the observed periods
used_channels <- character(0)
for (sc in names(spend_grid)) {
v <- spend_grid[[sc]][obs_idx]
if (all(!is.finite(v)) || isTRUE(all(v == v[1])) || isTRUE(sd(v) == 0) || is.na(sd(v))) {
dropped_names <- c(dropped_names, channel_names[[sc]])
dropped_reason <- c(dropped_reason, "constant — no spend variation to learn from")
} else if (all(v <= 0)) {
dropped_names <- c(dropped_names, channel_names[[sc]])
dropped_reason <- c(dropped_reason, "no positive spend recorded")
} else {
used_channels <- c(used_channels, sc)
}
}
if (length(used_channels) < 1) {
stop(sprintf("None of the mapped spend columns(%s) vary over time, so none of them can explain movement in '%s'.",
paste(channel_names[sp_cols], collapse = ", "), outcome_name))
}Step 7: Build the non-marketing controls — trend and seasonality
tt <- seq_len(n_grid) - 1L
ctrl <- matrix(as.numeric(tt), ncol = 1)
ctrl_names <- "trend"
seasonal_terms <- character(0)
for (si in seq_along(season_periods)) {
P <- season_periods[si]
if (P <= 1 || n_obs < 1.5 * P) next
K <- if (n_obs >= 3 * P) 2L else 1L
for (k in seq_len(K)) {
ctrl <- cbind(ctrl, sin(2 * pi * k * tt / P), cos(2 * pi * k * tt / P))
nm <- sprintf("seas_%s_%d", season_labels[si], k)
ctrl_names <- c(ctrl_names, paste0(nm, "_sin"), paste0(nm, "_cos"))
seasonal_terms <- c(seasonal_terms, sprintf("%s cycle, harmonic %d", season_labels[si], k))
}
}
colnames(ctrl) <- ctrl_names
n_params <- 1L + ncol(ctrl) + length(used_channels)
min_needed <- max(24L, as.integer(3 * n_params))
if (n_obs < min_needed) {
stop(sprintf("Only %d %s periods of '%s' are usable. Fitting carryover, saturation and a baseline for %d channel(s) needs at least %d periods.",
n_obs, cadence, outcome_name, length(used_channels), min_needed))
}Step 8: Fit carryover and saturation
Two stages, and the decay rates are FITTED in both, never assumed. (a) A coordinate-wise GRID search: each channel's (decay, saturation scale) pair is scored against the whole model with the other channels held at their current values, repeated until no channel improves. (b) A continuous Nelder-Mead REFINEMENT started from the grid solution, because the grid is deliberately coarse and the criterion surface between grid points is not flat. The criterion in both stages is the Akaike information criterion, counting the two fitted shape parameters per channel as parameters. When the series is long enough, the last fifth of it is additionally held back and scored out-of-sample as an INDEPENDENT check — reported, never selected on.
y_obs <- y_grid[obs_idx]
h <- max(6L, as.integer(ceiling(0.2 * n_obs)))
use_holdout <- (n_obs >= 30L) && ((n_obs - h) >= (3L * n_params))
train_i <- seq_len(n_obs); test_i <- if (use_holdout) (n_obs - h + 1L):n_obs else integer(0)
selection_rule <- "aic"
selection_label <- sprintf(
"the Akaike information criterion over all %d fitted %s periods, with the shape parameters counted as parameters%s",
n_obs, cadence,
if (use_holdout) sprintf(", then checked against the last %d %ss held out of the fit", h, period_word) else "")Saturation scales are proposed as multiples of the channel's own mean positive carried spend, then FILTERED: a scale so small that every active period sits on the flat top of the curve turns the regressor into a near constant, and a near-constant regressor's coefficient is not separable from the intercept — it explodes while the intercept absorbs the offset, inflating the channel's contribution with no change in fit. A candidate must therefore let the channel traverse at least half of its response curve across the observed periods; if none does, the widest-spanning candidate is kept and the channel is flagged as weakly identified.
MIN_TRANSFORM_SPAN <- 0.5
kappa_candidates <- function(a) {
p <- a[a > 0]
if (!length(p)) return(list())
m <- mean(p)
if (!is.finite(m) || m <= 0) return(list())
ks <- sort(unique(m * c(0.25, 0.4, 0.6, 0.85, 1.2, 1.8, 3.0)))
out <- vector("list", length(ks)); spans <- numeric(length(ks))
for (i in seq_along(ks)) {
z <- saturate_negexp(a, ks[i])
zo <- z[obs_idx]
spans[i] <- if (length(zo) && all(is.finite(zo))) diff(range(zo)) else 0
out[[i]] <- list(kappa = ks[i], z = z)
}
ok <- which(spans >= MIN_TRANSFORM_SPAN)
if (!length(ok)) {
w <- safe_which_max(spans)
if (is.na(w)) return(list())
ok <- w
}
out[ok]
}
Zcur <- matrix(0, nrow = n_grid, ncol = length(used_channels),
dimnames = list(NULL, used_channels))
par <- list()
n_rejected_scales <- 0L
for (sc in used_channels) {
a0 <- adstock_geometric(spend_grid[[sc]], 0)
k0 <- kappa_candidates(a0)
if (!length(k0)) {
par[[sc]] <- list(theta = 0, kappa = max(1e-9, mean(a0[a0 > 0])))
Zcur[, sc] <- saturate_negexp(a0, par[[sc]]$kappa)
} else {
pick <- k0[[max(1L, as.integer(ceiling(length(k0) / 2)))]]
par[[sc]] <- list(theta = 0, kappa = pick$kappa)
Zcur[, sc] <- pick$z
}
}
n_shape <- 2L * length(used_channels)
score_design <- function(Z) {
if (any(!is.finite(Z))) return(Inf)
X <- cbind(`(Intercept)` = 1, ctrl, Z)[obs_idx, , drop = FALSE]
f <- tryCatch(stats::lm.fit(X, y_obs), error = function(e) NULL)
if (is.null(f)) return(Inf)
r <- as.numeric(f$residuals); nn <- length(r); rss <- sum(r^2)
if (!is.finite(rss) || rss <= 0) return(Inf)
nn * log(rss / nn) + 2 * (f$rank + 1L + n_shape)
}
best_score <- score_design(Zcur)
n_evals <- 0L
n_passes <- 0L
for (pass in seq_len(3L)) {
n_passes <- pass
improved <- FALSE
for (sc in used_channels) {
keep_col <- Zcur[, sc]
best_col <- keep_col
best_par <- par[[sc]]
for (th in theta_grid) {
a <- adstock_geometric(spend_grid[[sc]], th)
cands <- kappa_candidates(a)
n_rejected_scales <- n_rejected_scales + (7L - length(cands))
for (cd in cands) {
Zcur[, sc] <- cd$z
s <- score_design(Zcur)
n_evals <- n_evals + 1L
if (is.finite(s) && s < best_score - 1e-9) {
best_score <- s; best_col <- cd$z
best_par <- list(theta = th, kappa = cd$kappa); improved <- TRUE
}
}
}
Zcur[, sc] <- best_col
par[[sc]] <- best_par
}
if (!improved) break
}Stage (b): continuous refinement from the grid solution. Decay is optimised on a logistic scale bounded below 0.95 (a decay of 1 would mean spend never stops working) and the saturation scale on a log scale. The same minimum-span rule applies as a hard barrier, so the optimiser cannot walk into the region where a channel's coefficient stops being separable from the intercept.
THETA_MAX <- 0.95
span_of <- function(z) {
zo <- z[obs_idx]
if (!length(zo) || !all(is.finite(zo))) return(NA_real_)
diff(range(zo))
}Channels whose widest available candidate already falls short of the minimum span cannot be held to it — the barrier for those is their own best achievable span, so the refinement can still improve their fit instead of being frozen at an infeasible start. They stay flagged.
span_floor <- setNames(sapply(used_channels, function(sc) {
g <- span_of(Zcur[, sc])
if (is.na(g)) 0 else min(MIN_TRANSFORM_SPAN, g)
}), used_channels)
span_ok <- function(z, sc) {
sp <- span_of(z)
!is.na(sp) && sp >= span_floor[[sc]] - 1e-9
}
to_v <- function(pl) unlist(lapply(used_channels, function(sc) {
th <- min(max(pl[[sc]]$theta, 1e-4), THETA_MAX - 1e-4)
c(log(th / (THETA_MAX - th)), log(max(pl[[sc]]$kappa, 1e-9)))
}))
from_v <- function(v) {
out <- list()
for (i in seq_along(used_channels)) {
e <- exp(v[2 * i - 1])
out[[used_channels[i]]] <- list(theta = THETA_MAX * e / (1 + e),
kappa = exp(v[2 * i]))
}
out
}
n_refine <- 0L
obj_v <- function(v) {
if (any(!is.finite(v))) return(1e12)
pl <- from_v(v)
Z <- Zcur
for (sc in used_channels) {
z <- saturate_negexp(adstock_geometric(spend_grid[[sc]], pl[[sc]]$theta), pl[[sc]]$kappa)
if (!span_ok(z, sc)) return(1e12)
Z[, sc] <- z
}
n_refine <<- n_refine + 1L
s <- score_design(Z)
if (!is.finite(s)) 1e12 else s
}
v0 <- to_v(par)
if (all(is.finite(v0)) && is.finite(obj_v(v0))) {
op <- tryCatch(stats::optim(v0, obj_v, method = "Nelder-Mead",
control = list(maxit = 800, reltol = 1e-10)),
error = function(e) NULL)
if (!is.null(op) && is.finite(op$value) && op$value < best_score - 1e-9) {
par <- from_v(op$par)
for (sc in used_channels) {
Zcur[, sc] <- saturate_negexp(
adstock_geometric(spend_grid[[sc]], par[[sc]]$theta), par[[sc]]$kappa)
}
best_score <- op$value
}
}Step 9: Final fit on every observed period, with the chosen transforms
fit_frame <- as.data.frame(cbind(ctrl, Zcur)[obs_idx, , drop = FALSE])
ch_alias <- setNames(paste0("ch_", seq_along(used_channels)), used_channels)
names(fit_frame) <- c(ctrl_names, unname(ch_alias[used_channels]))
fit_frame$.y <- y_obs
final_fit <- lm(.y ~ ., data = fit_frame)
aliased_channels <- character(0)
cf <- coef(final_fit)
bad <- names(cf)[is.na(cf)]
bad_ch <- names(ch_alias)[ch_alias %in% bad]
if (length(bad_ch) > 0) {
aliased_channels <- unname(channel_names[bad_ch])
used_channels <- setdiff(used_channels, bad_ch)
if (length(used_channels) < 1) {
stop(sprintf("Every mapped spend column is an exact linear combination of the others, so no channel-level effect on '%s' can be separated.",
outcome_name))
}
Zcur <- Zcur[, used_channels, drop = FALSE]
ch_alias <- setNames(paste0("ch_", seq_along(used_channels)), used_channels)
fit_frame <- as.data.frame(cbind(ctrl, Zcur)[obs_idx, , drop = FALSE])
names(fit_frame) <- c(ctrl_names, unname(ch_alias[used_channels]))
fit_frame$.y <- y_obs
final_fit <- lm(.y ~ ., data = fit_frame)
}
sm <- summary(final_fit)
cmat <- sm$coefficients
r2 <- as.numeric(sm$r.squared)
adj_r2 <- as.numeric(sm$adj.r.squared)
dfres <- final_fit$df.residual
tcrit <- if (dfres > 0) qt(0.975, dfres) else NA_real_
resid_v <- as.numeric(residuals(final_fit))
dw <- if (sum(resid_v^2) > 0) sum(diff(resid_v)^2) / sum(resid_v^2) else NA_real_Step 10: Contribution, share, return and marginal return per channel
total_outcome <- sum(y_obs)
ch_rows <- list()
for (i in seq_along(used_channels)) {
sc <- used_channels[i]
ali <- unname(ch_alias[sc])
zi <- Zcur[obs_idx, sc]
sum_z <- sum(zi)
est <- if (ali %in% rownames(cmat)) cmat[ali, "Estimate"] else NA_real_
se <- if (ali %in% rownames(cmat)) cmat[ali, "Std. Error"] else NA_real_
pv <- if (ali %in% rownames(cmat)) cmat[ali, "Pr(>|t|)"] else NA_real_
contrib <- est * sum_z
c_lo <- (est - tcrit * se) * sum_z
c_hi <- (est + tcrit * se) * sum_z
spend_tot <- sum(spend_grid[[sc]][obs_idx])
roi <- if (spend_tot > 0) contrib / spend_tot else NA_real_
r_lo <- if (spend_tot > 0) c_lo / spend_tot else NA_real_
r_hi <- if (spend_tot > 0) c_hi / spend_tot else NA_real_
# Marginal return: lift the channel's whole spend path by 1 percent and
# re-run its own fitted adstock + saturation. Captures carryover and
# curvature exactly rather than approximating the derivative.
bump <- spend_grid[[sc]] * 1.01
z_bump <- saturate_negexp(adstock_geometric(bump, par[[sc]]$theta), par[[sc]]$kappa)
d_spend <- 0.01 * spend_tot
kfac <- if (d_spend > 0) (sum(z_bump[obs_idx]) - sum_z) / d_spend else NA_real_
mroi <- est * kfac
m_lo <- (est - tcrit * se) * kfac
m_hi <- (est + tcrit * se) * kfac
sat90 <- 2.302585 * par[[sc]]$kappa * (1 - par[[sc]]$theta)
avg_spend <- mean(spend_grid[[sc]][obs_idx])
ss_now <- if (par[[sc]]$theta < 1) avg_spend / (1 - par[[sc]]$theta) else avg_spend
pct_ceiling <- 100 * saturate_negexp(ss_now, par[[sc]]$kappa)
ch_rows[[length(ch_rows) + 1]] <- data.frame(
channel = unname(channel_names[[sc]]),
contribution = round(contrib, 1),
share_pct = round(100 * contrib / total_outcome, 2),
spend = round(spend_tot, 1),
return_per_spend = round(roi, 4),
marginal_return = round(mroi, 4),
decay_rate = round(par[[sc]]$theta, 2),
saturation_point = round(sat90, 1),
pct_of_ceiling = round(pct_ceiling, 1),
p_value = signif(pv, 3),
ci_low = round(r_lo, 4),
ci_high = round(r_hi, 4),
m_low = round(m_lo, 4),
m_high = round(m_hi, 4),
contrib_low = round(c_lo, 1),
contrib_high = round(c_hi, 1),
avg_spend = round(avg_spend, 1),
semantic = sc,
stringsAsFactors = FALSE
)
}
channel_all <- do.call(rbind, ch_rows)
channel_all <- channel_all[order(-abs(channel_all$contribution)), , drop = FALSE]
rownames(channel_all) <- NULL
channel_all$significance <- ifelse(is.na(channel_all$p_value), "",
ifelse(channel_all$p_value < 0.001, "p below 0.001",
ifelse(channel_all$p_value < 0.01, "p below 0.01",
ifelse(channel_all$p_value < 0.05, "p below 0.05",
"not distinguishable from zero"))))
marketing_total <- sum(channel_all$contribution, na.rm = TRUE)
marketing_share <- 100 * marketing_total / total_outcome
baseline_share <- 100 - marketing_share
degenerate_decomposition <- isTRUE(marketing_share > 100) || isTRUE(baseline_share < 0)Step 11: Collinearity — can the per-channel split be believed?
base_fit <- rep(coef(final_fit)[["(Intercept)"]], n_obs)
for (nm in ctrl_names) {
cv <- coef(final_fit)[[nm]]
if (!is.null(cv) && is.finite(cv)) base_fit <- base_fit + cv * ctrl[obs_idx, nm]
}
vifs <- setNames(rep(NA_real_, length(used_channels)), used_channels)
if (length(used_channels) >= 2) {
for (sc in used_channels) {
others <- cbind(ctrl[obs_idx, , drop = FALSE],
Zcur[obs_idx, setdiff(used_channels, sc), drop = FALSE])
aux <- tryCatch(summary(lm(Zcur[obs_idx, sc] ~ others)), error = function(e) NULL)
r2j <- if (!is.null(aux)) as.numeric(aux$r.squared) else NA_real_
vifs[sc] <- if (is.na(r2j)) NA_real_ else if (r2j >= 1 - 1e-10) 9999 else 1 / (1 - r2j)
}
} else {
vifs[used_channels] <- 1
}
pair_max <- setNames(rep(NA_real_, length(used_channels)), used_channels)
if (length(used_channels) >= 2) {
S <- sapply(used_channels, function(sc) spend_grid[[sc]][obs_idx])
CM <- suppressWarnings(cor(S))
for (sc in used_channels) {
v <- CM[sc, setdiff(used_channels, sc)]
v <- v[is.finite(v)]
pair_max[sc] <- if (length(v)) max(abs(v)) else NA_real_
}
}
base_corr <- setNames(rep(NA_real_, length(used_channels)), used_channels)
for (sc in used_channels) {
v <- suppressWarnings(cor(spend_grid[[sc]][obs_idx], base_fit))
base_corr[sc] <- if (is.finite(v)) v else NA_real_
}
vif_max <- if (all(is.na(vifs))) NA_real_ else max(vifs, na.rm = TRUE)
collin_verdict <- if (is.na(vif_max)) "not assessable" else
if (vif_max >= 10) "unreliable" else if (vif_max >= 5) "fragile" else "stable"
collin_phrase <- switch(collin_verdict,
unreliable = "the per-channel split is NOT trustworthy",
fragile = "the per-channel split is fragile and should be read as indicative",
stable = "the per-channel split is stable",
"the per-channel split could not be assessed")
curve_span <- setNames(sapply(used_channels, function(sc) {
z <- Zcur[obs_idx, sc]
if (!length(z) || !all(is.finite(z))) NA_real_ else diff(range(z))
}), used_channels)
# At or below the barrier means the saturation scale was decided by the
# identification constraint rather than by the data.
weak_span_channels <- unname(channel_names[used_channels[
is.finite(curve_span) & curve_span <= MIN_TRANSFORM_SPAN + 1e-6]])
collinearity_df <- data.frame(
channel = unname(channel_names[used_channels]),
vif = round(pmin(vifs, 9999), 2),
max_pair_correlation = round(pair_max, 3),
baseline_correlation = round(base_corr, 3),
curve_span = round(curve_span, 3),
verdict = ifelse(is.na(vifs), "not assessable",
ifelse(vifs >= 10, "not separable from the other channels",
ifelse(vifs >= 5, "partly entangled with the other channels",
ifelse(is.finite(curve_span) & curve_span <= MIN_TRANSFORM_SPAN + 1e-6,
"separable, but its contribution level is weakly identified",
"separately identified")))),
stringsAsFactors = FALSE
)
collinearity_df <- collinearity_df[order(-collinearity_df$vif), , drop = FALSE]
rownames(collinearity_df) <- NULL
neg_channels <- channel_all$channel[is.finite(channel_all$contribution) &
channel_all$contribution < 0]
demand_followers <- collinearity_df$channel[is.finite(collinearity_df$baseline_correlation) &
abs(collinearity_df$baseline_correlation) >= 0.5]Step 12: Holdout accuracy at the chosen parameters
holdout_rmse <- NA_real_; holdout_mape <- NA_real_
if (use_holdout) {
Xf <- cbind(`(Intercept)` = 1, ctrl, Zcur)[obs_idx, , drop = FALSE]
ftr <- tryCatch(stats::lm.fit(Xf[train_i, , drop = FALSE], y_obs[train_i]),
error = function(e) NULL)
if (!is.null(ftr)) {
b <- ftr$coefficients; b[is.na(b)] <- 0
pr <- as.numeric(Xf[test_i, , drop = FALSE] %*% b)
if (all(is.finite(pr))) {
holdout_rmse <- sqrt(mean((y_obs[test_i] - pr)^2))
den <- abs(y_obs[test_i])
okd <- den > 0
if (any(okd)) holdout_mape <- 100 * mean(abs(y_obs[test_i][okd] - pr[okd]) / den[okd])
}
}
}Step 13: Decomposition over time (stacked view)
comp_mat <- data.frame(period = format(grid_dates[obs_idx], "%Y-%m-%d"),
stringsAsFactors = FALSE)
comp_mat[["Baseline(no marketing)"]] <- base_fit
for (sc in used_channels) {
est <- cmat[unname(ch_alias[sc]), "Estimate"]
comp_mat[[unname(channel_names[[sc]])]] <- est * Zcur[obs_idx, sc]
}
n_comp <- ncol(comp_mat) - 1L
max_chart_rows <- 1200L
max_periods <- max(8L, as.integer(floor(max_chart_rows / max(1L, n_comp))))
bucket_size <- 1L
if (n_obs > max_periods) {
bucket_size <- as.integer(ceiling(n_obs / max_periods))
grp <- ((seq_len(n_obs) - 1L) %/% bucket_size) + 1L
lab <- tapply(comp_mat$period, grp, function(x) x[1])
num <- comp_mat[, -1, drop = FALSE]
agg2 <- rowsum(as.matrix(num), group = grp)
comp_mat <- data.frame(period = as.character(lab), agg2,
check.names = FALSE, stringsAsFactors = FALSE)
}
decomp_df <- do.call(rbind, lapply(setdiff(names(comp_mat), "period"), function(cn) {
data.frame(period = comp_mat$period, component = cn,
contribution = round(as.numeric(comp_mat[[cn]]), 1),
stringsAsFactors = FALSE)
}))
rownames(decomp_df) <- NULLStep 14: Response curves (sustained-spend steady state)
curve_rows <- list()
for (sc in used_channels) {
est <- cmat[unname(ch_alias[sc]), "Estimate"]
th <- par[[sc]]$theta; kp <- par[[sc]]$kappa
xmax <- max(spend_grid[[sc]][obs_idx])
if (!is.finite(xmax) || xmax <= 0) next
xs <- seq(0, 1.5 * xmax, length.out = 25)
a_ss <- if (th < 1) xs / (1 - th) else xs
curve_rows[[length(curve_rows) + 1]] <- data.frame(
channel = unname(channel_names[[sc]]),
spend_level = round(xs, 1),
response = round(est * saturate_negexp(a_ss, kp), 1),
stringsAsFactors = FALSE)
}
curves_df <- do.call(rbind, curve_rows)
rownames(curves_df) <- NULLStep 16: Metrics + json answer
top_i <- safe_which_max(channel_all$contribution)
top_channel <- if (is.na(top_i)) "not estimable" else channel_all$channel[top_i]
top_share <- if (is.na(top_i)) NA_real_ else channel_all$share_pct[top_i]
roi_i <- safe_which_max(channel_all$return_per_spend)
best_roi_channel <- if (is.na(roi_i)) "not estimable" else channel_all$channel[roi_i]
best_roi <- if (is.na(roi_i)) NA_real_ else channel_all$return_per_spend[roi_i]
metrics <- list(
`Periods Analysed` = n_obs,
`Channels Modelled` = length(used_channels),
`Model R-squared` = round(r2, 3),
`Marketing-Attributed Share` = fmt_pct(marketing_share),
`Baseline Share` = fmt_pct(baseline_share),
`Largest Contributor` = top_channel,
`Largest Contributor Share` = fmt_pct(top_share),
`Best Return per Spend` = best_roi_channel,
`Highest Channel VIF` = if (is.na(vif_max)) NA_real_ else round(min(vif_max, 9999), 2),
`Collinearity Verdict` = collin_verdict,
`Decomposition Verdict` = if (degenerate_decomposition) "degenerate" else "coherent",
`Parameter Selection` = selection_rule
)
neg_sentence <- if (length(neg_channels) > 0) paste0(
" ", paste(neg_channels, collapse = " and "),
if (length(neg_channels) == 1) " carries a NEGATIVE estimated effect" else " carry NEGATIVE estimated effects",
" and is reported that way rather than clipped to zero — in observational spend data that usually means the spend rose when the outcome was falling, not that the advertising destroyed demand.") else ""
collin_sentence <- if (collin_verdict == "unreliable") sprintf(
" Channel spends move together too closely(highest variance inflation factor %s), so %s: read the combined marketing number, not the per-channel split.",
formatC(min(vif_max, 9999), format = "f", digits = 1), collin_phrase)
else if (collin_verdict == "fragile") sprintf(
" Channel spends partly move together(highest variance inflation factor %s), so %s.",
formatC(vif_max, format = "f", digits = 1), collin_phrase)
else sprintf(" Channel spends vary independently enough(highest variance inflation factor %s) that %s.",
if (is.na(vif_max)) "not assessable" else formatC(vif_max, format = "f", digits = 1),
collin_phrase)
degenerate_sentence <- if (degenerate_decomposition) sprintf(
" WARNING: this decomposition is degenerate — it credits marketing with %s of total %s and leaves the non-marketing baseline at %s. No model can honestly claim marketing produced more than the whole outcome, so the split below should not be used; look for an outlier period, a control the model is missing, or a channel whose spend barely varies.",
fmt_pct(marketing_share), outcome_name, fmt_pct(baseline_share)) else ""
span_sentence <- if (length(weak_span_channels) > 0) sprintf(
" %s never move far along their own response curve in this data(the spend range covers less than half of it), so their contribution LEVEL leans on the model's assumed curve shape rather than on anything observed; treat those totals as the softest numbers in the report.",
paste(weak_span_channels, collapse = " and ")) else ""
demand_sentence <- if (length(demand_followers) > 0) sprintf(
" %s spend tracks the fitted baseline closely(correlation of %s or more), so %s credited with demand that would have arrived anyway.",
paste(demand_followers, collapse = " and "),
formatC(max(abs(base_corr[is.finite(base_corr)])), format = "f", digits = 2),
if (length(demand_followers) == 1) "part of its contribution is likely" else "part of their contribution is likely")
else " No channel's spend tracks the fitted baseline closely enough to suggest it is simply following demand."
json_output <- list(
answer = paste0(
"Media mix model of ", outcome_name, " over ", fmt_num(n_obs), " ", cadence,
" periods with ", length(used_channels), " channel(s): the model reproduces ",
fmt_pct(100 * r2), " of the movement, and attributes ", fmt_pct(marketing_share),
" of total ", outcome_name, " to marketing and ", fmt_pct(baseline_share),
" to the non-marketing baseline, trend and seasonality. ",
"The largest contributor is ", top_channel, " at ", fmt_pct(top_share),
" of the outcome; the best return per unit of spend is ", best_roi_channel,
if (is.na(best_roi)) "" else paste0(" at ", formatC(best_roi, format = "f", digits = 3),
" units of ", outcome_name, " per unit of spend"),
". Carryover and saturation were fitted per channel by grid search on ",
selection_label, ".", degenerate_sentence, collin_sentence, span_sentence,
demand_sentence, neg_sentence,
" These are associations fitted to observed spend, not experimental effects: the honest test of the causal claim is a holdout or geo experiment, which this analysis is not."
),
cards = lapply(
c("tldr", "overview", "preprocessing", "decomposition", "channel_contribution",
"roi_intervals", "response_curves", "collinearity", "model_fit", "methods"),
function(cid) list(id = cid, metrics = metrics)
)
)
list(
initial_rows = initial_rows, final_rows = n_obs,
rows_removed = max(0L, initial_rows - n_obs),
date_name = date_name, outcome_name = outcome_name,
channel_names = channel_names, used_channels = used_channels,
dropped_df = dropped_df, aliased_channels = aliased_channels,
cadence = cadence, period_word = period_word, n_periods = n_obs,
n_grid = n_grid, n_missing_periods = n_missing_periods,
n_dup_rows = n_dup_rows, n_blank_spend = n_blank_spend,
n_missing_outcome = n_missing_outcome,
par = par, theta_grid = theta_grid, n_evals = n_evals, n_passes = n_passes,
selection_rule = selection_rule, selection_label = selection_label,
holdout_h = h, holdout_rmse = holdout_rmse, holdout_mape = holdout_mape,
channel_all = channel_all, channel_df = channel_df, roi_df = roi_df,
decomp_df = decomp_df, curves_df = curves_df,
collinearity_df = collinearity_df, fit_df_out = fit_df_out,
methods_df = methods_df, bucket_size = bucket_size,
seasonal_terms = seasonal_terms, seasonal_note = seasonal_note,
n_params = n_params, r2 = r2, adj_r2 = adj_r2, dw = dw,
total_outcome = total_outcome, marketing_total = marketing_total,
marketing_share = marketing_share, baseline_share = baseline_share,
vifs = vifs, vif_max = vif_max, collin_verdict = collin_verdict,
curve_span = curve_span, weak_span_channels = weak_span_channels,
degenerate_decomposition = degenerate_decomposition,
degenerate_sentence = degenerate_sentence,
span_sentence = span_sentence, n_rejected_scales = n_rejected_scales,
min_transform_span = MIN_TRANSFORM_SPAN,
collin_phrase = collin_phrase, base_corr = base_corr,
neg_channels = neg_channels, demand_followers = demand_followers,
top_channel = top_channel, top_share = top_share,
best_roi_channel = best_roi_channel, best_roi = best_roi,
collin_sentence = collin_sentence, demand_sentence = demand_sentence,
neg_sentence = neg_sentence,
metrics = metrics, json_output = json_output
)
}