Chapter 8 - Conditional Manatees

Chapter 8 - Conditional Manatees#

Key Takeaway#

Intro#

A quiet assumption runs through most statistical analysis is that our data is a faithful, unbiased sample of the thing we’re trying to understand. Unfortunately, the reality is never quite so simple. Data enters our sample through a process, be it a survey methodology, a survival filter, or a reporting system, and it can shape our conclusions if we’re not paying attention. But if we’re being intellectually honest and think critically about our data, what we’ll find is that our data is oftens conditional on how it ended up in our sample in the first place.

For example in Florida, manatees are often discovered to bear propeller scars on their backs as they’re creatures who often feed near the surface of the water as they need to come up for oxygen every 15 minutes or so. The intuitive response is to require propeller guards on boats. However, autopsies reveal that collisions with blunt parts of the hull called the keel actually cause far more manatee deaths than propellers do. The reason propeller scars are so common among surviving manatees is precisely because while a propeller strike is gruesome, cruel, and surely painful, it is also less fatal than a keel strike. Therefore the manatees killed by keels simply aren’t in the sample of living, but scarred, animals we observe in the wild.

With the manatee cause of death, the evidence was misleading because it is conditional on survival which exemlifies one of the key themes of this chapter: That every inference is conditional on something and that statistical modeling is largely the art of making that conditionality explicit and principled.

Conditioning is one of the most important principles of statistical inference and appears at multiple levels:

Data, like the manatee scars, enter our sample conditionally, such as whether it lived or died;
Posterior distributions are conditional on data;
And all model-based inference is conditional on the model.

For every linear model we’ve built so far, we assumed that every predictor had an independent association with the mean of the outcome. But what if we wanted to express how the association between predictor and outcome can be conditional? For example with our milk example on the relationship between milk energy and brain size in Chapter 5, what if we wanted our predictors to be conditional on the taxonomic group they belong to between apes, monkeys, and prosimians? This idea of the importance of a predictor being dependent on another predictor is what’s known as interaction or moderation. Interaction is a kind of conditioning which allow parameters (really their posterior distributions) to be conditional on further aspects of the data and is something we commonly see in multilevel models for example.

In this chapter, we will be exploring how interaction works in terms of specifying, interpreting, and plotting them for our understanding.

import warnings

import arviz as az
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns

from scipy import stats

pd.options.mode.chained_assignment = None
warnings.filterwarnings("ignore", category=FutureWarning)

%config Inline.figure_format = 'retina'
az.style.use("arviz-darkgrid")
az.rcParams["stats.hdi_prob"] = 0.89  # set credible interval for entire notebook
az.rcParams["stats.information_criterion"] = "waic"  # set information criterion to use in `compare`
az.rcParams["stats.ic_scale"] = "deviance"  # set information criterion scale
np.random.seed(0)

8.1 Building an Interaction#

Consider a regression of log GDP per capita against terrain ruggedness, run separately for African and non-African nations. As an aside, we’re using the “log” of GDP per capita rather than the raw values themselves because it tends to be an exponential property as wealth generates wealth. Outside Africa, the slope is negative as rugged terrain hinders transport, restricts market access, and ultimately drags down economic productivity which makes intuitive sense as to why ‘ruggedness’ brings down GDP. But within Africa, however, we see the opposite where the slope is positive. So we need to ask: Why would the rugged terrain be associated with higher GDP in Africa?

The leading hypothesis on this question is historical: Rugged regions were harder for slavers to raid and invade. The Atlantic and Indian Ocean slave trades devastated accessible coastal populations for generations, leaving persistent economic scars. Regions protected by terrain suffered less and that protection may still show up in economic data today. Whether this is truly causal is uncertain. Ruggedness correlates to other important factors such as coastlines, distance to sea, and other confounds but to see the reversal of slope is intellectually curious and demands an explanation.

Let’s generate a hypothesis of the potential causal DAG at hand: Terrain ruggedness ® and which Continent © a country belongs to both influence GDP (G), with unobserved confounders (U) lurking in the background. Note that the double circle is meant to communicate that the variable is an unobserved, or latent variable. The DAG doesn’t tell us how R and C combine to produce G. That’s what the statistical model has to figure out.

# Build the DAG
G = nx.DiGraph()
G.add_nodes_from(["R", "G", "C", "U"])
G.add_edges_from([
    ("R", "G"),  # R → G
    ("C", "G"),  # C → G
    ("U", "G"),  # U → G
    ("U", "R"),  # U → R
])

# Node positions — match the original diagram layout
pos = {
    "R": (0.0,  0.0),
    "G": (1.0,  0.0),
    "C": (2.0,  0.0),
    "U": (0.5, -1.0),
}

fig, ax = plt.subplots(figsize=(6, 4), facecolor="white")
ax.set_facecolor("white")

# Observed nodes (single circle)
nx.draw_networkx_nodes(G, pos, nodelist=["R", "G", "C"], ax=ax,
                       node_color="white", edgecolors="black",
                       node_size=900, linewidths=1.8)

# Unobserved node U (double circle — drawn twice at different sizes)
nx.draw_networkx_nodes(G, pos, nodelist=["U"], ax=ax,
                       node_color="white", edgecolors="black",
                       node_size=1050, linewidths=1.8)
nx.draw_networkx_nodes(G, pos, nodelist=["U"], ax=ax,
                       node_color="white", edgecolors="black",
                       node_size=600, linewidths=1.2)

nx.draw_networkx_labels(G, pos, ax=ax,
                        font_size=13, font_weight="bold", font_color="black")

nx.draw_networkx_edges(G, pos, ax=ax,
                       edge_color="black", arrows=True,
                       arrowstyle="-|>", arrowsize=18, width=1.5,
                       node_size=1050,
                       min_source_margin=22, min_target_margin=22)

ax.axis("off")
plt.tight_layout()

The obvious temptation might be to split the data and fit separate regressions, one for Africa and another for non-African countries. However, there’s four reasons we could argue as to why this isn’t the best idea:

First, parameters like σ (the error standard deviation) plausibly don’t differ between Africa and the rest of the world. Splitting the data forces you to estimate σ separately for each subset, sacrificing precision — you’ve made an accidental assumption that the variance differs when there’s no scientific reason to believe it does.

Second, you need the continent variable inside the model to make probability statements about its predictive value. If you split the data, you never quantify how much uncertainty there is about whether “African” vs. “not African” matters at all.

Third, model comparison via information criteria (WAIC, PSIS) requires that the models being compared are fit to the same data. Splitting makes that impossible.

Fourth, once we move into multilevel models in Chapter 13, we’ll eventually find that models can actually borrow information across groups. What we learn about ruggedness outside Africa should, in principle, inform our estimates inside Africa and vice versa. Splitting the data walls off that borrowing entirely.

The key strategy we need to exploit here is to make the model split the data, not us.

Code 8.1#

d = pd.read_csv("https://raw.githubusercontent.com/vanislekahuna/Statistical-Rethinking-PyMC/refs/heads/main/Data/rugged.csv", delimiter=";")

# make log version of the outcome
d["log_gdp"] = np.log(d["rgdppc_2000"])

# extract countries with GDP data
dd = d.dropna(subset=["log_gdp"])

# rescale variables
dd["log_gdp_std"] = dd["log_gdp"] / dd["log_gdp"].mean()
dd["rugged_std"] = dd["rugged"] / dd["rugged"].max()

dd.head()

