Executive Summary

This analysis implements a comprehensive Marketing Mix Model (MMM) for pharmaceutical marketing channels using:

Panel regression with fixed/random effects
Adstock transformation to capture carryover effects
ROI and marginal ROI calculations from regression coefficients
Two optimization scenarios: Fixed budget reallocation and unconstrained budget expansion

1. Setup and Package Loading

# Install and load required packages
required_packages <- c("dplyr", "tidyr", "scales", "gt", "plotly", "plm", "lmtest")

suppressPackageStartupMessages({
  for(pkg in required_packages) {
    if(!require(pkg, character.only = TRUE, quietly = TRUE)) {
      install.packages(pkg, quiet = TRUE)
      library(pkg, character.only = TRUE, quietly = TRUE)
    }
  }
})

# Verify plm is loaded
if(!"plm" %in% loadedNamespaces()) {
  stop("Package 'plm' failed to load. Please install manually: install.packages('plm')")
}

set.seed(123)

2. Data Generation

2.1 Generate Panel Data

generate_mmm_data <- function(n_periods = 52, n_regions = 20) {
  # Create panel structure: regions over time (increased to 20 regions)
  regions <- paste0("Region_", LETTERS[1:n_regions])
  
  data <- expand.grid(
    Week = 1:n_periods,
    Region = regions
  )
  
  # Get number of rows for random generation
  n_rows <- nrow(data)
  
  # Generate marketing spend by channel
  data <- data %>%
    mutate(
      # Regional fixed effects
      region_effect = as.numeric(factor(Region)) * 50000,
      
      # Marketing channels with realistic spending patterns
      Digital_Spend = rnorm(n_rows, 80000, 15000) + region_effect * 0.3,
      TV_Spend = rnorm(n_rows, 150000, 25000) + region_effect * 0.5,
      Print_Spend = rnorm(n_rows, 40000, 8000) + region_effect * 0.2,
      Sales_Rep_Spend = rnorm(n_rows, 120000, 20000) + region_effect * 0.4,
      Events_Spend = rnorm(n_rows, 60000, 12000) + region_effect * 0.25,
      
      # Control variables
      Seasonality = sin(2 * pi * Week / 52) * 100000 + 500000,
      Competitor_Spend = rnorm(n_rows, 200000, 30000),
      Price = rnorm(n_rows, 100, 10)
    ) %>%
    # Ensure non-negative spending
    mutate(across(ends_with("_Spend"), ~pmax(., 0)))
  
  return(data)
}

# Generate base data
mmm_data <- generate_mmm_data(n_periods = 52, n_regions = 20)

cat("Data dimensions:", nrow(mmm_data), "rows x", ncol(mmm_data), "columns\n")

## Data dimensions: 1040 rows x 11 columns

cat("Time periods:", length(unique(mmm_data$Week)), "weeks\n")

## Time periods: 52 weeks

cat("Regions:", length(unique(mmm_data$Region)), "regions\n")

## Regions: 20 regions

2.2 Adstock Transformation

# Geometric adstock transformation
adstock_geometric <- function(x, decay_rate = 0.5) {
  adstocked <- numeric(length(x))
  adstocked[1] <- x[1]
  
  for(i in 2:length(x)) {
    adstocked[i] <- x[i] + decay_rate * adstocked[i-1]
  }
  
  return(adstocked)
}

# Apply adstock to panel data by group
apply_adstock_panel <- function(data, channels, decay_rates) {
  
  for(i in seq_along(channels)) {
    channel <- channels[i]
    decay <- decay_rates[i]
    adstock_col_name <- paste0(channel, "_Adstock")
    
    # Split by region, apply adstock, combine
    data_split <- split(data, data$Region)
    
    for(region in names(data_split)) {
      region_data <- data_split[[region]]
      region_data <- region_data[order(region_data$Week), ]
      
      # Apply adstock to this region's data
      adstock_values <- adstock_geometric(region_data[[channel]], decay)
      region_data[[adstock_col_name]] <- adstock_values
      
      data_split[[region]] <- region_data
    }
    
    # Combine back
    data <- do.call(rbind, data_split)
    rownames(data) <- NULL
  }
  
  return(data)
}

# Define channels and decay rates
channels <- c("Digital_Spend", "TV_Spend", "Print_Spend", 
              "Sales_Rep_Spend", "Events_Spend")

