Awesome Marketing Science Guide

# Install# pip install lightweight-mmm
from lightweight_mmm import lightweight_mmmfrom lightweight_mmm import preprocessingimport jax.numpy as jnp
# Prepare datamedia_data = jnp.array(media_data_raw)  # Shape: (n_time_periods, n_media_channels)target = jnp.array(revenue_data)  # Shape: (n_time_periods,)
# Scale datamedia_data_scaled, scaler = preprocessing.CustomScaler().fit_transform(media_data)
# Initialize and train modelmmm = lightweight_mmm.LightweightMMM()mmm.fit(    media=media_data_scaled,    target=target,    number_warmup=1000,    number_samples=1000,    media_prior=jnp.ones(n_channels))
# Get channel contributionscontribution = mmm.get_posterior_metrics()["contribution_per_channel"]

# Install# pip install pymc-marketing
import pymc as pmfrom pymc_marketing.mmm import MMM
# Initialize modelmmm = MMM(    date_column="date",    channel_columns=["tv", "radio", "digital"],    adstock_max_lag=8,    validate_data=True)
# Fit modelmmm.fit(    data=marketing_df,    target_column="sales",    tune=1000,    draws=1000)
# Optimize budgetoptimized = mmm.optimize_budget(    total_budget=100000,    budget_bounds={"tv": (10000, 50000), "radio": (5000, 30000)})
# Get ROI curvesroi_curves = mmm.compute_channel_roi()

# Install# pip install matched-markets
from matched_markets.methodology import trimmed_match
# Prepare datageo_data = {    'geos': ['geo1', 'geo2', 'geo3', ...],    'response': [100, 150, 120, ...],    'spend': [50, 75, 60, ...]}
# Design testdesign = trimmed_match.TrimmedMatch(    data=geo_data,    treatment_geos=['geo1', 'geo2'],    exclude_cooldown_period=7)
# Run matchingdesign.fit()
# Estimate treatment effectiroas = design.estimate_treatment_effect(    post_treatment_response=post_data['response'],    post_treatment_spend=post_data['spend'])
print(f"Incremental ROAS: {iroas.point_estimate:.2f}")print(f"95% CI: [{iroas.ci_lower:.2f}, {iroas.ci_upper:.2f}]")

# Install# pip install pychattr
from pychattr.channel_attribution import MarkovModel
# Journey data format: each row is a conversion path# format: 'channel1 > channel2 > channel3'journeys = [    'paid_search > display > direct',    'social > email > direct',    'direct',    'paid_search > direct']
conversions = [1, 1, 0, 1]revenues = [100, 150, 0, 80]
# Fit Markov modelmarkov = MarkovModel(    paths=journeys,    conversions=conversions,    revenues=revenues,    order=1  # First-order Markov chain)
markov.fit()
# Get attribution resultsattribution = markov.attribution()print(attribution)
# Compare to heuristic modelsfrom pychattr.channel_attribution import HeuristicModel
heuristic = HeuristicModel(    paths=journeys,    conversions=conversions,    revenues=revenues)
# Get last-touch, first-touch, linearresults = heuristic.attribution(    heuristic_type=['last_touch', 'first_touch', 'linear'])

# Install# pip install CausalPy
import pandas as pdimport causalpy as cp
# Synthetic control examplesc = cp.pymc_models.SyntheticControl(    data=data,    treatment_time=70,    formula="actual ~ 0 + a + b + c + d + e + f + g",    model=cp.pymc_models.WeightedSumFitter(        sample_kwargs={"draws": 2000, "target_accept": 0.95}    ))
# Get resultsresult = sc.fit()
# Plotresult.plot()
# Get causal impactimpact = result.summary()print(f"Average Treatment Effect: {impact['causal_impact'].mean():.2f}")

# Install# pip install dowhy
import dowhyfrom dowhy import CausalModel
# Define causal modelmodel = CausalModel(    data=df,    treatment='ad_spend',    outcome='revenue',    common_causes=['seasonality', 'competitor_activity'],    instruments=['budget_shock'])
# Identify causal effectidentified_estimand = model.identify_effect(proceed_when_unidentifiable=True)
# Estimate causal effectestimate = model.estimate_effect(    identified_estimand,    method_name="backdoor.propensity_score_matching")
print(f"Causal Effect: {estimate.value}")
# Refute resultsrefute = model.refute_estimate(    identified_estimand,    estimate,    method_name="random_common_cause")

# Install# pip install econml
from econml.dml import LinearDMLfrom sklearn.ensemble import GradientBoostingRegressor
# Double ML for heterogeneous treatment effectsdml = LinearDML(    model_y=GradientBoostingRegressor(),    model_t=GradientBoostingRegressor())
# Fit modeldml.fit(    Y=outcomes,  # Revenue    T=treatment,  # Ad spend    X=features   # Customer features, context)
# Get treatment effect for new datatreatment_effect = dml.effect(X_test)
# Get confidence intervalslb, ub = dml.effect_interval(X_test, alpha=0.05)

