Alternative Response Curve Visualisation - Paired Plots : Request for Feedback

[louismagowan]

Hi there! Would love to hear some thoughts / feedback on this one.

Response curves in MMM are quite a tricky concept and it doesn't feel like the topic has been fully "solved" yet. Balancing actionability and ease of understanding for stakeholders with reliability and faithfulness of insights is difficult.

One solution that we do at my company (see here for video of it) is to create paired plots. for our response curves.

1. The top plot contains adstocked input variable against the channel contribution (i.e. the channel's "response")
2. The bottom plot contains the raw input variable, plotted over time
3. The colour coding of the plot points links adstocked input to the raw version that generated it

What do you guys think? Is there anything glaringly wrong or incorrect with this approach? Is it a useful and correct way of conveying insights from MMM?

Some rough/untidy code to reproduce those graphs is demoed below

import pandas as pd
from pymc_marketing.mmm import (
    DelayedAdstock,
    MichaelisMentenSaturation,
    MMM,
)
from pymc_marketing.mmm.transformers import delayed_adstock
from pymc_marketing.mmm.utils import michaelis_menten
import arviz as az
from scipy.optimize import curve_fit
import numpy as np

# Initial plotting configs
import matplotlib.pyplot as plt
az.style.use("arviz-darkgrid")
plt.rcParams["figure.dpi"] = 100



# Usual PyMC-Marketing model code
data_url = "https://raw.githubusercontent.com/pymc-labs/pymc-marketing/main/data/mmm_example.csv"
data = pd.read_csv(data_url, parse_dates=["date_week"])

X = data.drop("y",axis=1)
X[["x1", "x2"]] = X[["x1", "x2"]]  * 1000000 # Multiply spend by 1mil for pretend-reasons, just to demo CAC being in 100s
y = data["y"]
mmm = MMM(
    adstock=DelayedAdstock(l_max=8),
    saturation=MichaelisMentenSaturation(),
    date_column="date_week",
    channel_columns=["x1", "x2"],
    control_columns=[
        "event_1",
        "event_2",
        "t",
    ],
    yearly_seasonality=2,
)
mmm.fit(X, y)

# Generate posterior predictive samples
posterior_preds = mmm.sample_posterior_predictive(X)

# Get original scale posterior predictions for full dataset
posterior_preds = mmm.sample_posterior_predictive(
                        X,
                         var_names=["y",
                                    "channel_contributions", 
                                    "control_contributions", 
                                    "yearly_seasonality_contribution", 
                                    ],
                         progressbar=False,
                         extend_idata=False,
                         original_scale=True,
                         )

# Get the various contributions in the original scale
channel_contribs = pd.DataFrame(posterior_preds.channel_contributions.mean(dim="sample").values, columns=posterior_preds.channel_contributions.channel.values, index=posterior_preds.date.values)
control_contribs = pd.DataFrame(posterior_preds.control_contributions.mean(dim="sample").values, columns=posterior_preds.control_contributions.control.values, index=posterior_preds.date.values)
yearly_seasonality_contribs = pd.DataFrame(posterior_preds.yearly_seasonality_contribution.mean(dim="sample").values, columns=["yearly_seasonality"], index=posterior_preds.date.values)

# Merge into one df
all_contribs = pd.concat([channel_contribs, control_contribs, yearly_seasonality_contribs], axis=1)

# Add the intercept, prediction and true values
intercept = mmm.get_target_transformer().inverse_transform(mmm.idata.posterior.intercept.values.reshape(-1, 1)).mean()
all_contribs["intercept"] = intercept
all_contribs["prediction"] = all_contribs.sum(axis=1)
all_contribs["true"] = y.values

# Create dataframe to store adstocked spends inside
adstocked_spend_by_channel = pd.DataFrame(index=posterior_preds.coords["date"],
                               columns=posterior_preds.coords["channel"])

# Generate adstocked spend for each channel
for channel in mmm.channel_columns:

    # Get the mean adstock parameters for that channel
    channel_mean_adstock_alpha = mmm.idata.posterior.sel(channel=channel)['adstock_alpha'].mean().values
    channel_mean_adstock_theta = mmm.idata.posterior.sel(channel=channel)['adstock_theta'].mean().values

    # Apply the adstock transformation
    channel_adstocked_spend = delayed_adstock(
        X[channel],
        alpha=channel_mean_adstock_alpha,
        theta=channel_mean_adstock_theta,
        l_max=mmm.adstock.l_max,
        normalize=False # Don't normalize, we plot raw spend in pair plot later
        ).eval()

    # Store the adstocked spend in the dataframe
    adstocked_spend_by_channel[channel] = channel_adstocked_spend

# Create response curves for each channel
for channel in mmm.channel_columns:

    # Create plotting dataset
    plot_df = pd.DataFrame({
        "adstocked_input": adstocked_spend_by_channel[channel],
        "target": all_contribs[channel],
        "raw_input": X[channel].values
        }).reset_index().rename(columns={"index": "date_week",
                                        "target": "y"})

    # Add extra date columns for categorical plotting
    plot_df["half_year"] = np.where(plot_df["date_week"].dt.month <= 6, 'H1', 'H2')
    # date_marker col is what we use to colour our data points with
    plot_df["date_marker"] = plot_df["date_week"].dt.year.astype(str)


    # For recent years - give extra granularity in date_marker
    recent_years_mask = plot_df.date_marker.isin(["2020", "2021"])
    plot_df.loc[recent_years_mask, "date_marker"] = (
        plot_df.loc[recent_years_mask, "half_year"] + " " + 
        plot_df.loc[recent_years_mask, "date_marker"]
        )

    # Generate the Michaelis-Menten curve

    # Adjust these values based on your data
    alpha_initial_estimate = plot_df["y"].max()
    lam_initial_estimate = plot_df["adstocked_input"].median()  # or some other appropriate value
    initial_guess = [alpha_initial_estimate, lam_initial_estimate]
    # Add bounds to the curve_fit function
    lower_bounds = [0, 0]  # alpha and lambda must be positive
    upper_bounds = [np.inf, np.inf]  # or set some reasonable upper limits
    fitted_params, _ = curve_fit(michaelis_menten, plot_df["adstocked_input"], plot_df["y"], 
                                p0=initial_guess, bounds=(lower_bounds, upper_bounds), maxfev=5000)

    # Unpack optimal values for M-M params
    alpha_limit, lam_constant = fitted_params

    # Generate x_fit and y_fit - creating values for M-M curve
    plot_df["x_fit"] = np.linspace(0, plot_df["adstocked_input"].max(), plot_df.shape[0])
    plot_df["y_fit"] = michaelis_menten(plot_df.x_fit, alpha_limit, lam_constant)

    ### Plot Response curve
    # Setting up the figure and subplots with specified figsize and horizontal spacing
    fig, axs = plt.subplots(2, 1, figsize=(10, 10), gridspec_kw={'hspace': 0.3})

    # Extracting unique years for color mapping across the plots
    markers = plot_df['date_marker'].unique()
    colors = plt.cm.viridis(np.linspace(0, 1, len(markers)))

    ##### First Plot: Adstocked Spend vs. Sales
    ax = axs[0]
    # Looping through the unique years to create a scatter plot for each year
    for i, marker in enumerate(markers):
        # Filter data for the current year
        yearly_data = plot_df[plot_df['date_marker'] == marker]
        
        # Creating a scatter plot for each year's data
        ax.scatter(yearly_data['adstocked_input'], yearly_data["y"], alpha=0.9, color=colors[i], label=marker)
        
    # Plotting the fitted response curve
    ax.plot(plot_df.x_fit, plot_df.y_fit, color="red", alpha=0.6, label="\nFitted Response\nCurve")


    # Setting labels and title for the first plot
    ax.set_xlabel(f'Adstocked {channel}')
    ax.set_ylabel('Incremental Target')
    ax.set_title(f'Response Curve')

    ##### Second Plot: Raw Spend Over Time
    ax = axs[1]  # Reference to the second subplot
    for i, year in enumerate(markers):
        yearly_data = plot_df[plot_df['date_marker'] == year]
        ax.scatter(yearly_data["date_week"], yearly_data["raw_input"], label=year, color=colors[i])

    # Setting labels and title for the second plot
    ax.set_xlabel('Week')
    ax.set_ylabel(f'Raw {channel}')
    ax.set_title(f'Raw {channel} Over Time')

    # Setting the overall figure title
    # Calculate overall CAC for channel
    cac = plot_df["raw_input"].sum() / plot_df["y"].sum() # Multiply CAC for pretend-reasons, just to demo CAC being in 100s
    fig.suptitle(f"{channel} | CAC: {cac:.2f}", x=0.53, y=1.02, color="tab:purple", fontname="monospace", fontweight="bold")

    # Creating a common legend for the figure, placed to the right of the plots
    handles, labels = axs[0].get_legend_handles_labels()  # Gathering legend handles and labels from the first subplot
    legend = fig.legend(handles, labels, loc='upper right', bbox_to_anchor=(1.3, 0.6), title='Legend', frameon=True)

[louismagowan]

cc @iraur

[louismagowan]

@juanitorduz

[juanitorduz]

Dear @louismagowan, thank you for the input! At the moment we have two types of response curves:

• Direct ones where we simply plot spend vs mean posterior contribution. For example

• Total plots where we vary the spend on the whole timeline equally and plot the resulting contribution. For example

---

One thing we are missing from the direct plots is that we are not considering properly the adstock effect, which explains the contribution points when there is no spend.

It would be great to have these direct plots with the adsctock included as your first plot @louismagowan ! Moreover, it would be great to work with the whole posterior instead of taking the posterior mean. Does it make sense?

[louismagowan]

Thanks for the feedback @juanitorduz! 😁

And good point about the posterior dist rather than the mean. It sounds like there might be some material for a pull request in here then

1. Wrap adstocking inputs into a function
2. Add adstocked variation (bool toggle) into the response curve functions.

Also, what did you think about the colour-matching over time? What do you think are the advantages of a cost-share approach you have vs the implementation I added?

Thanks! ☺️

[juanitorduz]

What did you think about the colour-matching over time?

I think we can make it as optional as I know people that have a time reading these colors.

What do you think are the advantages of a cost-share approach you have vs the implementation I added?

I think these approaches solve different questions (global vs local) so it is nice to have them both (and more!). The more diagnostics, the better! Different audiences need different plots ;)

[louismagowan]

Ok clear - thank you! 😁

[iraur]

Hi @juanitorduz - I actually have a follow-up question on Louis' screenshot above which is: why would the plotted points not follow a perfect Michaelis-Menten curve in this case?

In the model we first adstock and then saturate our spend via Michaelis-Menten for each channel and so, if I understand well, the weekly contributions should sit nice and neatly on a Michaelis-Menten curve when plotted on adstocked spend. Whereas what we see in the plot above is a scattering around the curve as if some noise has been added. Any comments on that would be highly appreciated!

[juanitorduz]

Could it be an effect of taking the mean? We have uncertainty both in the x-axis (adstock alpha parameter) and the y-axis (saturation parameters)

(sorry about the late reply)

[PK1706]

I don't think that it is only due to taking the mean. I ran into the same problem when using the median as also discussed here.
Just as an example here is a plot of the Adstocked Spend against the Response (both calculated with the median):

A further explanation of this would be highly appreciated as I have tried to wrap my head around it for a long time.