# Decay rates: Digital=0.3 (fast decay), TV=0.7 (long memory), 
# Print=0.5, Sales Rep=0.6 (personal touch lasts), Events=0.4
decay_rates <- c(0.3, 0.7, 0.5, 0.6, 0.4)

# Apply adstock transformation
mmm_data <- apply_adstock_panel(mmm_data, channels, decay_rates)

# Verify adstock columns were created
cat("\nAdstock columns created:\n")

## 
## Adstock columns created:

adstock_cols <- grep("_Adstock$", names(mmm_data), value = TRUE)
print(adstock_cols)

## [1] "Digital_Spend_Adstock"   "TV_Spend_Adstock"       
## [3] "Print_Spend_Adstock"     "Sales_Rep_Spend_Adstock"
## [5] "Events_Spend_Adstock"

2.3 Generate Revenue

# Generate revenue based on adstocked variables with realistic effects
mmm_data <- mmm_data %>%
  mutate(
    # True data generating process
    Revenue = 200000 + 
      # Adstocked channel effects (coefficients represent true ROI components)
      0.85 * Digital_Spend_Adstock +
      0.65 * TV_Spend_Adstock +
      1.10 * Print_Spend_Adstock +
      0.75 * Sales_Rep_Spend_Adstock +
      0.55 * Events_Spend_Adstock +
      # Control variables
      0.8 * Seasonality +
      -0.3 * Competitor_Spend +
      2000 * Price +
      # Regional fixed effects
      region_effect +
      # Random error
      rnorm(nrow(.), 0, 50000)
  ) %>%
  mutate(Revenue = pmax(Revenue, 0))

cat("Revenue summary:\n")

## Revenue summary:

summary(mmm_data$Revenue)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 1258001 2548738 3472784 3492740 4458112 5601914

3. Panel Regression Models

cat("\n========================================")

## 
## ========================================

cat("\n MARKETING MIX MODEL ESTIMATION")

## 
##  MARKETING MIX MODEL ESTIMATION

cat("\n========================================\n")

## 
## ========================================

# Prepare panel data object
mmm_data$Region <- as.character(mmm_data$Region)
mmm_data$Week <- as.integer(mmm_data$Week)

panel_data <- plm::pdata.frame(mmm_data, index = c("Region", "Week"))

# Model formula
formula_pooled <- Revenue ~ Digital_Spend_Adstock + TV_Spend_Adstock + 
  Print_Spend_Adstock + Sales_Rep_Spend_Adstock + Events_Spend_Adstock +
  Seasonality + Competitor_Spend + Price

# Model 1: Pooled OLS (baseline)
model_pooled <- plm::plm(formula_pooled, data = panel_data, model = "pooling")

# Model 2: Fixed Effects Model
model_fe <- plm::plm(formula_pooled, data = panel_data, model = "within")

# Model 3: Random Effects Model
model_re <- tryCatch({
  plm::plm(formula_pooled, data = panel_data, model = "random", random.method = "swar")
}, error = function(e) {
  message("Swamy-Arora method failed, using alternative random effects estimation")
  plm::plm(formula_pooled, data = panel_data, model = "random", random.method = "amemiya")
})

# Model 4: Mixed Effects Model
cat("\n--- Estimating Mixed Effects Model ---\n")

## 
## --- Estimating Mixed Effects Model ---

model_mixed <- tryCatch({
  plm::pvcm(Revenue ~ Digital_Spend_Adstock + TV_Spend_Adstock + 
              Print_Spend_Adstock + Sales_Rep_Spend_Adstock + Events_Spend_Adstock +
              Seasonality + Competitor_Spend + Price,
            data = panel_data, 
            model = "random",
            effect = "individual")
}, error = function(e) {
  message("Mixed effects model failed, will use fixed or random effects instead")
  NULL
})

# Hausman test
hausman_test <- tryCatch({
  plm::phtest(model_fe, model_re)
}, error = function(e) {
  message("Hausman test could not be computed, defaulting to Fixed Effects model")
  list(p.value = 0.01, statistic = NA, method = "Manual selection")
})

cat("\n--- Hausman Test (Fixed vs Random Effects) ---\n")

## 
## --- Hausman Test (Fixed vs Random Effects) ---

print(hausman_test)

## 
##  Hausman Test
## 
## data:  formula_pooled
## chisq = 632.7, df = 8, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

