Predicted Target ≠ Sum of Contributions OOS period (Case Study Notebook)

[RamiFisherTW]

(Reopen of issue 1826)

I'm experiencing a discrepancy between direct predictions and reconstructed predictions from model components when using sample_posterior_predictive.
Steps to Reproduce:

1. Generate test predictions with component decomposition:

y_pred_test = mmm.sample_posterior_predictive(
    X_pred=X_test[0:1],
    include_last_observations=True,
    original_scale=True,
    var_names=["y", "channel_contribution", "intercept", "fourier_contribution", 
               "control_contribution", "yearly_seasonality_contribution"],
    extend_idata=False,
    progressbar=True,
    random_seed=rng
)

2. Decompose components and reconstruct total prediction:

# Extract component contributions
control_contrib_pred_test = y_pred_test['control_contribution'].mean(dim=["sample"]).sum()
channel_contrib_pred_test = y_pred_test['channel_contribution'].mean(dim=["sample"]).sum()
intercept_pred_test = y_pred_test['intercept'].mean(dim=["sample"]) * mmm.get_target_transformer()["scaler"].scale_[0]
seasonality_pred_test = y_pred_test['yearly_seasonality_contribution'].mean(dim=["sample"]).sum()

# Reconstruct total prediction
total_prediction_test = control_contrib_pred_test + channel_contrib_pred_test + intercept_pred_test + seasonality_pred_test

Issue:

Reconstructed sum: 67,879,827.31
Direct prediction (y_pred_test['y'].mean(dim=["sample"]).sum()): 72,591,532.16
Difference: ~4,711,704.85

When using sample_response_distribution with the same allocation strategy:

response = mmm.sample_response_distribution(
    allocation_strategy=allocation_strategy,  # Using X_test[0] budget
    time_granularity=model_granularity,
    num_periods=num_periods,
    noise_level=0.0,
)

I get the same channel contribution value (27,098,727.90), which seems consistent.
Questions:

Why is there a discrepancy between direct prediction and the sum of components?
Are there any missing components in the reconstruction? - (It seems to be the intercept)
Why does budget_optimize only return channel contribution as the optimized value when I set 'total_contribution' as the response?

Any insights would be appreciated!

[juanitorduz]

Can you please try with

y_out_of_sample_with_adstock = mmm.sample_posterior_predictive(
X_out_of_sample, extend_idata=False, include_last_observations=True
)

As in the last part of the example notebook

https://www.pymc-marketing.io/en/latest/notebooks/mmm/mmm_example.html

Note the include last observations parameter

[RamiFisherTW]

It is the same code that I have wrote in the description

y_pred_test = mmm.sample_posterior_predictive(
X_pred=X_test[0:1],
include_last_observations=True,
original_scale=True,
var_names=["y", "channel_contribution", "intercept", "fourier_contribution",
"control_contribution", "yearly_seasonality_contribution"],
extend_idata=False,
progressbar=True,
random_seed=rng
)

Can you explain what should I change?

[RamiFisherTW]

Can you test the Case Study notebook on your side and check if you can get it right where y and reconstructed are equal?

[RamiFisherTW]

@juanitorduz Any update?

[juanitorduz]

@RamiFisherTW please understand open source developers have other responsibilities (we do not work on weekends) and try to look into all the issues as soon as we can.

This is in my to-do list. Thank you for understanding.

[juanitorduz]

Hi @RamiFisherTW ! I thing I found the issue: when we do posterior sampling, we are not scaling the intercept correctly. This should work (I will create a fix)

import xarray as xr

y_pred_test = mmm.sample_posterior_predictive(
    X_test,
    include_last_observations=True,
    original_scale=True,
    var_names=[
        "channel_contribution",
        "control_contribution",
        "fourier_contribution",
        "intercept",
        "mu",
        "y",
        "yearly_seasonality_contribution",
    ],
    extend_idata=False,
    progressbar=False,
    random_seed=rng,
)

# ! We compute this by hand.
intercept_original_scale = xr.apply_ufunc(
    mmm.target_transformer.inverse_transform,
    y_pred_test["intercept"].unstack(),
    input_core_dims=[["chain", "draw"]],
    output_core_dims=[["chain", "draw"]],
    vectorize=True,
) # <- This is not being handled correctly

test_decomposition_without_intercept = (
    y_pred_test["channel_contribution"].sum(dim="channel")
    + y_pred_test["yearly_seasonality_contribution"]
    + y_pred_test["control_contribution"].sum(dim="control")
)


test_error = (
    y_pred_test["y"] - test_decomposition_without_intercept
).unstack().transpose(..., "date") - intercept_original_scale


az.plot_hdi(x=X_test[date_column].to_numpy(), y=test_error, smooth=False, hdi_prob=0.94)

[RamiFisherTW]

Thanks!
According to error plot there is still very high absolute error , what have I got wrong?