Figure 8.2. Separate linear regressions inside and outside of Africa, for log-GDP against terrain ruggedness.#

fig, axs = plt.subplots(1, 2, figsize=(12, 6), sharey=True)

sns.regplot(
    x=dd.loc[dd["cont_africa"] == 1]["rugged_std"],
    y=dd.loc[dd["cont_africa"] == 1]["log_gdp_std"],
    scatter_kws={"color": "b"},
    line_kws={"color": "k"},
    ax=axs[0],
)

sns.regplot(
    x=dd.loc[dd["cont_africa"] == 0]["rugged_std"],
    y=dd.loc[dd["cont_africa"] == 0]["log_gdp_std"],
    scatter_kws={"edgecolor": "k", "facecolor": "w"},
    line_kws={"color": "k"},
    ax=axs[1],
)

axs[0].set_ylabel("log GDP (as proportion of mean)")
axs[1].set_ylabel("")
axs[0].set_title("African nations")
axs[1].set_title("Non-African nations")

# label countries
for _, africa in dd.loc[(dd["rugged_std"] > 0.7) & (dd["cont_africa"] == 1)].iterrows():
    axs[0].text(
        africa["rugged_std"],
        africa["log_gdp_std"] - 0.02,
        africa["country"],
        ha="center",
    )

for _, non_africa in dd.loc[(dd["rugged_std"] > 0.7) & (dd["cont_africa"] == 0)].iterrows():
    axs[1].text(
        non_africa["rugged_std"] + 0.03,
        non_africa["log_gdp_std"],
        non_africa["country"],
        va="center",
    )

for ax in axs:
    ax.set_xlim(-0.1, 1.1)
    ax.set_xlabel("ruggedness (standardised)")


#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=0.5,
    y=-0.06,
    t="Figure 8.2. Separate linear regressions inside and outside of Africa, for log-GDP against terrain ruggedness. The slope is positive inside Africa, \n \
    but negative outside. How can we recover this reversal of the slope, using the combined data?",
    ma="left"
  );

8.1.1 Building the Rugged Model#

Before introducing interaction, let’s first fit a simple model regressing log GDP on ruggedness using all nations pooled, essentially ignoring continent. We’ll rescale the variables thoughtfully so log GDP is expressed as a proportion of the international mean (i.e. 1.0 = average, 0.8 = 80% of average, 1.5 = 50% above average), and ruggedness is divided by its maximum observed value (with 0 being a country with perfectly flat to 1 representing Lesotho, the most rugged country in the sample). These choices make it easier to assign interpretable priors.

We’ll use the same geocentric skeleton as we’ve done with past Simple Linear Bayesian Regression models:

\[\log(y_i) \sim \text{Normal}(\mu_i, \sigma)\]

\[\mu_i = \alpha + \beta(r_i - \bar{r})\]

The y-intercept α is the expected log GDP when ruggedness is at its sample mean. Since the outcome is scaled to have a mean of 1, the prior α ~ Normal(1, 0.1) is sensible as it’s basically saying the intercept is almost certainly between 0.8 and 1.2.

For the slope β, the initial guess of Normal(0, 1) signals we don’t have any bias towards there being a positive or negative relationship between ruggedness and GDP. However as we’ll see in Figure 8.3, our standard deviation for this prior is far too wide. Prior predictive simulation will reveal that this prior allows slopes so extreme that they predict GDP values far outside anything ever observed. A slope of 0.3 seems much more plausible as you can see with the graph on the right so let’s tighten our slope prior there to β ~ Normal(0, 0.3) instead.

The key lesson we need to learn here is that prior predictive simulation is a discipline, not an optional extra. Even when data are plentiful enough to wash out bad priors, the exercise of thinking through priors forces you to understand what your model is actually saying.

The fitted model (m_8_1t) shows essentially no overall association between ruggedness and log GDP which makes sense, because the positive Africa effect and the negative everywhere-else effect cancel out when pooled.

Code 8.2#

with pm.Model() as m_8_1:
    a = pm.Normal("a", 1, 1)
    b = pm.Normal("b", 0, 1)

    mu = a + b * (dd["rugged_std"].values - dd["rugged_std"].mean())
    sigma = pm.Exponential("sigma", 1)

    log_gdp_std = pm.Normal("log_gdp_std", mu, sigma, shape=dd.shape[0])

Code 8.3#

with m_8_1:
    m_8_1_prior = pm.sample_prior_predictive()

# Figure 8.3 is below

Code 8.4#

beta_prior = az.extract_dataset(m_8_1_prior.prior)["b"].to_numpy()
np.sum(np.abs(beta_prior > 0.6)) / len(beta_prior)

Code 8.5#

with pm.Model() as m_8_1t:
    a = pm.Normal("a", 1, 0.1)
    b = pm.Normal("b", 0, 0.3)

    mu = a + b * (dd["rugged_std"].values - dd["rugged_std"].mean())
    sigma = pm.Exponential("sigma", 1)

    log_gdp_std = pm.Normal("log_gdp_std", mu, sigma, observed=dd["log_gdp_std"].values)

    # m_8_1t_trace = pm.sample()
    m_8_1t_trace = pm.sample(idata_kwargs={"log_likelihood": True})

    m_8_1t_prior = pm.sample_prior_predictive()

Figure 8.3. Simulating in search of reasonable priors for terrain ruggedness.#

# Figure 8.3

_, (ax1, ax2) = plt.subplots(1, 2, figsize=[7, 4], constrained_layout=True)

n = 100
rugged_plot = np.linspace(-0.1, 1.1, n)

# Prior 1
prior = m_8_1_prior.prior.sel(draw=slice(None, None, int(len(m_8_1_prior.prior.draw) / n)))
reglines = prior["a"].T.to_numpy() + rugged_plot * prior["b"].T.to_numpy()
for regline in reglines:
    ax1.plot(
        rugged_plot,
        regline,
        color="k",
        lw=1,
        alpha=0.3,
    )
ax1.set_title("a ~ Normal(1, 1)\nb ~ Normal(0, 1)")

# Prior 2
prior_t = m_8_1t_prior.prior.sel(draw=slice(None, None, int(len(m_8_1t_prior.prior.draw) / n)))
reglines_t = prior_t["a"].T.to_numpy() + rugged_plot * prior_t["b"].T.to_numpy()

for regline in reglines_t:
    ax2.plot(
        rugged_plot,
        regline,
        color="k",
        lw=1,
        alpha=0.3,
    )
ax2.set_title("a ~ Normal(1, 0.1)\nb ~ Normal(0, 0.3)")

for ax in (ax1, ax2):
    ax.set_xlabel("ruggedness")
    ax.set_ylabel("log GDP (prop of mean)")
    ax.axhline(0.7, ls="dashed", color="k", lw=1)
    ax.axhline(1.3, ls="dashed", color="k", lw=1)
    ax.set_ylim(0.5, 1.5)


#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=0.5,
    y=-0.06,
    t="Figure 8.3. Simulating in search of reasonable priors for the terrain ruggedness example. \n \
    The dashed horizontal lines indicate the minimum and maximum observed GDP values. \n \
    Left: The first guess with very vague priors. \n \
    Right: The improved model with much more plausible priors.",
    ma="left"
  );

Code 8.6#

az.summary(m_8_1t_trace, kind="stats", round_to=2)

8.1.2 Adding an Indicator Variable Isn’t Enough#

The next step is to let African and non-African nations have different intercepts. Usually the conventional approach would be to add another predictor via a 0/1 indicator variable (which forces asymmetric uncertainty on the priors) and evolving it from a ‘Simple’ to a ‘Multiple Linear Regression’ model.