# Select best model
if(hausman_test$p.value < 0.05) {
  best_model <- model_fe
  model_type <- "Fixed Effects"
  cat("\n*** Fixed Effects Model Selected (p < 0.05) ***\n")
} else {
  best_model <- model_re
  model_type <- "Random Effects"
  cat("\n*** Random Effects Model Selected (p >= 0.05) ***\n")
}

## 
## *** Fixed Effects Model Selected (p < 0.05) ***

3.1 Model Summary Table

# Extract model summary
model_summary_data <- summary(best_model)

# Extract coefficients table
coef_table <- as.data.frame(coef(model_summary_data))
coef_table$Variable <- rownames(coef_table)
rownames(coef_table) <- NULL

coef_table <- coef_table %>%
  select(Variable, Estimate, `Std. Error`, `t-value`, `Pr(>|t|)`)

# Create gt table
gt_model_summary <- coef_table %>%
  gt() %>%
  tab_header(
    title = paste(model_type, "Model: Regression Coefficients"),
    subtitle = "Marketing Mix Model with Adstock Transformation"
  ) %>%
  fmt_number(
    columns = c(Estimate, `Std. Error`, `t-value`, `Pr(>|t|)`),
    decimals = 4
  ) %>%
  cols_label(
    Variable = "Variable",
    Estimate = "Coefficient",
    `Std. Error` = "Std. Error",
    `t-value` = "t-statistic",
    `Pr(>|t|)` = "p-value"
  ) %>%
  data_color(
    columns = `Pr(>|t|)`,
    colors = scales::col_numeric(
      palette = c("#006400", "#90EE90", "#FFD700", "#FF6B6B"),
      domain = c(0, 0.1),
      reverse = TRUE
    )
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = everything(), rows = `Pr(>|t|)` < 0.05)
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  ) %>%
  tab_source_note(
    source_note = paste0(
      "Model: ", model_type, " | R² = ", sprintf("%.4f", model_summary_data$r.squared[1])
    )
  )

gt_model_summary

Variable	Coefficient	Std. Error	t-statistic	p-value
Fixed Effects Model: Regression Coefficients
Marketing Mix Model with Adstock Transformation
Digital_Spend_Adstock	0.8621	0.0993	8.6785	0.0000
TV_Spend_Adstock	0.6502	0.0264	24.5867	0.0000
Print_Spend_Adstock	1.0468	0.1512	6.9216	0.0000
Sales_Rep_Spend_Adstock	0.7564	0.0572	13.2171	0.0000
Events_Spend_Adstock	0.5430	0.1168	4.6495	0.0000
Seasonality	0.8165	0.0231	35.2858	0.0000
Competitor_Spend	−0.3000	0.0528	−5.6831	0.0000
Price	1,911.4950	168.3908	11.3515	0.0000
Model: Fixed Effects \| R² = 0.9560

3.2 Model Comparison

model_comparison <- data.frame(
  Model = c("Pooled OLS", "Fixed Effects", "Random Effects", 
            if(!is.null(model_mixed)) "Mixed Effects" else NULL, "Selected Model"),
  R_Squared = c(
    summary(model_pooled)$r.squared[1],
    summary(model_fe)$r.squared[1],
    summary(model_re)$r.squared[1],
    if(!is.null(model_mixed)) summary(model_mixed)$r.squared[1] else NULL,
    summary(best_model)$r.squared[1]
  ),
  Adj_R_Squared = c(
    summary(model_pooled)$r.squared[2],
    summary(model_fe)$r.squared[2],
    summary(model_re)$r.squared[2],
    if(!is.null(model_mixed)) summary(model_mixed)$r.squared[2] else NULL,
    summary(best_model)$r.squared[2]
  ),
  N_Obs = c(
    nobs(model_pooled), nobs(model_fe), nobs(model_re),
    if(!is.null(model_mixed)) nobs(model_mixed) else NULL, nobs(best_model)
  )
)

gt_models <- model_comparison %>%
  gt() %>%
  tab_header(title = "Panel Regression Model Comparison") %>%
  fmt_number(columns = c(R_Squared, Adj_R_Squared), decimals = 4) %>%
  fmt_number(columns = N_Obs, decimals = 0) %>%
  tab_style(
    style = cell_fill(color = "#E8F4F8"),
    locations = cells_body(rows = Model == "Selected Model")
  )

gt_models