import numpy as npfrom scipy import stats
def cuped_variance_reduction(    pre_experiment_metric,    experiment_metric_control,    experiment_metric_treatment):    """    Implement CUPED variance reduction using pre-experiment data    Based on: https://exp-platform.com/Documents/2013-02-CUPED-ImprovingSensitivityOfControlledExperiments.pdf    """        # Combine control and treatment    Y = np.concatenate([experiment_metric_control, experiment_metric_treatment])    X = np.concatenate([pre_experiment_metric[:len(experiment_metric_control)],                        pre_experiment_metric[-len(experiment_metric_treatment):]])        # Compute theta (optimal coefficient)    theta = np.cov(X, Y)[0, 1] / np.var(X)        # Adjust metrics    Y_adjusted = Y - theta * (X - np.mean(X))        # Split back    n_control = len(experiment_metric_control)    Y_adj_control = Y_adjusted[:n_control]    Y_adj_treatment = Y_adjusted[n_control:]        # Compute variance reduction    original_var = np.var(Y)    adjusted_var = np.var(Y_adjusted)    variance_reduction = 1 - (adjusted_var / original_var)        # Run t-test    t_stat, p_value = stats.ttest_ind(Y_adj_treatment, Y_adj_control)        return {        'treatment_effect': np.mean(Y_adj_treatment) - np.mean(Y_adj_control),        'p_value': p_value,        'variance_reduction': variance_reduction,        'theta': theta    }

"""Unified measurement approach combining MMM, attribution, and experiments"""import pandas as pdfrom pymc_marketing.mmm import MMMfrom lightweight_mmm import preprocessing
class UnifiedMeasurement:    def __init__(self):        self.mmm = None        self.attribution_model = None        self.geo_test_results = {}        def fit_mmm(self, aggregate_data):        """Fit media mix model on aggregate time series"""        self.mmm = MMM(            date_column="date",            channel_columns=["tv", "radio", "digital", "social"],            adstock_max_lag=8        )                self.mmm.fit(            data=aggregate_data,            target_column="revenue",            tune=1000,            draws=1000        )                return self.mmm        def calibrate_with_experiments(self, experiment_results):        """        Use geo experiment results to calibrate MMM        Implements approach from Meta's Robyn        """        for channel, result in experiment_results.items():            # Get experiment lift            exp_lift = result['incremental_revenue'] / result['spend']                        # Get MMM prediction for same period            mmm_lift = self.mmm.compute_channel_roi()[channel]                        # Calculate calibration factor            calibration = exp_lift / mmm_lift                        # Store for adjustment            self.geo_test_results[channel] = {                'experiment_roi': exp_lift,                'model_roi': mmm_lift,                'calibration_factor': calibration            }        def optimize_budget(self, total_budget, constraints):        """Optimize budget allocation using calibrated model"""                # Get base optimization        base_opt = self.mmm.optimize_budget(            total_budget=total_budget,            budget_bounds=constraints        )                # Adjust with experimental calibration        calibrated_allocation = {}        for channel, amount in base_opt.items():            if channel in self.geo_test_results:                factor = self.geo_test_results[channel]['calibration_factor']                calibrated_allocation[channel] = amount * factor            else:                calibrated_allocation[channel] = amount                # Normalize to budget        total = sum(calibrated_allocation.values())        calibrated_allocation = {            k: v * (total_budget / total)             for k, v in calibrated_allocation.items()        }                return calibrated_allocation

"""Complete geo experiment workflow"""
def design_geo_test(historical_data, treatment_geos, test_duration_weeks):    """Design matched-market test"""    from matched_markets.methodology import trimmed_match        # Find best control geos    design = trimmed_match.TrimmedMatch(        data=historical_data,        treatment_geos=treatment_geos,        pre_treatment_periods=52  # 1 year lookback    )        design.fit()        # Power analysis    power_results = design.power_analysis(        effect_sizes=np.arange(0.05, 0.30, 0.05),        test_duration=test_duration_weeks    )        return design, power_results
def analyze_geo_test(design, post_treatment_data):    """Analyze completed test"""        # Estimate treatment effect    result = design.estimate_treatment_effect(        post_treatment_response=post_treatment_data['revenue'],        post_treatment_spend=post_treatment_data['spend']    )        # Calculate iROAS    iroas = result.incremental_revenue / result.incremental_spend        # Check statistical significance    is_significant = result.p_value < 0.05        return {        'iroas': iroas,        'incremental_revenue': result.incremental_revenue,        'p_value': result.p_value,        'confidence_interval': (result.ci_lower, result.ci_upper),        'is_significant': is_significant    }

# Example: Grid search for adstock and saturation parametersfrom lightweight_mmm import optimize_media
# Define parameter rangesadstock_range = {    'tv': (0.3, 0.8),    'digital': (0.1, 0.5),    'radio': (0.2, 0.6)}
saturation_range = {    'tv': (0.5, 1.5),    'digital': (0.3, 1.0),    'radio': (0.4, 1.2)}
# Run optimizationbest_params = optimize_media.optimize_hyperparameters(    media_data=media_data,    target=target,    adstock_ranges=adstock_range,    saturation_ranges=saturation_range,    n_iter=100)