\[\mu_i = \alpha + \beta(r_i - \bar{r}) + \gamma A_i\]

Instead, let’s try something different by utilizing the index variable approach. Similar to what we did in Chapter 5, nations in Africa will get an index value of 1 while all others get a value of 2. The new model with an index variable becomes:

\[\mu_i = \alpha_{\text{cid}[i]} + \beta(r_i - \bar{r})\]

Now we have two intercepts, α₁ and α₂, with the same prior on each. Again, this is structurally identical to the indicator-variable approach in Chapter 5 but makes it straightforward to assign symmetric, sensible priors.

Code 8.7#

cid = pd.Categorical(dd["cont_africa"])

cid

Code 8.8#

with pm.Model() as m_8_2:
    a = pm.Normal("a", 1, 0.1, shape=cid.categories.size)
    b = pm.Normal("b", 0, 0.3)

    mu = a[np.array(cid)] + b * (dd["rugged_std"].values - dd["rugged_std"].mean())
    sigma = pm.Exponential("sigma", 1)

    log_gdp_std = pm.Normal("log_gdp_std", mu, sigma, observed=dd["log_gdp_std"].values)

    m_8_2_trace = pm.sample(idata_kwargs={"log_likelihood": True})

Reviewing WAIC#

Widely Applicable Information Criterion (WAIC) is a measure of a model’s expected out-of-sample predictive accuracy, computed entirely from the training data. The way it works is that for each observation, WAIC asks: How surprised would this model be, on average across its posterior, if it saw this data point again in a fresh sample?

Specifically, it takes the log of the average likelihood across posterior samples for each observation and takes the average of these to get to the log pointwise predictive density (lppd). WAIC then applies a penalty for model complexity, which is the effective number of parameters (pWAIC), estimated as the variance of the individual log-likelihoods across posterior samples.

More flexible models overfit in ways that show up as higher variance in those per-point likelihoods. Subtracting this penalty from the lppd (and multiplying by −2 to get the deviance scale) gives WAIC. Note: A lower WAIC = better expected predictive performance.

# We'll flatten chains and draws into a single "samples" axis for clarity. Final shape we want: (n_samples, n_observations)
# Shape is (n_chains, n_draws, n_obs) right now
log_posterior_likelihood_m1 = m_8_1t_trace.log_likelihood["log_gdp_std"].values

# We need to stack chains and draws together
print(f"The log posterior likelihood's shape (n_chains, n_draws, n_obs): {log_posterior_likelihood_m1.shape}")
n_chains, n_draws, n_obs = log_posterior_likelihood_m1.shape

# Shape is now (n_samples, n_obs)
log_likelihood_m1 = log_posterior_likelihood_m1.reshape(n_chains * n_draws, n_obs)

# Each row = one posterior draw, each column = one observation
print(f"n_samples : {n_chains * n_draws}")
print(f"n_obs     : {n_obs}")
print(f"log_likelihood_m1 new shape: {log_likelihood_m1.shape}")


# ─────────────────────────────────────────────────────────────────────────────
# STEP 1 — lppd  (log pointwise predictive density)
# For each observation i, average the likelihood across all posterior samples, then take the log of that average.
def log_mean_exp(x, axis=0):
    """
    For each observation i, we're subtracting the maximum value of the array from every element,
    computing  the exponential of these new, smaller numbers,
    calculating the average of these exponentiated values,
    taking the natural logarithm of that average,
    then adding the maximum value back to the final result to restore the original scale of your data.
    `squeeze()` simply removes any flat, empty dimensions (like turning a matrix of [[5]] into just the number 5)
    """
    x_max = x.max(axis=axis, keepdims=True)

    return (np.log(np.exp(x - x_max).mean(axis=axis)) + x_max.squeeze())

# lppd_per_obs: log of the average likelihood that the model assigns to each data point, averaging over all posterior samples
lppd_per_obs = log_mean_exp(log_likelihood_m1, axis=0)  # shape: (n_obs,)

# lppd = sum across all observations
lppd = lppd_per_obs.sum()

print(f"\n── lppd ──────────────────────────────────")
print(f"lppd per obs (first 5): {np.round(lppd_per_obs[:5], 4)}")
print(f"Total lppd             : {lppd:.4f}")


# ─────────────────────────────────────────────────────────────────────────────
# STEP 2 — pWAIC  (effective number of parameters / complexity penalty)
# For each observation i, compute the VARIANCE of the log-likelihood across all posterior samples.
# A model that is overfit will show high variance in how confidently it predicts each point,
# the posterior disagreement is the signal.
pWAIC_per_obs = log_likelihood_m1.var(axis=0, ddof=1)  # shape: (n_obs,)

pWAIC = pWAIC_per_obs.sum()
print(f"\n── pWAIC ─────────────────────────────────")
print(f"pWAIC per obs (first 5): {np.round(pWAIC_per_obs[:5], 4)}")
print(f"Total pWAIC            : {pWAIC:.4f}")


# ─────────────────────────────────────────────────────────────────────────────
# STEP 3 — WAIC
#
# WAIC = -2 * (lppd - pWAIC)
#
# The -2 scaling puts it on the deviance scale (lower = better),
# consistent with AIC and DIC. The penalty pWAIC is subtracted from
# lppd because a more complex model has higher variance in its
# log-likelihoods and therefore a larger penalty — shrinking the
# effective predictive accuracy down from the raw lppd.

WAIC = -2 * (lppd - pWAIC)

print(f"\n── WAIC ──────────────────────────────────")
print(f"lppd  : {lppd:.4f}")
print(f"pWAIC : {pWAIC:.4f}")
print(f"WAIC  : {WAIC:.4f}  (lower = better out-of-sample prediction)")


# ─────────────────────────────────────────────────────────────────────────────
# STEP 4 — Standard error of WAIC
#
# Because WAIC is a sum of n_obs independent pointwise terms, its
# standard error can be estimated as the standard deviation of those
# per-observation contributions, scaled by sqrt(n_obs).
# This gives you a sense of how uncertain the WAIC estimate itself is.
# ─────────────────────────────────────────────────────────────────────────────

waic_per_obs = -2 * (lppd_per_obs - pWAIC_per_obs)   # per-obs WAIC contribution
waic_se      = np.sqrt(n_obs * waic_per_obs.var(ddof=1))

print(f"WAIC SE: {waic_se:.4f}")


# ─────────────────────────────────────────────────────────────────────────────
# SANITY CHECK — compare against ArviZ's built-in WAIC
# Should match closely (minor differences due to ArviZ's internal conventions)
# ─────────────────────────────────────────────────────────────────────────────

az_waic = az.waic(m_8_1t_trace, var_name="log_gdp_std")
print(f"\n── ArviZ WAIC (sanity check) ─────────────")
print(az_waic)

Code 8.9#

Now let’s read the comparison for models m_8_1 and m_8_2 (An example of Code 8.9):

         WAIC  pWAIC  dWAIC  weight    SE   dSE
m8.2   -252.3    4.3    0.0     1.0  15.3    NA
m8.1   -188.7    2.7   63.5     0.0  13.3  15.15

Reading this row by row:

WAIC: m8.2 scores −252.3, m8.1 scores −188.7. On the deviance scale, more negative is better (lower deviance = better fit). m8.2 is clearly preferred.
pWAIC (effective parameters): m8.2 uses ~4.3 effective parameters (it has two intercepts, one slope, and σ), while m8.1 uses ~2.7. The extra complexity in m8.2 is justified by its predictive gain.
dWAIC (difference from best model): m8.1 is 63.5 WAIC units worse than m8.2. This is a large gap.
weight: Akaike-style model weights, which can be interpreted roughly as the probability that a given model will make the best predictions on new data drawn from the same process. m8.2 gets essentially all the weight (1.0); m8.1 gets none (0.0).
SE: The standard error of each model’s WAIC estimate, based on the variability of per-observation contributions. m8.2’s SE is 15.3.
dSE: The standard error of the difference in WAIC between the two models, which accounts for correlation in per-observation scores. Here dSE = 15.15, but the difference itself is 63.5 — more than four standard errors. That’s very strong evidence that the continent variable is doing real predictive work.

az.compare({"m_8_1t": m_8_1t_trace, "m_8_2": m_8_2_trace})

Code 8.10#

What does m_8_2_trace actually buy us? Looking at az.summary() output, α₁ (Africa) ≈ 0.88 and α₂ (non-Africa) ≈ 1.05. African nations have reliably lower average log GDP. The 89% posterior interval for the contrast α₁ − α₂ runs from about −0.20 to −0.14: solidly negative. So m_8_2_trace correctly captures that African nations are, on average, poorer.

But here is the crucial limitation: m_8_2_trace has only adjusted the intercept. It fits parallel regression lines — the slope with respect to ruggedness is the same for both groups. That means it still cannot reproduce the reversal of slope we saw in Figure 8.2. It’s a better model than m_8_1t_trace, but it’s still the wrong model for the question we care about.

az.summary(m_8_2_trace, kind="stats", round_to=2)

Code 8.11#

m_8_2_posterior = az.extract_dataset(m_8_2_trace.posterior)
diff_a0_a1 = m_8_2_posterior["a"][1, :] - m_8_2_posterior["a"][0, :]
az.hdi(diff_a0_a1.to_numpy())

Figure 8.4. Including an indicator for African nations has not effect on the slope.#

Code 8.12#

fig, ax = plt.subplots()

# extract posterior samples of parameters
a0 = m_8_2_posterior["a"][0, :].to_numpy()
a1 = m_8_2_posterior["a"][1, :].to_numpy()
b = m_8_2_posterior["b"].to_numpy()

rugged_plot = np.linspace(-0.1, 1.1)

ax.scatter(
    dd.loc[cid == 0, "rugged_std"],
    dd.loc[cid == 0, "log_gdp_std"],
    label="Not Africa",
    facecolor="w",
    lw=1,
    edgecolor="k",
)
pred0 = a0 + rugged_plot.reshape(-1, 1) * b
ax.plot(rugged_plot, pred0.mean(1), color="grey")
az.plot_hdi(rugged_plot, pred0.T, color="grey", hdi_prob=0.97, ax=ax)

ax.scatter(
    dd.loc[cid == 1, "rugged_std"],
    dd.loc[cid == 1, "log_gdp_std"],
    label="Africa",
    color="b",
)
pred1 = a1 + rugged_plot.reshape(-1, 1) * b
ax.plot(rugged_plot, pred1.mean(1), color="b")
az.plot_hdi(rugged_plot, pred1.T, color="b", hdi_prob=0.97, ax=ax, fill_kwargs={"alpha": 0.2})

ax.legend(frameon=True)

ax.set_xlim(-0.1, 1.1)
ax.set_ylim(0.7, 1.3)
ax.set_xlabel("ruggedness (standardised)")
ax.set_ylabel("log GDP (as proportion of mean)")

#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=1.32,
    y=.75,
    t="Figure 8.4. Including an indicator for \n \
    African nations has no effect on the slope. \n \
    African nations are shown in blue. \n \
    Non-African nations are shown in black. \n \
    Regression means for each subset of nations are \n \
    shown in corresponding colors, along with \n \
    97% intervals shown by shading.",
    ma="left"
  )

8.1.3 Adding an Interaction Does Work#

To recover the slope reversal, we need to make the slope itself conditional on continent since this is the interaction. In index-variable form:

\[\mu_i = \alpha_{\text{cid}[i]} + \beta_{\text{cid}[i]}(r_i - \bar{r})\]

Now both the intercept and the slope differ by continent. There is a conventional alternative using an indicator variable and an explicit interaction term (β + γAᵢ), but that formulation makes it harder to assign sensible priors since the prior on γ effectively says the slope inside Africa is more uncertain than the slope outside, which is an arbitrary asymmetry. The index approach avoids this by placing the same Normal(0, 0.3) prior on both slopes.

Code 8.13#

with pm.Model() as m_8_3:
    a = pm.Normal("a", 1, 0.1, shape=cid.categories.size)
    b = pm.Normal("b", 0, 0.3, shape=cid.categories.size)

    mu = a[np.array(cid)] + b[np.array(cid)] * (dd["rugged_std"].values - dd["rugged_std"].mean())
    sigma = pm.Exponential("sigma", 1)

    log_gdp_std = pm.Normal("log_gdp_std", mu, sigma, observed=dd["log_gdp_std"].values)

    m_8_3_trace = pm.sample(idata_kwargs={"log_likelihood": True})

Code 8.14#

After fitting model m8.3, the key output is:

b[0] -0.14   (non-Africa)
b[1]  0.13   (Africa)

The slope is reversed: positive inside Africa (+0.13), negative outside (−0.14). This is exactly the pattern we set out to recover.

az.summary(m_8_3_trace, kind="stats", round_to=2)

Code 8.15#

The PSIS comparison of all three models makes the case decisively: m_8_3 captures more than 95% of the model weight. The interaction is strongly favored for prediction. However, we should be careful to note that predictive superiority is not the same as causal identification. Model comparison tells you what features matter for prediction in the sample, not what the true causal structure is.

az.compare({"m_8_1t": m_8_1t_trace, "m_8_2": m_8_2_trace, "m_8_3": m_8_3_trace}, ic="loo")

Code 8.16#

m_8_3_loo = az.loo(m_8_3_trace, pointwise=True)

plt.plot(m_8_3_loo.loo_i)

8.1.4 Plotting the Interaction#

Once an interaction is in play, plotting requires care. The approach here is straightforward: Make two separate bivariate plots, one for African nations and one for non-African nations. In each plot, the regression line and its uncertainty interval are computed from the posterior, fixing the continent index to the appropriate value and varying ruggedness over its observed range.

Figure 8.5 shows the result: The slope clearly reverses direction between the two panels. By achieving this within a single unified model, we can also quantify the statistical evidence for the reversal — something impossible if we had simply split the data.

Code 8.17#

m_8_3_posterior = az.extract_dataset(m_8_3_trace.posterior)

Figure 8.5. Posterior predictions for the terrain ruggedness model#

_, axs = plt.subplots(
    1,
    2,
    figsize=[8, 4],
    sharey=True,
    constrained_layout=True,
)

ax1, ax0 = axs

# extract posterior samples of parameters
a0 = m_8_3_posterior["a"][0, :].to_numpy()
a1 = m_8_3_posterior["a"][1, :].to_numpy()
b0 = m_8_3_posterior["b"][0, :].to_numpy()
b1 = m_8_3_posterior["b"][1, :].to_numpy()

rugged_plot = np.linspace(-0.1, 1.1)