Model	R_Squared	Adj_R_Squared	N_Obs
Panel Regression Model Comparison
Pooled OLS	0.9966	0.9966	1,040
Fixed Effects	0.9560	0.9548	1,040
Random Effects	0.9956	0.9956	1,040
Selected Model	0.9560	0.9548	1,040

4. ROI and mROI Calculations

# Extract coefficients
coefficients <- coef(best_model)
channel_coefs <- coefficients[grep("_Spend_Adstock", names(coefficients))]

# Calculate totals
spend_summary <- mmm_data %>%
  summarise(
    Total_Digital_Spend = sum(Digital_Spend),
    Total_TV_Spend = sum(TV_Spend),
    Total_Print_Spend = sum(Print_Spend),
    Total_Sales_Rep_Spend = sum(Sales_Rep_Spend),
    Total_Events_Spend = sum(Events_Spend),
    Total_Digital_Adstock = sum(Digital_Spend_Adstock),
    Total_TV_Adstock = sum(TV_Spend_Adstock),
    Total_Print_Adstock = sum(Print_Spend_Adstock),
    Total_Sales_Rep_Adstock = sum(Sales_Rep_Spend_Adstock),
    Total_Events_Adstock = sum(Events_Spend_Adstock)
  )

# Create ROI dataframe
roi_data <- data.frame(
  Channel = c("Digital", "TV", "Print", "Sales_Rep", "Events"),
  Coefficient = c(
    channel_coefs["Digital_Spend_Adstock"],
    channel_coefs["TV_Spend_Adstock"],
    channel_coefs["Print_Spend_Adstock"],
    channel_coefs["Sales_Rep_Spend_Adstock"],
    channel_coefs["Events_Spend_Adstock"]
  ),
  Total_Spend = c(
    spend_summary$Total_Digital_Spend,
    spend_summary$Total_TV_Spend,
    spend_summary$Total_Print_Spend,
    spend_summary$Total_Sales_Rep_Spend,
    spend_summary$Total_Events_Spend
  ),
  Total_Adstock = c(
    spend_summary$Total_Digital_Adstock,
    spend_summary$Total_TV_Adstock,
    spend_summary$Total_Print_Adstock,
    spend_summary$Total_Sales_Rep_Adstock,
    spend_summary$Total_Events_Adstock
  ),
  Decay_Rate = decay_rates
)

rownames(roi_data) <- NULL

roi_data <- roi_data %>%
  mutate(
    Total_Contribution = Coefficient * Total_Adstock,
    ROI = Total_Contribution / Total_Spend,
    mROI_Immediate = Coefficient,
    mROI_Steady_State = Coefficient / (1 - Decay_Rate),
    Efficiency = ROI * (1 + Decay_Rate)
  ) %>%
  arrange(desc(ROI))

4.1 ROI and mROI Table

gt_roi <- roi_data %>%
  select(Channel, Total_Spend, Total_Contribution, ROI, 
         mROI_Immediate, mROI_Steady_State, Decay_Rate) %>%
  gt() %>%
  tab_header(
    title = "Marketing Mix Model: ROI & Marginal ROI Analysis",
    subtitle = "Based on Panel Regression with Adstock Transformation"
  ) %>%
  fmt_currency(columns = c(Total_Spend, Total_Contribution), decimals = 0) %>%
  fmt_number(columns = c(ROI, mROI_Immediate, mROI_Steady_State), decimals = 3) %>%
  fmt_number(columns = Decay_Rate, decimals = 2) %>%
  data_color(
    columns = ROI,
    colors = scales::col_numeric(
      palette = c("#FF6B6B", "#FFD93D", "#6BCF7F", "#4D96FF"),
      domain = NULL
    )
  )

gt_roi

Channel	Total_Spend	Total_Contribution	ROI	mROI_Immediate	mROI_Steady_State	Decay_Rate
Marketing Mix Model: ROI & Marginal ROI Analysis
Based on Panel Regression with Adstock Transformation
TV	$429,749,589	$889,599,331	2.070	0.650	2.167	0.70
Print	$150,589,946	$309,150,484	2.053	1.047	2.094	0.50
Sales_Rep	$342,961,943	$629,862,865	1.837	0.756	1.891	0.60
Digital	$247,390,231	$302,196,735	1.222	0.862	1.232	0.30
Events	$198,776,067	$177,577,485	0.893	0.543	0.905	0.40

5. Budget Optimization