"""Prior elicitation using domain expertiseBased on: https://github.com/louismagowan/mmm-prior-elicitation"""
def elicit_channel_priors(historical_roas_estimates):    """    Convert historical ROAS ranges into prior distributions    """    priors = {}        for channel, roas_range in historical_roas_estimates.items():        # Convert ROAS range to lognormal parameters        lower, upper = roas_range        mean_log = (np.log(lower) + np.log(upper)) / 2        std_log = (np.log(upper) - np.log(lower)) / 4  # ~95% in range                priors[channel] = {            'distribution': 'lognormal',            'mu': mean_log,            'sigma': std_log        }        return priors
# Example usageexpert_estimates = {    'tv': (1.5, 3.0),      # TV ROAS between 1.5x and 3.0x    'digital': (2.0, 5.0),    'radio': (1.0, 2.5)}
priors = elicit_channel_priors(expert_estimates)

# Check for common issuesimport arviz as az
# 1. Check R-hat (should be < 1.01)def check_convergence(mmm_trace):    summary = az.summary(mmm_trace)    problematic = summary[summary['r_hat'] > 1.01]    if len(problematic) > 0:        print("Convergence issues for:", problematic.index.tolist())        print("Solution: Increase warmup/draws or reparameterize")    return problematic
# 2. Check for scaling issuesdef check_data_scaling(media_data, target):    for i, channel in enumerate(media_data.columns):        ratio = target.std() / media_data.iloc[:, i].std()        if ratio > 100 or ratio < 0.01:            print(f"{channel}: scaling issue (ratio: {ratio:.2f})")            print("Solution: Normalize or standardize data")
# 3. Increase samplingmmm.fit(    media=media_data,    target=target,    number_warmup=2000,  # Increased from 1000    number_samples=2000,    number_chains=4)

def diagnose_low_power(design, target_power=0.8):    """    Diagnose why geo test has low power    """    diagnostics = {}        # Check match quality    diagnostics['match_quality'] = design.match_quality_score    if diagnostics['match_quality'] < 0.7:        print("Poor match quality. Consider:")        print("- Increasing lookback window")        print("- Adding more control geos")        print("- Using time-based regression instead of matching")        # Check historical variance    diagnostics['outcome_cv'] = design.historical_cv    if diagnostics['outcome_cv'] > 0.3:        print("High outcome variance. Consider:")        print("- Longer test duration")        print("- Larger treatment effect")        print("- CUPED-style variance reduction")        # Check sample size    diagnostics['n_geos'] = len(design.treatment_geos) + len(design.control_geos)    if diagnostics['n_geos'] < 10:        print("Small sample size. Consider:")        print("- Aggregating smaller geos")        print("- Using synthetic control instead")        return diagnostics

def compare_attribution_models(journey_data):    """    Compare multiple attribution approaches    """    from pychattr.channel_attribution import MarkovModel, HeuristicModel        results = {}        # Heuristic models    heuristic = HeuristicModel(        paths=journey_data['paths'],        conversions=journey_data['conversions'],        revenues=journey_data['revenues']    )        results['last_touch'] = heuristic.attribution(heuristic_type='last_touch')    results['first_touch'] = heuristic.attribution(heuristic_type='first_touch')    results['linear'] = heuristic.attribution(heuristic_type='linear')        # Markov models    for order in [1, 2]:        markov = MarkovModel(            paths=journey_data['paths'],            conversions=journey_data['conversions'],            revenues=journey_data['revenues'],            order=order        )        markov.fit()        results[f'markov_order_{order}'] = markov.attribution()        # Compare    comparison = pd.DataFrame(results)    print("\nAttribution Comparison:")    print(comparison)        # Flag large discrepancies    for channel in comparison.index:        channel_results = comparison.loc[channel]        if channel_results.std() / channel_results.mean() > 0.5:            print(f"\nWarning: Large variation for {channel}")            print("Consider: combining with MMM or incrementality tests")        return comparison

from google.cloud import bigqueryimport pandas as pd
def load_marketing_data_for_mmm(project_id, lookback_days=730):    """    Load aggregated marketing data from BigQuery    """    client = bigquery.Client(project=project_id)        query = f"""    SELECT        DATE(timestamp) as date,        SUM(CASE WHEN channel = 'TV' THEN spend ELSE 0 END) as tv_spend,        SUM(CASE WHEN channel = 'Digital' THEN spend ELSE 0 END) as digital_spend,        SUM(CASE WHEN channel = 'Radio' THEN spend ELSE 0 END) as radio_spend,        SUM(revenue) as revenue,        SUM(conversions) as conversions    FROM `{project_id}.marketing.daily_summary`    WHERE DATE(timestamp) >= DATE_SUB(CURRENT_DATE(), INTERVAL {lookback_days} DAY)    GROUP BY date    ORDER BY date    """        df = client.query(query).to_dataframe()    return df