ax0.scatter(
    dd.loc[cid == 0, "rugged_std"],
    dd.loc[cid == 0, "log_gdp_std"],
    label="Not Africa",
    facecolor="w",
    lw=1,
    edgecolor="k",
)
# calculating predicted manually because this is a pain with categorical variabiles in PyMC
pred0 = a0 + rugged_plot.reshape(-1, 1) * b0
ax0.plot(rugged_plot, pred0.mean(1), color="grey")
az.plot_hdi(rugged_plot, pred0.T, color="grey", hdi_prob=0.97, ax=ax0)
ax0.set_title("Non-African Nations")

ax1.scatter(
    dd.loc[cid == 1, "rugged_std"],
    dd.loc[cid == 1, "log_gdp_std"],
    label="Africa",
    color="b",
)
# calculating predicted manually because this is a pain with categorical variabiles in PyMC
pred1 = a1 + rugged_plot.reshape(-1, 1) * b1
ax1.plot(rugged_plot, pred1.mean(1), color="k")
az.plot_hdi(
    rugged_plot,
    pred1.T,
    color="blue",
    hdi_prob=0.97,
    ax=ax1,
    fill_kwargs={"alpha": 0.2},
)
ax1.set_title("African Nations")


ax.set_xlim(-0.1, 1.1)
ax0.set_xlabel("ruggedness (standardised)")
ax1.set_xlabel("ruggedness (standardised)")
ax0.set_ylabel("")
ax1.set_ylabel("log GDP (as proportion of mean)")

#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=0.4,
    y=-0.06,
    t="Figure 8.5. Posterior predictions for the terrain ruggedness model, \n \
    including the interaction between Africa and ruggedness. \n \
    Shaded regions are 97% posterior intervals of the mean.",
    ma="left"
  );

8.2. Symmetry of interactions.#

Within the GDP-ruggedness interaction question, there’s two equally valid questions we can ask about the interaction between the two variables:

How much does the association between ruggedness and log GDP depend on whether the nation is in Africa?
How much does the association between being in Africa and GDP depend on ruggedness?

To most humans, these sound like different questions as we tend to think of one variable as “the predictor” and the other as “the moderator.” However, a statistical model treats them identically.

Here the interaction term is symmetric as you cannot specify a model in which ruggedness moderates the effect of continent without simultaneously specifying a model in which continent moderates the effect of ruggedness. That may have been a mouthful so let’s consider the following mathematical expression:

\[ \mu_i = \underbrace{(2 - \text{CID}_i)(\alpha_1 + \beta_1(r_i - \bar{r}))}_{\text{CID}[i]=1} + \underbrace{(\text{CID}_i - 1)(\alpha_2 + \beta_2(r_i - \bar{r}))}_{\text{CID}[i]=2} \]

The equation may appear complex at first glance, but it is simply a conditional model that activates one of two parameter sets based on the value of $CID_i$. When $CID_i = 1$, only the Africa parameters (α₁, β₁) contribute while the second term drops to zero. When $CID_i = 2$, the reverse occurs where only the second term active. This means that reassigning a nation to Africa changes the prediction in a way that depends on that nation’s ruggedness… Unless its ruggedness happens to equal the mean, r̄, wherein which case the effect is uniform regardless.

Let’s plot the reverse interpretation of this equation where the association of being in Africa with log GDP dpeends upon terrain ruggedness. We can demonstrate this by plotting the difference in expected log GDP between a hypothetical nation inside Africa and the same nation outside Africa, holding its ruggedness constant (Figure 8.6).

Figure 8.6. The other side of the interaction between ruggedness and continent.#

Code 8.18#

fig, ax = plt.subplots(figsize=(6, 5))

rugged_plot = np.linspace(-0.1, 1.1)

delta = pred1 - pred0  # using 'pred' from above

ax.plot(rugged_plot, delta.mean(1), c="k")
az.plot_hdi(rugged_plot, delta.T, ax=ax, color="grey")

ax.axhline(0, ls="dashed", zorder=1, color=(0, 0, 0, 0.5))
ax.text(0.01, 0.01, "Africa higher GDP")
ax.text(0.01, -0.03, "Africa lower GDP")

ax.set_xlabel("ruggedness")
ax.set_ylabel("expected difference log GDP")
ax.set_xlim(0, 1)

#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=1.34,
    y=.75,
    t="Figure 8.6. The other side of the interaction \n \
    between ruggedness and continent. The vertical \n \
    axis is the difference in expected proportional \n \
    log GDP for a nation in Africa and one outside \n \
    Africa. At low ruggedness, we expect “moving” \n \
    a nation to Africa to hurt its economy. \n \
    But at high ruggedness, the opposite is true. \n \
    The association between continent and economy \n \
    depends upon ruggedness, just as much as the \n \
    association between ruggedness and economy \n \
    depend upon continent.",
    ma="left"
  )

In this counter-factual plot, there are no actual data points here, only model predictions which reveals something not obvious from Figure 8.5 which is that African nations are almost always economically worse off than comparable non-African ones, except at very high ruggedness values, where the difference narrows to near zero.

The practical implication is that our causal interpretation of an interaction should be driven by the scientific question and what’s actually manipulable, not by the mathematical framing. We can more easily imagine changing a country’s terrain ruggedness (flattening hills, blasting tunnels) than changing its continent. The continent variable is, in any case, a proxy for historical and political factors which are the actual causal levers are those factors, not the geography per se.

8.3 Continuous Interactions#

In the previous section we allowed the slope between terrain ruggedness and GDP to differ depending on whether a nation was in Africa or not. We can describe this as an interaction where the effect of one variable is conditional on the value of another. However because the continent variable is a categorical distinction in our Africa example, the model had a natural and limited set of answers. These were the perfect conditions to illustrate categorical interactions where you could estimate one slope for African nations and another for everywhere else then compare them. Demanding as it was to set up, the interpretation stayed manageable because the conditioning variable carved the data into a small number of discrete groups.

Continuous interactions push this logic further where instead of a slope that takes one of two or three fixed values depending on a category, the slope now varies fluidly and continuously as a second predictor changes. The mathematics are essentially the same where you’re still making one parameter conditional on another, however there’s no discrete grouping to anchor the interpretation. You can’t just report two slopes and call it done. The effect of one predictor on the outcome is now a continuously moving target, shifting across the entire range of the other predictor simultaneously.

This is why a table of posterior means and standard deviations becomes nearly useless for understanding continuous interactions. The solution needs to be visual. Throughout this section we work through a concrete two-predictor example and introduce the triptych plot which is a panel of three figures, each showing the relationship between one predictor and the outcome while holding the other predictor fixed at a different representative value. Comparing the panels side by side reveals what no single figure could which is how the slope changes as the second variable moves across its range, precisely what a continuous interaction is in practice.

8.3.1 A Winter Flower#

The dataset we’ll load in the code cell below contains bloom sizes from tulips grown in a greenhouse experiment, crossing three levels of soil moisture (water: 1 = low, 2 = medium, 3 = high) with three levels of light exposure (shade: 1 = high light, 2 = medium, 3 = low light). There is also a bed variable indicating which section of the greenhouse a plant came from, but it is not used in the main models of this section. The bloom column will be our target variable.

Our hypothetical causal DAG structure can be simple: Water ($W$) and Shade ($S$) both influence Bloom ($B$):

$W → B ← S$

The scientific motivation for our interaction hypothesis is concrete: Photosynthesis requires both water and light. In the absence of light, additional water can’t do much and in the absence of water, additional light can’t do much either. The two inputs are complementary, not independent. Therefore, one way to model this interdependency between light and water is through interaction effects which is a good start. However, the issue with this approach is that we don’t know how each predictor influences the bloom of a flower because in theory.