5.1 Scenario 1: Fixed Budget

total_budget <- sum(roi_data$Total_Spend)

optimal_allocation_fixed <- roi_data %>%
  mutate(
    Weight = mROI_Steady_State / sum(mROI_Steady_State),
    Optimal_Budget = Weight * total_budget,
    Current_Budget = Total_Spend,
    Budget_Change = Optimal_Budget - Current_Budget,
    Change_Pct = (Budget_Change / Current_Budget) * 100
  ) %>%
  select(Channel, Current_Budget, Optimal_Budget, Budget_Change, 
         Change_Pct, ROI, mROI_Steady_State)

gt_allocation_fixed <- optimal_allocation_fixed %>%
  gt() %>%
  tab_header(
    title = "SCENARIO 1: Fixed Budget Allocation",
    subtitle = "Budget reallocation within current total spend"
  ) %>%
  fmt_currency(columns = c(Current_Budget, Optimal_Budget, Budget_Change), decimals = 0) %>%
  fmt_number(columns = c(Change_Pct, ROI, mROI_Steady_State), decimals = 2) %>%
  data_color(
    columns = Change_Pct,
    colors = scales::col_numeric(
      palette = c("#FF6B6B", "#FFFFFF", "#6BCF7F"),
      domain = c(-100, 100)
    )
  )

gt_allocation_fixed

Channel	Current_Budget	Optimal_Budget	Budget_Change	Change_Pct	ROI	mROI_Steady_State
SCENARIO 1: Fixed Budget Allocation
Budget reallocation within current total spend
TV	$429,749,589	$358,096,024	−$71,653,565	−16.67	2.07	2.17
Print	$150,589,946	$345,916,033	$195,326,087	129.71	2.05	2.09
Sales_Rep	$342,961,943	$312,446,323	−$30,515,620	−8.90	1.84	1.89
Digital	$247,390,231	$203,490,478	−$43,899,753	−17.75	1.22	1.23
Events	$198,776,067	$149,518,918	−$49,257,149	−24.78	0.89	0.90

5.2 Scenario 2: Unconstrained Budget

# Function to calculate optimal spend
calculate_optimal_spend <- function(current_spend, current_contribution, elasticity = 0.7, target_mroi = 1.0) {
  current_mroi <- current_contribution / current_spend
  
  if(elasticity == 1) {
    if(current_mroi > target_mroi) {
      return(current_spend * 10)
    } else {
      return(current_spend)
    }
  }
  
  spend_ratio <- (target_mroi / (elasticity * current_mroi))^(1/(elasticity - 1))
  optimal_spend <- current_spend * spend_ratio
  optimal_spend <- min(optimal_spend, current_spend * 5)
  
  return(optimal_spend)
}

# Calculate unconstrained optimal budget
optimal_allocation_unconstrained <- roi_data %>%
  mutate(
    Current_Budget = Total_Spend,
    Current_mROI = Total_Contribution / Total_Spend,
    Elasticity = 0.7,
    Optimal_Budget_Unconstrained = mapply(
      calculate_optimal_spend,
      current_spend = Total_Spend,
      current_contribution = Total_Contribution,
      elasticity = Elasticity,
      target_mroi = 1.0
    ),
    Spend_Ratio = Optimal_Budget_Unconstrained / Total_Spend,
    Optimal_Contribution = Total_Contribution * (Spend_Ratio^Elasticity),
    Budget_Change = Optimal_Budget_Unconstrained - Current_Budget,
    Budget_Change_Pct = (Budget_Change / Current_Budget) * 100,
    Contribution_Change = Optimal_Contribution - Total_Contribution,
    Incremental_ROI = Contribution_Change / Budget_Change
  ) %>%
  select(Channel, Current_Budget, Current_mROI, Optimal_Budget_Unconstrained,
         Budget_Change, Budget_Change_Pct, Contribution_Change, Incremental_ROI)

gt_allocation_unconstrained <- optimal_allocation_unconstrained %>%
  gt() %>%
  tab_header(
    title = "SCENARIO 2: Unconstrained Budget Allocation",
    subtitle = "Expand spending until marginal ROI equals 1.0"
  ) %>%
  fmt_currency(columns = c(Current_Budget, Optimal_Budget_Unconstrained, 
                            Budget_Change, Contribution_Change), decimals = 0) %>%
  fmt_number(columns = c(Current_mROI, Budget_Change_Pct, Incremental_ROI), decimals = 3)