Code 8.19#

tulips = "https://raw.githubusercontent.com/vanislekahuna/Statistical-Rethinking-PyMC/refs/heads/main/Data/tulips.csv"

d = pd.read_csv(tulips, delimiter=";")
d.head()

8.3.2. The models.#

Here, we’ll be comparing the following two models, m_8_4 and m_8_5.

With the model m_8_4, the aim is to illustrate the independent effect of water and shade but no interaction effects (i.e. main effects only):

$B_i \sim Normal(\mu_i,\sigma)$

$\mu_i=\alpha+\beta_W(W_i−\overline{W})+\beta_S(S_i−\overline{S})$

Code 8.20#

The data cleaning/transformation steps we’ll take before modelling will be to rescale our variables before fitting them. Our target, blooms, will be divided by their maximum observed value so that the outcome scales from 0 to 1, thereby preserving zero as a meaningful boundary. The predictors water and shade are then each transformed to be mean-centered, running from −1 to 1).

d["blooms_std"] = d["blooms"] / d["blooms"].max()
d["water_cent"] = d["water"] - d["water"].mean()
d["shade_cent"] = d["shade"] - d["shade"].mean()

Code 8.21#

a = stats.norm.rvs(0.5, 1, 10000)
sum((a < 0) | (a > 1)) / len(a)

Code 8.22#

Now with m_8_4, we need to be careful about our prior specification to complete the rest of the model. The $\alpha$ prior for the y-intercept reflects the model’s intuition that for when water and shade are both at their means, blooms should be around halfway to the maximum. And a standard deviation of 0.25 should put only 5% of the mass outside the valid distribution:

$α \sim Normal(0.5, 0.25)$

a = stats.norm.rvs(0.5, 0.25, 10000)
sum((a < 0) | (a > 1)) / len(a)

Code 8.23#

Slopes for water and shade should both get a prior of $\text{Normal}(0, 0.25)$ on on the reasoning that a slope of $0.5$, which would carry blooms from zero to maximum across the full range of either predictor, is already at the extreme of plausibility and should therefore sit about two standard deviations out. Remember we’re trying to fit weakly informative priors here so that we discourage overfitting.

The entire model for m_8_4 looks like the following statistical notation:

$B_i \sim \text{Normal}(\mu_i,\sigma)$

$\mu_i=\alpha+\beta_W(W_i−\overline{W})+\beta_S(S_i−\overline{S})$

$α \sim \text{Normal}(0.5, 0.25)$

$\beta_w \sim \text{Normal}(0.5, 0.25)$

$\beta_s \sim \text{Normal}(0.5, 0.25)$

$\sigma \sim \text{Exponential}(1)$

with pm.Model() as m_8_4:
    a = pm.Normal("a", 0.5, 0.25)
    bw = pm.Normal("bw", 0, 0.25)
    bs = pm.Normal("bs", 0, 0.25)

    mu = a + bw * d["water_cent"].values + bs * d["shade_cent"].values
    sigma = pm.Exponential("sigma", 1)

    blooms_std = pm.Normal("blooms_std", mu, sigma, observed=d["blooms_std"].values)

    m_8_4_trace = pm.sample()
    m_8_4_post = az.extract_dataset(m_8_4_trace.posterior)

\[B_i \sim Normal(\mu_i,\sigma)\]

\[\mu_i=\alpha+\beta_W(W_i−\overline{W})+\beta_S(S_i−\overline{S}) +\beta_{WS}W_iS_i \]

Code 8.24#

Moving on to our second model, m_8_5 which we’ll use illustrate the interaction between water and shade, remember that the goal is to let the effect of water on bloom depend on how much shade the plant is receiving and vice versa. The biological reality is that water and light are complementary inputs to photosynthesis. At low light, additional water can’t do much because the plant doesn’t have enough light to use it. At high light, water becomes the binding constraint. Neither variable acts independently of the other. The question here is how to encode this mathematically.

To recap, the non-interaction model for m_8_4 defined the mean ($\mu_i$) as:

$\mu_i = \alpha + \beta_W W_i + \beta_S S_i$

where the slope $\beta_W$ is a fixed constant representing the rate of change when shade is at its mean value, meaning that water always has the same effect regardless of the shade level.

To make that slope conditional on shade, we need to give it its own linear model by defining $\gamma_{W,i}$ as the effective slope of water for observation $i$:

$\gamma_{W,i} = \beta_W + \beta_{WS} S_i$

Substituting this back into the equation for $\mu_i$ will result in the following mean:

$\mu_i = \alpha + \gamma_{W,i} W_i + \beta_S S_i$

$\mu_i = \alpha + (\beta_W + \beta_{WS} S_i) W_i + \beta_S S_i$

And if we expand the bracket to give the conventional form of a continuous interaction, we’ll get the following:

$\mu_i = \alpha + \beta_W W_i + \beta_S S_i + \beta_{WS} W_i S_i$

The new term $ \beta_{WS} W_i S_i $ is now the product of the two predictors. What makes it a continuous interaction, rather than a categorical one, is that instead of estimating a discrete set of slopes for each group, the slope of water now shifts fluidly across the entire observed range of shade and the rate of that shift is governed by $\beta_{WS}$ parameter.

Although we constructed this by making the water slope conditional on shade, the algebra is identical to making the shade slope conditional on water. You cannot have one without the other as the linear interaction is always symmetric. Specifying that the effect of water depends on shade automatically implies that the effect of shade depends on water, with exactly the same interaction coefficient $\beta_{WS}$ governing both.

This brings us to the full model specification:

$B_i \sim \text{Normal}(\mu_i, \sigma)$

$\mu_i = \alpha + \beta_W W_i + \beta_S S_i + \beta_{WS} W_i S_i$

$α \sim \text{Normal}(0.5, 0.25)$

$\beta_w \sim \text{Normal}(0.5, 0.25)$

$\beta_s \sim \text{Normal}(0.5, 0.25)$

$\sigma \sim \text{Exponential}(1)$

The priors on $\alpha$, $\beta_W$, and $\beta_S$ carry over from m_8_4. The challenge is the prior on $\beta_{WS}$, because interaction parameters don’t have a directly interpretable natural scale. The way to reason about it is to ask: What would the strongest plausible interaction look like? If shade at its maximum value ($S_i = 1$) were to completely eliminate the effect of water, then:

\[\gamma_{W,i} = \beta_W + \beta_{WS} \cdot 1 = 0 \implies \beta_{WS} = -\beta_W\]

That’s the ceiling of what the interaction could plausibly do as it would need to exactly cancel the main effect of water. Assigning $\beta_{WS}$ the same prior standard deviation as $\beta_W$, namely $\text{Normal}(0, 0.25)$, places that extreme scenario at roughly two standard deviations out. It’s weakly informative as it doesn’t rule out strong interactions, but it doesn’t treat them as the default either.

The interaction parameter $β_{WS}$ gets the same prior, motivated by the argument that the strongest conceivable interaction is one in which shade at its maximum completely cancels the effect of water, which would require $β_{WS} ≈ −β_W$.

The continuous interaction term $β_WS$, $W_i,$ $S_i$ can be read as: The slope of water on blooms changes by $β_WS$ for every one-unit increase in shade. Equivalently (by the symmetry of interactions), the slope of shade on blooms changes by $β_WS$ for every one-unit increase in water.