gt_allocation_unconstrained

Channel	Current_Budget	Current_mROI	Optimal_Budget_Unconstrained	Budget_Change	Budget_Change_Pct	Contribution_Change	Incremental_ROI
SCENARIO 2: Unconstrained Budget Allocation
Expand spending until marginal ROI equals 1.0
TV	$429,749,589	2.070	$1,479,581,400	$1,049,831,812	244.289	$1,224,088,383	1.166
Print	$150,589,946	2.053	$504,315,853	$353,725,907	234.893	$411,300,735	1.163
Sales_Rep	$342,961,943	1.837	$792,330,379	$449,368,435	131.026	$502,037,676	1.117
Digital	$247,390,231	1.222	$146,802,489	−$100,587,742	−40.660	−$92,478,893	0.919
Events	$198,776,067	0.893	$41,569,064	−$157,207,002	−79.087	−$118,193,107	0.752

# Store for compatibility
optimal_allocation <- optimal_allocation_fixed

6. Visualizations

# Add optimal budget for plotting
roi_data_plot <- roi_data %>%
  left_join(optimal_allocation_fixed %>% select(Channel, Optimal_Budget), by = "Channel")

6.1 ROI vs mROI Bubble Chart

p1 <- plot_ly(roi_data_plot, 
              x = ~ROI, y = ~mROI_Steady_State, text = ~Channel, 
              type = 'scatter', mode = 'markers+text',
              marker = list(size = ~sqrt(Total_Spend)/100, color = ~Decay_Rate,
                           colorscale = 'Viridis', showscale = TRUE,
                           colorbar = list(title = "Decay Rate"),
                           line = list(color = 'white', width = 2)),
              textposition = 'top center',
              textfont = list(size = 12, color = 'black'),
              hoverinfo = 'text',
              hovertext = ~paste0('<b>', Channel, '</b><br>',
                                 'ROI: ', round(ROI, 3), '<br>',
                                 'mROI (SS): ', round(mROI_Steady_State, 3), '<br>',
                                 'Current Spend: $', formatC(Total_Spend, format = "f", digits = 0, big.mark = ","), '<br>',
                                 'Optimal Budget: $', formatC(Optimal_Budget, format = "f", digits = 0, big.mark = ","))) %>%
  layout(title = list(text = "ROI vs Marginal ROI (Steady State)"),
         xaxis = list(title = "ROI", zeroline = TRUE),
         yaxis = list(title = "mROI (Steady State)", zeroline = TRUE),
         plot_bgcolor = 'white', paper_bgcolor = 'white')

p1

6.2 Channel Contributions

p2 <- plot_ly(roi_data %>% arrange(desc(Total_Contribution)), 
              x = ~Total_Contribution, y = ~reorder(Channel, Total_Contribution),
              type = 'bar',
              marker = list(color = ~ROI,
                           colorscale = list(c(0, "#FF6B6B"), c(1, "#6BCF7F")),
                           showscale = TRUE)) %>%
  layout(title = "Total Revenue Contribution by Channel",
         xaxis = list(title = "Revenue Contribution ($)"),
         yaxis = list(title = ""))

p2

6.3 Budget Allocation Comparison

p4 <- plot_ly(optimal_allocation_fixed, y = ~Channel) %>%
  add_trace(x = ~Current_Budget, name = 'Current', type = 'bar', orientation = 'h',
            marker = list(color = '#FFA07A')) %>%
  add_trace(x = ~Optimal_Budget, name = 'Optimal', type = 'bar', orientation = 'h',
            marker = list(color = '#6BCF7F')) %>%
  layout(title = "SCENARIO 1: Budget Reallocation",
         xaxis = list(title = "Budget ($)"),
         barmode = 'group')

p4

6.4 Unconstrained Budget Expansion

p6 <- plot_ly(optimal_allocation_unconstrained, y = ~Channel) %>%
  add_trace(x = ~Current_Budget, name = 'Current', type = 'bar', orientation = 'h',
            marker = list(color = '#FFA07A')) %>%
  add_trace(x = ~Optimal_Budget_Unconstrained, name = 'Optimal (Unconstrained)',
            type = 'bar', orientation = 'h', marker = list(color = '#4D96FF')) %>%
  layout(title = "SCENARIO 2: Budget Expansion",
         xaxis = list(title = "Budget ($)"),
         barmode = 'group')

p6

7. Key Insights

best_roi_channel <- roi_data$Channel[which.max(roi_data$ROI)]
best_mroi_channel <- roi_data$Channel[which.max(roi_data$mROI_Steady_State)]

cat("\n========================================")

## 
## ========================================

cat("\n KEY INSIGHTS")

## 
##  KEY INSIGHTS

cat("\n========================================\n")

## 
## ========================================

cat(sprintf("\n1. Model: %s (R² = %.4f)\n", model_type, summary(best_model)$r.squared[1]))

## 
## 1. Model: Fixed Effects (R² = 0.9560)

cat(sprintf("\n2. Best ROI Channel: %s (%.3f)\n", best_roi_channel, max(roi_data$ROI)))

## 
## 2. Best ROI Channel: TV (2.070)

cat(sprintf("\n3. Highest mROI: %s (%.3f)\n", best_mroi_channel, max(roi_data$mROI_Steady_State)))

## 
## 3. Highest mROI: TV (2.167)

cat(sprintf("\n4. Total Budget: $%s\n", formatC(total_budget, format = "f", digits = 0, big.mark = ",")))

## 
## 4. Total Budget: $1,369,467,776

cat(sprintf("\n5. Total Contribution: $%s\n", formatC(sum(roi_data$Total_Contribution), format = "f", digits = 0, big.mark = ",")))

## 
## 5. Total Contribution: $2,308,386,900

cat(sprintf("\n6. Portfolio ROI: %.3f\n", sum(roi_data$Total_Contribution) / total_budget))

## 
## 6. Portfolio ROI: 1.686

cat("\n========================================\n")

## 
## ========================================

Session Info

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] lmtest_0.9-40 zoo_1.8-14    plm_2.6-7     plotly_4.11.0 ggplot2_4.0.0
## [6] gt_1.1.0      scales_1.4.0  tidyr_1.3.1   dplyr_1.1.4  
## 
## loaded via a namespace (and not attached):
##  [1] sandwich_3.1-1     sass_0.4.10        generics_0.1.4     xml2_1.4.0        
##  [5] lattice_0.22-5     digest_0.6.37      magrittr_2.0.4     evaluate_1.0.5    
##  [9] grid_4.5.1         RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0    
## [13] Formula_1.2-5      httr_1.4.7         purrr_1.1.0        crosstalk_1.2.2   
## [17] viridisLite_0.4.2  lazyeval_0.2.2     jquerylib_0.1.4    Rdpack_2.6.4      
## [21] cli_3.6.5          rlang_1.1.6        rbibutils_2.4      miscTools_0.6-28  
## [25] withr_3.0.2        cachem_1.1.0       yaml_2.3.10        parallel_4.5.1    
## [29] tools_4.5.1        bdsmatrix_1.3-7    maxLik_1.5-2.1     vctrs_0.6.5       
## [33] R6_2.6.1           lifecycle_1.0.4    fs_1.6.6           htmlwidgets_1.6.4 
## [37] MASS_7.3-65        pkgconfig_2.0.3    pillar_1.11.1      bslib_0.9.0       
## [41] gtable_0.3.6       Rcpp_1.1.0         glue_1.8.0         data.table_1.17.8 
## [45] collapse_2.1.5     xfun_0.53          tibble_3.3.0       tidyselect_1.2.1  
## [49] rstudioapi_0.17.1  knitr_1.50         farver_2.1.2       nlme_3.1-168      
## [53] htmltools_0.5.8.1  rmarkdown_2.30     compiler_4.5.1     S7_0.2.0

Marketing Mix Modeling with Adstock & Panel Regression

David Guy PhD: DataBridge Statistical Consulting Inc.

2025-12-12

Executive Summary

1. Setup and Package Loading

2. Data Generation

2.1 Generate Panel Data

2.2 Adstock Transformation

2.3 Generate Revenue

3. Panel Regression Models

3.1 Model Summary Table

3.2 Model Comparison

4. ROI and mROI Calculations

4.1 ROI and mROI Table

5. Budget Optimization

5.1 Scenario 1: Fixed Budget

5.2 Scenario 2: Unconstrained Budget

6. Visualizations

6.1 ROI vs mROI Bubble Chart

6.2 Channel Contributions

6.3 Budget Allocation Comparison

6.4 Unconstrained Budget Expansion

7. Key Insights

Session Info