with pm.Model() as m_8_5:
    a = pm.Normal("a", 0.5, 0.25)
    bw = pm.Normal("bw", 0, 0.25)
    bs = pm.Normal("bs", 0, 0.25)
    bws = pm.Normal("bws", 0, 0.25)

    mu = (
        a
        + bw * d["water_cent"].values
        + bs * d["shade_cent"].values
        + bws * d["water_cent"].values * d["shade_cent"].values
    )
    sigma = pm.Exponential("sigma", 1)

    blooms_std = pm.Normal("blooms_std", mu, sigma, observed=d["blooms_std"].values)

    m_8_5_trace = pm.sample()
    m_8_5_post = az.extract_dataset(m_8_5_trace.posterior)

8.3.3. Plotting posterior predictions.#

With a triptych plot, the strategy here is to fix shade at each of its three values (−1, 0, +1) and, in each panel, plot the predicted bloom as a function of water.

For m_8_4 where there was no interaction, the slope of water on blooms is identical across all three shade panels. Without an interaction, the model has no mechanism to let one variable’s slope depend on the other so the lines in all three panels are then parallel but just shifted vertically.

For m_8_5 with the interaction, the slope of water on bloom visibly flattens as shade increases (i.e., as light decreases). At low shade (high light), water has a steep positive effect. At high shade (low light), the water effect nearly disappears. This matches the biological reality where light is the binding constraint at high shade levels so really, more water doesn’t help.

The triptych format makes this conditional dependence immediately readable in a way that a table of four parameters never could.

Figure 8.7. Triptych plots of posterior predicted blooms across water and shade treatments.#

Code 8.25#

_, axs = plt.subplots(2, 3, figsize=[9, 5], sharey=True, sharex=True, constrained_layout=True)

n_lines = 20
pred_x = np.array([-1, 1])

for i, shade in enumerate([-1, 0, 1]):
    ind = d.shade_cent == shade
    for ax in axs[:, i]:
        ax.scatter(d.loc[ind, "water_cent"], d.loc[ind, "blooms_std"])
    # top row, m_8_4
    ax = axs[0, i]
    ax.set_title(f"m_8_4 post: shade = {shade:.0f}", fontsize=11)
    pred_y = (
        m_8_4_post["a"][:n_lines].to_numpy()
        + m_8_4_post["bw"][:n_lines].to_numpy() * pred_x.reshape(-1, 1)
        + m_8_4_post["bs"][:n_lines].to_numpy() * shade
    )
    ax.plot(pred_x, pred_y, lw=1, color=(0, 0, 0, 0.4))

    # bottom row, m_8_5
    ax = axs[1, i]
    ax.set_title(f"m_8_5 post: shade = {shade:.0f}", fontsize=11)
    pred_y = (
        m_8_5_post["a"][:n_lines].to_numpy()
        + m_8_5_post["bw"][:n_lines].to_numpy() * pred_x.reshape(-1, 1)
        + m_8_5_post["bs"][:n_lines].to_numpy() * shade
        + m_8_5_post["bws"][:n_lines].to_numpy() * pred_x.reshape(-1, 1) * shade
    )
    ax.plot(pred_x, pred_y, lw=1, color=(0, 0, 0, 0.4))

for ax in axs.flat:
    if ax.get_subplotspec().is_first_col():
        ax.set_ylabel("blooms")
    if ax.get_subplotspec().is_last_row():
        ax.set_xlabel("water")


#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=0.45,
    y=-0.06,
    t="Figure 8.7. Triptych plots of posterior predicted blooms across water and shade treatments. \n \
    Top row: Without an interaction between water and shade. \n \
    Bottom row: With an interaction between water and shade. \n \
    Each plot shows 20 posterior lines for each level of shade.",
    ma="left"
  );

8.3.4 Plotting Prior Predictions#

Applying the same triptych structure to prior predictive draws (before seeing any data) serves as a sanity check on the priors. Since the lines are scattered in the weakly-informative prior, it’s hard to make any sense of what’s going on here so instead we’ll follow the changes for one of the slopes which has been bolded in Figure 8.8.

For the no-interaction model, bolded lines from the same prior draw have the same slope in all three panels which correctly reflects that without an interaction, the water slope cannot vary with shade. For the interaction model, the same bolded lines show systematically changing slopes across panels. Overall, the priors are weakly informative as they stay mostly within the valid outcome space and impose no strong directional bias, which is a reasonable starting point when prior scientific knowledge is limited.

Code 8.26#

with m_8_4:
    m_8_4_priors = pm.sample_prior_predictive(var_names=["a", "bw", "bs"])
    m_8_4_priors = az.extract_dataset(m_8_4_priors.prior)

with m_8_5:
    m_8_5_priors = pm.sample_prior_predictive(var_names=["a", "bw", "bs", "bws"])
    m_8_5_priors = az.extract_dataset(m_8_5_priors.prior)

_, axs = plt.subplots(2, 3, figsize=[9, 5], sharey=True, sharex=True, constrained_layout=True)

n_lines = 20
pred_x = np.array([-1, 1])

for i, shade in enumerate([-1, 0, 1]):
    # top row, m_8_4
    ax = axs[0, i]
    ax.set_title(f"m8.4 prior: shade = {shade:.0f}", fontsize=11)
    pred_y = (
        m_8_4_priors["a"][:n_lines].to_numpy()
        + m_8_4_priors["bw"][:n_lines].to_numpy() * pred_x.reshape(-1, 1)
        + m_8_4_priors["bs"][:n_lines].to_numpy() * shade
    )
    ax.plot(pred_x, pred_y, lw=1, color=(0, 0, 0, 0.4))
    ax.plot(pred_x, pred_y[:, 0], lw=2, color="k")

    # bottom row, m_8_5
    ax = axs[1, i]
    ax.set_title(f"m8.5 prior: shade = {shade:.0f}", fontsize=11)
    pred_y = (
        m_8_5_priors["a"][:n_lines].to_numpy()
        + m_8_5_priors["bw"][:n_lines].to_numpy() * pred_x.reshape(-1, 1)
        + m_8_5_priors["bs"][:n_lines].to_numpy() * shade
        + m_8_5_priors["bws"][:n_lines].to_numpy() * pred_x.reshape(-1, 1) * shade
    )
    ax.plot(pred_x, pred_y, lw=1, color=(0, 0, 0, 0.4))
    ax.plot(pred_x, pred_y[:, 0], lw=2, color="k")

for ax in axs.flat:
    ax.set_ylim(-0.5, 1.5)
    ax.axhline(1, ls="dashed", color=(0, 0, 0, 0.6))
    ax.axhline(0, ls="dashed", color=(0, 0, 0, 0.6))
    if ax.get_subplotspec().is_first_col():
        ax.set_ylabel("blooms")
    if ax.get_subplotspec().is_last_row():
        ax.set_xlabel("water")

#####################
### CODE ADDITION ###
#####################
plt.suptitle(
    x=0.45,
    y=-0.06,
    t="Figure 8.8. Triptych plots of prior predicted blooms across water and and shade treatments. \n \
    Top row: Without an interaction between water and shade. \n \
    Bottom row: With an interaction between water and shade. \n \
    Each plot shows 20 prior lines for each level of shade.",
    ma="left"
  );

Section 8.4 - Summary#

The central theme of this chapter is conditioning which is arguably the most important principle in statistical inference. Every inference you make is conditional on something, whether you’re aware of it or not, and appears in cases like:
- The data being conditional on how it entered our dataset in the first place;
- The posterior distributions which is conditional on the data itself;
- And any model-based inference is conditional on the model its drawing its inference from.
With the example of the Florida manatee, the data about it’s mortality were systematically misleading because they were conditional on survival and determined whether they entered our dataset or not. The manatees that didn’t survive look fundamentally different from those that did but were structurally absent from the samples of important studies on this topic.
Simple linear models from previous chapters assume that each predictor has an independent and fixed association with the outcome. However in this chapter, we’ll extend that framework to allow the association between one predictor and the outcome to be conditional on another predictor, which is known as an interaction (or moderation).
In the example where we regressed log GDP per capita against terrain ruggedness, we found that the relationship produced a negative slope for countries outside Africa which we can intuit by telling ourselves a story of how rugged terrain hinders transport and market access. However when we regressed the same variables for countries within Africa, we instead find they have a positive slope/relationship which historians and anthropologists have since attributed to rugged terrain historically protecting African regions from the slave trade. The same predictor having the opposite effects on a target variable depending on continent is a classic example here of interaction.
For practicing data scientists, the questions we might ask ourselves to perfectly quantify this relationship is why not just split the data by continent and fit two separate models? Well there’s four reasons why it might be a bad idea:
1. Parameters like standard deviation ($\sigma$) don’t depend on continent so we risk losing precision by estimating on smaller datasets;
2. We can’t make probability statements about the continent variable unless it’s explicitly inside the model;
3. Model comparison techniques like WAIC or PSIS requires all models to use the same data;
4. Other methods like Multilevel models (which we touch on in Chapter 13) can borrow information across groups so if we instead split the data, this would make these class of models impossible to use.
One way we can work with interactions is by making use of an indicator variable for both the slope and y-intercept variables to truly account for the “continental conditioning” and differences between African vs non-African countries:

$\mu_i = \alpha_{\text{cid}[i]} + \beta_{\text{cid}[i]}(r_i - \bar{r})$

Where:

$\text{cid}[i]$ is the index variable taking a value of 1 for African nations and 2 for non-African ones.

With both the intercept and slope conditioned by continent, the posterior slopes confirm the reversal: Within Africa, the slope is $\beta_1 \approx +0.13$ while non-African countries have a slope of $\beta_2 \approx -0.14$ which matches the intuition that social scientists developed about ruggedness having a negative effect on GDP, asides from African countries. Model m_8_3 captures more than 95% of PSIS model weight across all three models.
A subtle but important point about interactions are that they are symmetric, meaning the ruggedness-GDP can pose two equally valid conclusions from the model: One being that the slope is conditional on continent while the other valid conclusion is that the intercept is conditional on ruggedness. Both interpretations are simultaneously true and there is no logical basis inside the model for preferring one over the other.
The practical takeaway here is that causal interpretation of an interaction should be driven by what is scientifically manipulable, not by how the math happens to be framed. In this case, continent is a proxy for historical and political factors which are the real causal levers, not geography itself.
The Africa vs non-Africa distinction we made is an example of a categorical interactions where we could estimate a slope per discrete group to keep the interpretation of model outputs manageable. Continuous interactions on the other hand are far more complex because the slope now varies fluidly across the entire range of a second continuous variable. No single number summarises the effect so a table of posterior means and standard deviations becomes nearly useless for understanding what the model is saying.
With the tulip dataset which documents the bloom sizes grown under three levels of soil moisture ($W$ = water, centered: −1 to 1) and three levels of light exposure ($S$ = shade, centered: −1 to 1), we can build a continuous interaction model where we can have the slope of water conditional on shade by giving it it’s own linear model ($\gamma_{W,i}$) within the linear model:

$\gamma_{W,i} = \beta_W + \beta_{WS} S_i$

The science on this is clear that water and light are complementary inputs to photosynthesis so neither variable’s effect should be interpreted as independent of one another. Our m_8_5 model uses the following mathematical expression where the product term $\beta_{WS} W_i S_i$ is the full continuous interaction, the $\beta_{WS}$ parameter governs how the slope of water shifts across the entire observed range of shade (or vice versa), $ \beta_W $ is the effect of water when shade is at its mean, and $\beta_S$ is the effect of shade when water is at its mean:

$\mu_i = \alpha + \gamma_{W,i} W_i + \beta_S S_i $

$\mu_i = \alpha + (\beta_W + \beta_{WS} S_i) W_i + \beta_S S_i$

$\mu_i = \alpha + \beta_W W_i + \beta_S S_i + \beta_{WS} W_i S_i$

Just as in the ruggedness example, this continuous interaction is symmetric in that there’s two reasonable ways to interpret $\beta_{WS} W_i S_i$. We can say that the slope of water on bloom changes by $\beta_{WS}$ for every one-unit increase in shade OR equivalently (by the symmetry of interactions), that the slope of shade on bloom changes by $\beta_{WS}$ for every one-unit increase in water.
Because no single figure can show how a continuous interaction behaves, the solution is to fix shade at three representative values (−1, 0, +1) and plot predicted blooms against water separately through a triptych plot:
- In m_8_4 (no interaction): the slope of water is identical across all three panels as the lines are parallel, just shifted vertically by shade level.
- In m_8_5 (with interaction): the water slope visibly flattens as shade increases. At low shade (high light), water has a steep positive effect. At high shade (low light), the water slope nearly disappears, thus matching the biological intuition that light becomes the binding constraint when it is scarce.
- Note that the number of panels is not fixed at three, we can use as many representative values as the problem demands. But the principle where we hold one variable fixed, vary the other, and repeat should always remain the same.

https://www.telegraph.co.uk/multimedia/archive/02807/Manatee-and-Girl-9_2807045k.jpg

Source

%load_ext watermark
%watermark -n -u -v -iv -w

Chapter 8 - Conditional Manatees

Contents

Chapter 8 - Conditional Manatees#

Key Takeaway#

Intro#

8.1 Building an Interaction#

Code 8.1#

Figure 8.2. Separate linear regressions inside and outside of Africa, for log-GDP against terrain ruggedness.#

8.1.1 Building the Rugged Model#

Code 8.2#

Code 8.3#

Code 8.4#

Code 8.5#

Figure 8.3. Simulating in search of reasonable priors for terrain ruggedness.#

Code 8.6#

8.1.2 Adding an Indicator Variable Isn’t Enough#

Code 8.7#

Code 8.8#

Reviewing WAIC#

Code 8.9#

Code 8.10#

Code 8.11#

Figure 8.4. Including an indicator for African nations has not effect on the slope.#

Code 8.12#

8.1.3 Adding an Interaction Does Work#

Code 8.13#

Code 8.14#

Code 8.15#

Code 8.16#

8.1.4 Plotting the Interaction#

Code 8.17#

Figure 8.5. Posterior predictions for the terrain ruggedness model#

8.2. Symmetry of interactions.#

Figure 8.6. The other side of the interaction between ruggedness and continent.#

Code 8.18#

8.3 Continuous Interactions#

8.3.1 A Winter Flower#

Code 8.19#

8.3.2. The models.#

Code 8.20#

Code 8.21#

Code 8.22#

Code 8.23#

Code 8.24#

8.3.3. Plotting posterior predictions.#

Figure 8.7. Triptych plots of posterior predicted blooms across water and shade treatments.#

Code 8.25#

8.3.4 Plotting Prior Predictions#

Code 8.26#

Section 8.4 - Summary#