12. Transformations and regression

Sometimes its best to transform the data from its raw form before applying regression. Earlier looked at regression to the mean.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import arviz as az
import pymc as pm 
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.robust.scale import mad

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/__init__.py:50: FutureWarning: 
ArviZ is undergoing a major refactor to improve flexibility and extensibility while maintaining a user-friendly interface.
Some upcoming changes may be backward incompatible.
For details and migration guidance, visit: https://python.arviz.org/en/latest/user_guide/migration_guide.html
  warn(

def fit_and_plot_lm(data, predictors, outcome, add_constant=True, show_plot=True, scatter_kws=None, line_kws=None):
    """
    Fit a linear model using statsmodels, print summary, plot, and show formula.
    Args:
        data: pandas DataFrame
        predictors: list of predictor column names (str)
        outcome: outcome column name (str)
        add_constant: whether to add intercept (default True)
        show_plot: whether to plot (default True)
        scatter_kws: dict, kwargs for scatterplot
        line_kws: dict, kwargs for regression line
    """
    X = data[predictors].copy()
    if add_constant:
        X = sm.add_constant(X, prepend=False)
    y = data[outcome]
    model = sm.OLS(y, X)
    results = model.fit()
    print(results.summary())
    # Print formula
    params = results.params
    formula = f"{outcome} = " + " + ".join([f"{params[name]:.2f}*{name}" for name in predictors])
    if add_constant:
        formula = f"{outcome} = {params['const']:.2f} + " + " + ".join([f"{params[name]:.2f}*{name}" for name in predictors])
    print("Formula:", formula)
    # Print residual standard deviation and its uncertainty
    sigma = np.sqrt(results.mse_resid)
    sigma_se = sigma / np.sqrt(2 * results.df_resid)
    print(f"Residual std dev (σ): {sigma:.2f} ± {sigma_se:.2f}")
    # Print median absolute deviation of residuals
    print("MAD of residuals:", round(mad(results.resid), 2))
    # Plot if only one predictor
    if show_plot and len(predictors) == 1:
        x_name = predictors[0]
        ax = sns.scatterplot(data=data, x=x_name, y=outcome, **(scatter_kws or {}))
        x_vals = np.linspace(data[x_name].min(), data[x_name].max(), 100)
        y_vals = params['const'] + params[x_name] * x_vals if add_constant else params[x_name] * x_vals
        ax.plot(x_vals, y_vals, color='red', **(line_kws or {}))
        ax.set_title('Linear Regression Fit')
        plt.show()

def fit_and_plot_bayes(data, predictors, outcome,
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=True,
                       show_posterior=True, show_regression=True,
                       show_residuals=False,
                       n_regression_lines=100):
    """
    Fit a Bayesian linear regression with PyMC and produce diagnostic plots.

    The model is:
        y ~ Normal(mu, sigma)
        mu = intercept + sum_k (slope_k * x_k)
        intercept ~ Normal(intercept_mu, intercept_sigma)
        slope_k   ~ Normal(slope_mu, slope_sigma)     for each predictor k
        sigma     ~ HalfNormal(sigma_sigma)

    Supports one or multiple predictors. After sampling, prints a posterior
    summary and the mean regression formula, then optionally shows trace,
    forest, posterior density, regression-line, and residual plots.

    Args:
        data (pd.DataFrame): DataFrame containing predictor and outcome columns.
        predictors (str | list[str]): Predictor column name(s). A single string
            is treated as a one-predictor model; a list allows multiple predictors.
        outcome (str): Outcome column name.
        intercept_mu (float): Prior mean for the intercept Normal prior.
        intercept_sigma (float): Prior std dev for the intercept Normal prior.
        slope_mu (float): Prior mean for each slope's Normal prior.
        slope_sigma (float): Prior std dev for each slope's Normal prior.
        sigma_sigma (float): Scale for the HalfNormal prior on residual noise.
        samples (int): Number of posterior draws per chain.
        tune (int): Number of NUTS tuning (warm-up) steps.
        hdi_prob (float): Probability mass used for HDI summaries and plots.
        show_trace (bool): If True, plot MCMC traces and marginal densities.
        show_forest (bool): If True, plot a forest plot of posterior means + HDIs.
        show_posterior (bool): If True, plot posterior density for each parameter.
        show_regression (bool): If True, plot data with overlaid posterior
            regression lines (one subplot per predictor; other predictors held at mean).
        show_residuals (bool): If True, plot residuals (y - y_hat using posterior
            mean coefficients) vs fitted values and vs each predictor.
        n_regression_lines (int): Number of posterior draws to overlay on the
            regression plot when show_regression is True.

    Returns:
        arviz.InferenceData: The PyMC trace (posterior samples and sampling stats).
    """
    if isinstance(predictors, str):  # Allow a single predictor to be passed as a bare string.
        predictors = [predictors]  # Normalise to a list so downstream code can iterate uniformly.
    y = data[outcome].values  # Extract outcome values as a NumPy array for PyMC.

    with pm.Model() as model:  # Open a PyMC model context; variables below are registered to it.
        intercept = pm.Normal("intercept", mu=intercept_mu, sigma=intercept_sigma)  # Normal prior on the intercept.
        slopes = []  # Collect slope RVs (not strictly needed, but handy if you want to inspect them).
        mu = intercept  # Start building the linear predictor; will accumulate slope*x terms below.
        for pred in predictors:  # Loop over each predictor column name.
            s = pm.Normal(f"slope_{pred}", mu=slope_mu, sigma=slope_sigma)  # Normal prior on this predictor's slope.
            slopes.append(s)  # Keep a reference to the slope RV.
            mu = mu + s * data[pred].values  # Add this predictor's contribution to the linear predictor.
        sigma = pm.HalfNormal("sigma", sigma=sigma_sigma)  # HalfNormal prior on residual std dev (must be positive).
        likelihood = pm.Normal("y", mu=mu, sigma=sigma, observed=y)  # Likelihood: observed y conditional on mu and sigma.
        trace = pm.sample(samples, tune=tune, idata_kwargs={"log_likelihood": True})  # Run NUTS; store log-likelihood for LOO.

    summary = pm.summary(trace, hdi_prob=hdi_prob)  # Compute posterior summary stats (mean, sd, HDI, r_hat, ess).
    print(summary)  # Print the summary table for the user.

    # Build and print the mean-posterior regression formula for quick inspection.
    posterior = trace.posterior  # Shortcut to the posterior group of the InferenceData.
    intercept_mean = posterior["intercept"].values.flatten().mean()  # Posterior mean of the intercept (flatten chains+draws).
    formula = f"{outcome} = {intercept_mean:.2f}"  # Start the formula string with the intercept.
    for pred in predictors:  # Append a "+ slope*predictor" term for each predictor.
        slope_mean = posterior[f"slope_{pred}"].values.flatten().mean()  # Posterior mean slope for this predictor.
        formula += f" + {slope_mean:.2f}*{pred}"  # Append this term to the formula string.
    print(f"\nRegression formula: {formula}")  # Print the assembled formula.

    sigma_draws = posterior["sigma"].values.flatten()
    print(f"Residual std dev (σ): {sigma_draws.mean():.2f} ± {sigma_draws.std():.2f}")

    intercept_draws = posterior["intercept"].values.flatten()
    X_mat = data[predictors].values
    slope_mat = np.column_stack([posterior[f"slope_{pred}"].values.flatten() for pred in predictors])
    y_hat_all = intercept_draws[:, None] + slope_mat @ X_mat.T
    bayes_r2 = az.r2_score(y, y_hat_all)
    print(f"Bayesian R²: {bayes_r2['r2']:.3f} ± {bayes_r2['r2_std']:.3f}")

    loo = az.loo(trace, pointwise=True)  # PSIS-LOO; pointwise=True gives per-obs Pareto-k diagnostics.
    n_obs = len(y)
    print(f"LOO-ELPD: {loo.elpd_loo:.2f} ± {loo.se:.2f}  (p_loo={loo.p_loo:.1f}) <- should be lower than number of parameters ({len(predictors) + 2}) for reliable estimates")
    print("When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.")
    print(f"LOO log score (per obs): {loo.elpd_loo / n_obs:.3f} ± {loo.se / n_obs:.3f}")
    n_bad_k = int((loo.pareto_k.values > 0.7).sum())
    if n_bad_k > 0:
        print(f"  Warning: {n_bad_k} observations with Pareto-k > 0.7 (unreliable LOO estimates)")
    if loo.p_loo >= len(predictors) + 2:
        print(f"  *** WARNING: p_loo ({loo.p_loo:.1f}) >= number of parameters ({len(predictors) + 2}) — LOO estimates may be unreliable ***")

    # LOO R²: importance-sample the leave-one-out predicted means from the posterior draws
    log_like = trace.log_likelihood["y"].values       # (chains, draws, obs)
    ll_2d = log_like.reshape(-1, n_obs)               # (total_draws, obs)
    log_w = -ll_2d                                    # log IS weight = -log p(y_i | theta)
    log_w -= log_w.max(axis=0, keepdims=True)         # numerical stability before exp
    w = np.exp(log_w)
    w /= w.sum(axis=0, keepdims=True)                 # normalize per observation
    y_hat_loo = (w * y_hat_all).sum(axis=0)           # weighted LOO predicted means, shape (obs,)
    loo_r2 = 1 - np.sum((y - y_hat_loo) ** 2) / np.sum((y - y.mean()) ** 2)
    print(f"LOO R²: {loo_r2:.3f}")
    if loo_r2 < bayes_r2['r2']:
        print(f"  *** WARNING: LOO R² ({loo_r2:.3f}) < Bayesian R² ({bayes_r2['r2']:.3f}) — model may be overfitting ***")

    if show_trace:  # Optional: MCMC trace diagnostics.
        az.plot_trace(trace)  # Plot chain traces and marginal densities per parameter.
        plt.tight_layout()  # Tidy up subplot spacing.
        plt.show()  # Render the figure.

    if show_forest:  # Optional: forest plot of parameter posteriors.
        az.plot_forest(trace, hdi_prob=hdi_prob)  # Show posterior means and HDIs for all parameters.
        plt.show()  # Render the figure.

    if show_posterior:  # Optional: posterior density plots.
        az.plot_posterior(trace, hdi_prob=hdi_prob)  # Density plot per parameter, annotated with HDI.
        plt.show()  # Render the figure.

    if show_regression:  # Optional: overlay posterior regression lines on data.
        a_samples = posterior["intercept"].values.flatten()  # Flattened posterior draws of the intercept.
        slope_samples = {pred: posterior[f"slope_{pred}"].values.flatten() for pred in predictors}  # Flattened slope draws per predictor.
        idx = np.random.choice(len(a_samples), n_regression_lines, replace=False)  # Random subset of draw indices to plot.

        fig, axes = plt.subplots(1, len(predictors), figsize=(6 * len(predictors), 5))  # One subplot per predictor.
        if len(predictors) == 1:  # plt.subplots returns a single Axes when ncols=1; wrap it for uniform handling.
            axes = [axes]  # Make axes iterable in the single-predictor case.

        for ax, pred in zip(axes, predictors):  # Plot each predictor's partial regression view.
            x = data[pred].values  # Observed values of this predictor.
            ax.scatter(x, y, alpha=0.5)  # Scatter of y vs this predictor.
            x_grid = np.linspace(x.min(), x.max(), 100)  # Dense grid across this predictor for smooth lines.

            # For each posterior draw, compute the contribution of *other* predictors
            # held at their sample means so the plotted line isolates this predictor's effect.
            other_contribution = np.zeros(len(a_samples))  # Per-draw constant offset from held-fixed predictors.
            for other_pred in predictors:  # Sum contributions across all other predictors.
                if other_pred != pred:  # Skip the predictor currently on the x-axis.
                    other_contribution += slope_samples[other_pred] * data[other_pred].mean()  # slope_draw * mean(x_other).

            for i in idx:  # Overlay a thin gray line per sampled posterior draw.
                y_line = a_samples[i] + other_contribution[i] + slope_samples[pred][i] * x_grid  # Predicted y on the grid.
                ax.plot(x_grid, y_line, alpha=0.05, color="gray")  # Low alpha so overlap shades posterior uncertainty.

            # Posterior mean line for emphasis.
            mean_other = sum(slope_samples[op].mean() * data[op].mean() for op in predictors if op != pred)  # Mean offset from other predictors.
            y_mean = a_samples.mean() + mean_other + slope_samples[pred].mean() * x_grid  # Mean posterior prediction on the grid.
            ax.plot(x_grid, y_mean, color="red")  # Plot the mean line in red.
            ax.set_xlabel(pred)  # Label x-axis with the predictor name.
            ax.set_ylabel(outcome)  # Label y-axis with the outcome name.
            ax.set_title(f"{outcome} vs {pred} (others at mean)")  # Title notes that other predictors are held at mean.

        plt.tight_layout()  # Tidy spacing across subplots.
        plt.show()  # Render the figure.

    if show_residuals:  # Optional: residual diagnostics using posterior mean coefficients.
        a_mean = posterior["intercept"].values.flatten().mean()  # Posterior mean intercept (point estimate).
        slope_means = {pred: posterior[f"slope_{pred}"].values.flatten().mean() for pred in predictors}  # Posterior mean slope per predictor.
        y_hat = a_mean + sum(slope_means[pred] * data[pred].values for pred in predictors)  # Fitted values using mean coefficients.
        residuals = y - y_hat  # Residuals: observed minus fitted.

        fig, axes = plt.subplots(1, len(predictors) + 1, figsize=(6 * (len(predictors) + 1), 5))  # One subplot for fitted + one per predictor.

        axes[0].scatter(y_hat, residuals, alpha=0.5)  # Residuals vs fitted values — should look like noise around 0.
        axes[0].axhline(0, color="red", linestyle="--")  # Zero reference line to judge bias/structure.
        axes[0].set_xlabel("Fitted values")  # x-axis label.
        axes[0].set_ylabel("Residuals")  # y-axis label.
        axes[0].set_title("Residuals vs Fitted")  # Subplot title.

        for ax, pred in zip(axes[1:], predictors):  # Residuals vs each predictor — check for unmodeled structure.
            ax.scatter(data[pred].values, residuals, alpha=0.5)  # Scatter residuals against this predictor.
            ax.axhline(0, color="red", linestyle="--")  # Zero reference line.
            ax.set_xlabel(pred)  # x-axis label is the predictor name.
            ax.set_ylabel("Residuals")  # y-axis label.
            ax.set_title(f"Residuals vs {pred}")  # Subplot title.

        plt.tight_layout()  # Tidy subplot spacing.
        plt.show()  # Render the figure.

    return trace  # Return the InferenceData so the caller can do further analysis.

12.1 Linear transformations

Scaling of predictors and regression coefficients

The coeffient $\beta_j$ represents the average difference in y, comparing items that differ by one unit on the jth predictor, but are the same on all other predictors. In some cases, a difference of 1 unit in x is not the most relevant comparison.

Earlier predicting earnings from height:

Earnings = -85000 + 1600*height + error, with a residual standard error of 22,000.

The intercept of -85000 is not meaningful, as it corresponds to the predicted earnings of a person with a height of 0 inches. The slope of 1600 means that for every additional inch in height, earnings are predicted to increase by $1600.

Height can be expressed in many different units.

To understand the coefficients, we should better understand the variation in the predictor. One approach is to scale by the standard deviation of the predictor (3.8 inches for height).

Linear transformations of predictors X or outcome y dont affect the fit of the model. They help with

earnings = pd.read_csv('../ros_data/earnings.csv', skiprows=0)
# earnings['earn_k'] = earnings['earn'] / 1000
# earnings['c_height'] = earnings['height'] - 66 # Center height around 66 inches for better interpretability of the intercept

# calculate standard deviation of height
height_std = earnings['height'].std()
print(f"Standard deviation of height: {height_std:.2f} inches")

display(earnings.head())

earnings_clean = earnings.dropna(subset=['height', 'earn'])

earnings_model = fit_and_plot_bayes(earnings_clean, 'height', 'earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

print(earnings_model)

Standard deviation of height: 3.83 inches

	height	weight	male	earn	earnk	ethnicity	education	mother_education	father_education	walk	exercise	smokenow	tense	angry	age
0	74	210.0	1	50000.0	50.0	White	16.0	16.0	16.0	3	3	2.0	0.0	0.0	45
1	66	125.0	0	60000.0	60.0	White	16.0	16.0	16.0	6	5	1.0	0.0	0.0	58
2	64	126.0	0	30000.0	30.0	White	16.0	16.0	16.0	8	1	2.0	1.0	1.0	29
3	65	200.0	0	25000.0	25.0	White	17.0	17.0	NaN	8	1	2.0	0.0	0.0	57
4	63	110.0	0	50000.0	50.0	Other	16.0	16.0	16.0	5	6	2.0	0.0	0.0	91

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_height, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                  mean      sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept      -23.647  50.330  -120.873     76.484      0.507    0.543   
slope_height   321.562   2.501   316.625    326.353      0.025    0.026   
sigma         6718.789  24.245  6669.107   6764.250      0.239    0.280   

              ess_bulk  ess_tail  r_hat  
intercept       9837.0    6393.0    1.0  
slope_height   10383.0    6284.0    1.0  
sigma          10320.0    5888.0    1.0  

Regression formula: earn = -23.65 + 321.56*height
Residual std dev (σ): 6718.79 ± 24.24
Bayesian R²: 0.003 ± 0.000

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

LOO-ELPD: -27662.59 ± 1769.98  (p_loo=106.4) <- should be lower than number of parameters (3) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -15.233 ± 0.975
  Warning: 2 observations with Pareto-k > 0.7 (unreliable LOO estimates)

Inference data with groups:
	> posterior
	> log_likelihood
	> sample_stats
	> observed_data

Standardisation using z-scores

We can scale the coefficients by standardising the predictors. This is done by subtracting the mean and dividing by the standard deviation. The resulting variable is called a z-score.

For height above, this would be: z_height = (height - 66) / 3.8

Coefficients are then interpreted in units of standard deviations with respect to the corresponding predictor.

Why this is useful:

1. Coefficient now tells us 'what happens if height increases by one standard deviation?'. One SD is roughly the size of a typical difference between an average person and a randomly selected person.
2. Intercept is interpretable. With z-scores, the intercept corresponds to the predicted earnings of a person when all predictors are at their mean.
3. We can compare predictors: If we have two predictors, one with a coefficient of 0.5 and another with a coefficient of 1.0, we can say that the second predictor has a stronger relationship with the outcome variable than the first predictor, since a one standard deviation increase in the second predictor is associated with a larger change in the outcome variable compared to a one standard deviation increase in the first predictor. Cant do this with raw data.

Some statisticians prefer dividing by 2 standard deviations - this makes continuous predictors more comparable to binary predictors.

Standardisation is best done using our data if we have a large sample. If we have a small sample, we can use the standard deviation from a larger, more representative dataset.

12.2 Centering and standardizing for models with interactions

Similar challenges arise with models that have interaction terms. The interpretation of the main effects becomes more difficult when there are interactions, as the main effects represent the effect of a predictor when the other predictor is at its reference level (often 0).

kid_score = -7.05 + 45.94mom_hs + 0.92mom_iq + -0.43*mom_hs_iq

Coefficient on mom_hs is 45.94, which means that the predicted kid_score is 45.94 points higher for a child whose mother has a high school degree compared to a child whose mother does not have a high school degree, when the mother's IQ is 0. This is not a meaningful comparison, as an IQ of 0 is not realistic.

Centering by subtracting the mean

We can simplify interpretation by centering the predictors. This is done by subtracting the mean from each predictor.

Using a conventional centering point

We can also use an understandable value as the centering point. E.g., midpoint of the range of the predictor, or a meaningful value.

Standardizing by subtracting the mean and dividing by 2 standard deviations

Centering helped with interpretation, but there is still a scaling problem. Coefficients are not comparable across predictors. Mom completing high school is a complete change to a single point in IQ, which is not much of a change in IQ.

Why scaling by 2 standard deviations?

So that we can compare the coefficients of continuous predictors to binary predictors.

If we have a binary variable that is 0 50% of the time, its SD is 0.5. When standardising by 1SD, variable becomes -1 to 1, a 2 unit change so coefficient is half the difference between the two groups. When standardising by 2SD, variable becomes -0.5 to 0.5, a 1 unit change so coefficient is the difference between the two groups.

Multiplying each regression coefficient by 2 standard deviations of its predictor

For models with no interactions, we can simply multiply each regression coefficient by 2 standard deviations of its predictor.

kidiq = pd.read_csv('../ros_data/kidiq.csv', skiprows=0)
display(kidiq.head())

# interaction term between mom_hs and mom_iq
kidiq['mom_hs_iq'] = kidiq['mom_hs'] * kidiq['mom_iq']

fit1 = fit_and_plot_bayes(kidiq, ['mom_hs', 'mom_iq', 'mom_hs_iq'], 'kid_score',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=False, show_forest=False,
                       show_posterior=False, show_regression=False,
                       n_regression_lines=100)

# Centering by subtracting the mean
kidiq['c_mom_iq'] = kidiq['mom_iq'] - kidiq['mom_iq'].mean()  # Center mom_iq by subtracting its mean
kidiq['c_mom_hs'] = kidiq['mom_hs'] - kidiq['mom_hs'].mean()  # Center mom_hs by subtracting its mean
kidiq['c_mom_hs_iq'] = kidiq['c_mom_hs'] * kidiq['c_mom_iq']  # Interaction of centered predictors

fit2 = fit_and_plot_bayes(kidiq, ['c_mom_hs', 'c_mom_iq', 'c_mom_hs_iq'], 'kid_score',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=False, show_forest=False,
                       show_posterior=False, show_regression=False,
                       n_regression_lines=100)

# Standardising by subtracting the mean and dividing by 2 standard deviations
kidiq['s_mom_iq'] = (kidiq['mom_iq'] - kidiq['mom_iq'].mean()) / (2 * kidiq['mom_iq'].std())  # Standardize mom_iq
kidiq['s_mom_hs'] = (kidiq['mom_hs'] - kidiq['mom_hs'].mean()) / (2 * kidiq['mom_hs'].std())  # Standardize mom_hs
kidiq['s_mom_hs_iq'] = kidiq['s_mom_hs'] * kidiq['s_mom_iq']  # Interaction of standardized predictors

fit3 = fit_and_plot_bayes(kidiq, ['s_mom_hs', 's_mom_iq', 's_mom_hs_iq'], 'kid_score',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=False, show_forest=False,
                       show_posterior=False, show_regression=False,
                       n_regression_lines=100)

	kid_score	mom_hs	mom_iq	mom_work	mom_age
0	65	1	121.117529	4	27
1	98	1	89.361882	4	25
2	85	1	115.443165	4	27
3	83	1	99.449639	3	25
4	115	1	92.745710	4	27

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_mom_hs, slope_mom_iq, slope_mom_hs_iq, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 6 seconds.

                   mean      sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept        -7.154  13.202   -32.833     18.421      0.325    0.207   
slope_mom_hs     46.178  14.660    16.205     73.175      0.361    0.220   
slope_mom_iq      0.922   0.143     0.635      1.188      0.004    0.002   
slope_mom_hs_iq  -0.430   0.155    -0.714     -0.112      0.004    0.002   
sigma            18.035   0.624    16.834     19.234      0.011    0.009   

                 ess_bulk  ess_tail  r_hat  
intercept          1647.0    2375.0   1.01  
slope_mom_hs       1650.0    2532.0   1.01  
slope_mom_iq       1656.0    2481.0   1.01  
slope_mom_hs_iq    1642.0    2405.0   1.01  
sigma              3420.0    3711.0   1.00  

Regression formula: kid_score = -7.15 + 46.18*mom_hs + 0.92*mom_iq + -0.43*mom_hs_iq
Residual std dev (σ): 18.04 ± 0.62
Bayesian R²: 0.228 ± 0.030

Initializing NUTS using jitter+adapt_diag...

LOO-ELPD: -1872.48 ± 14.31  (p_loo=4.7) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -4.314 ± 0.033

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_c_mom_hs, slope_c_mom_iq, slope_c_mom_hs_iq, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                     mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept          87.618  0.909    85.906     89.472      0.010    0.009   
slope_c_mom_hs      2.835  2.465    -2.001      7.667      0.032    0.024   
slope_c_mom_iq      0.588  0.062     0.460      0.701      0.001    0.001   
slope_c_mom_hs_iq  -0.484  0.162    -0.797     -0.162      0.002    0.001   
sigma              18.021  0.625    16.807     19.256      0.007    0.007   

                   ess_bulk  ess_tail  r_hat  
intercept            8314.0    6668.0    1.0  
slope_c_mom_hs       5991.0    6361.0    1.0  
slope_c_mom_iq       7249.0    6052.0    1.0  
slope_c_mom_hs_iq    6537.0    6771.0    1.0  
sigma                8810.0    6513.0    1.0  

Regression formula: kid_score = 87.62 + 2.84*c_mom_hs + 0.59*c_mom_iq + -0.48*c_mom_hs_iq
Residual std dev (σ): 18.02 ± 0.62
Bayesian R²: 0.232 ± 0.031

Initializing NUTS using jitter+adapt_diag...

LOO-ELPD: -1872.60 ± 14.35  (p_loo=4.9) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -4.315 ± 0.033

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_s_mom_hs, slope_s_mom_iq, slope_s_mom_hs_iq, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                     mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept          87.619  0.908    85.884     89.402      0.010    0.009   
slope_s_mom_hs      2.376  1.978    -1.596      6.156      0.025    0.019   
slope_s_mom_iq     17.585  1.836    14.013     21.131      0.021    0.018   
slope_s_mom_hs_iq -11.892  4.040   -19.694     -3.913      0.051    0.040   
sigma              18.017  0.624    16.840     19.262      0.006    0.007   

                   ess_bulk  ess_tail  r_hat  
intercept            8692.0    6853.0    1.0  
slope_s_mom_hs       6205.0    6240.0    1.0  
slope_s_mom_iq       7743.0    6458.0    1.0  
slope_s_mom_hs_iq    6284.0    6180.0    1.0  
sigma                9370.0    6291.0    1.0  

Regression formula: kid_score = 87.62 + 2.38*s_mom_hs + 17.59*s_mom_iq + -11.89*s_mom_hs_iq
Residual std dev (σ): 18.02 ± 0.62
Bayesian R²: 0.231 ± 0.031
LOO-ELPD: -1872.62 ± 14.35  (p_loo=5.0) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -4.315 ± 0.033

12.3 Correlation and 'regression to the mean'

For a regression with constant term and one predictor, $y = a + bx + error$, if both x and y are standardised ($x = (x - mean(x)) / sd(x)$ and $y = (y - mean(y)) / sd(y)$), then the intercept a will be 0 and the slope b will be the correlation between x and y. Thus it will be between -1 and 1. The slope of regression with one predictor $b = \rho \sigma_y / \sigma_x$, where $\rho$ is the correlation between x and y, $\sigma_y$ is the standard deviation of y, and $\sigma_x$ is the standard deviation of x.

The principal component line and the regression line

With two standardised variables, the principal component line is the line that minimises the sum of squared distances from the data points to the line, while the regression line is the line that minimises the sum of squared vertical distances from the data points to the line. The principal component line is y = x. The regression line is flatter than the principal component line. Regression line gives unbiased predictions, in that it goes through the average of y given x.

Regression to the mean

When x and y are standardised, slope is always less than 1. So when x is 1SD above the mean, the predicted value of y is between 0 and 1 SD above the mean. This phenomenon is called regression to the mean - y is predicted to be closer to mean than x.

E.g., woman 10 inches taller than average and correlation of mother and daughter height is 0.5, then daughter is predicted to be 5 inches taller than average, which is still taller than average, but closer to the mean than the mother's height.

Similar for not perfectly correlated variables. E.g., we expect a team over two seasons who did very well in the first season to do less well in the second season, and a team that did very badly in the first season to do better in the second season, even if there is no change in the underlying quality of the teams.

Naive interpretation is that variable phenomena become more average over time. This is wrong due to the error, where some points fall closer to mean than others.

12.4 Logarithmic transformations

When linearity and additivity are not reasonable assumptions, we can use transformations of the predictors and/or the outcome variable to achieve a better fit. One common transformation is the logarithmic transformation.

Log transformations force the outcome variable to be positive, and they are multiplicative rather than additive.

Plain linear regression can predict negative values, which is not possible for some outcomes (e.g., earnings). Log transformation can help with this issue. Linear models assume effects add the same amount, but in real life, effects often multiply.

Logging outcome fixes both:

1. If we model log(y), predictions are on log scale, to get predictions on original scale, we need to exponentiate the predictions. Exponentiating a negative number gives a positive number, so we can never get negative predictions.
2. If we model log(y), the effect of a predictor is multiplicative rather than additive. A one unit increase in a predictor is associated with a percentage change in the outcome variable, rather than an absolute change.

$log(y) = b₀ + b₁x₁ + b₂x₂ + ε$

$y = e^{b₀ + b₁x₁ + b₂x₂ + ε}$

$y = e^{b₀} * e^{b₁x₁} * e^{b₂x₂} * e^ε$

Earnings and height

Take log of earnings and regress on height:

log_earn = 5.92 + 0.06*height

exponentiating the coefficients gives:

earnings = e^{5.92} * e^{0.06height} earnings = 372 * 1.06height

1 inch increase in height is associated with a 6% increase in earnings.

earnings = pd.read_csv('../ros_data/earnings.csv', skiprows=0)
# earnings['earn_k'] = earnings['earn'] / 1000
# earnings['c_height'] = earnings['height'] - 66 # Center height around 66 inches for better interpretability of the intercept

# Calculate log of earnings
earnings['log_earn'] = np.log(earnings['earn'])

display(earnings.head())

earnings_clean = earnings[earnings['earn'] > 0].dropna(subset=['height', 'log_earn'])

earnings_model = fit_and_plot_bayes(earnings_clean, 'height', 'log_earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

# --- Posterior predictive check (PPC) ---
# Extract flattened posterior draws (all chains concatenated) for each parameter.
posterior        = earnings_model.posterior
intercept_draws  = posterior["intercept"].values.flatten()
slope_draws      = posterior["slope_height"].values.flatten()
sigma_draws      = posterior["sigma"].values.flatten()

x     = earnings_clean["height"].values   # observed predictor
y_obs = earnings_clean["log_earn"].values  # observed outcome

rng    = np.random.default_rng(42)          # seeded RNG for reproducibility
n_sims = len(intercept_draws)               # total number of posterior draws (chains × samples)
yrep   = np.array([                         # build yrep: one simulated dataset per posterior draw
    rng.normal(                             #   draw y values from Normal(mu, sigma) —
        intercept_draws[i]                  #     mu = intercept (scalar, same for all obs in this draw)
        + slope_draws[i] * x,              #         + slope * height (vector, one value per observation)
        sigma_draws[i]                      #     sigma = residual std dev (scalar, same for all obs)
    )                                       #   result: array of length n_obs for this draw
    for i in range(n_sims)                  #   repeat for every posterior draw
])                                          # final shape: (n_sims, n_obs)

# Randomly select 100 simulated datasets to plot (same as R's `subset <- sample(n_sims, 100)`).
subset = rng.choice(n_sims, 100, replace=False)

# Wrap observed data and the 100 yrep draws into an ArviZ InferenceData object.
# posterior_predictive shape must be (chains, draws, obs); np.newaxis adds the chain dim.
idata_ppc = az.from_dict(
    observed_data={"y": y_obs},
    posterior_predictive={"y": yrep[np.newaxis, subset, :]}
)

# plot_ppc overlays the density of each simulated dataset (thin lines) against
# the observed data density (thick line) — equivalent to bayesplot::ppc_dens_overlay.
az.plot_ppc(idata_ppc, num_pp_samples=100, kind="kde")
plt.title("PPC: log(earn) ~ height")
plt.show()

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/pandas/core/arraylike.py:402: RuntimeWarning: divide by zero encountered in log
  result = getattr(ufunc, method)(*inputs, **kwargs)

	height	weight	male	earn	earnk	ethnicity	education	mother_education	father_education	walk	exercise	smokenow	tense	angry	age	log_earn
0	74	210.0	1	50000.0	50.0	White	16.0	16.0	16.0	3	3	2.0	0.0	0.0	45	10.819778
1	66	125.0	0	60000.0	60.0	White	16.0	16.0	16.0	6	5	1.0	0.0	0.0	58	11.002100
2	64	126.0	0	30000.0	30.0	White	16.0	16.0	16.0	8	1	2.0	1.0	1.0	29	10.308953
3	65	200.0	0	25000.0	25.0	White	17.0	17.0	NaN	8	1	2.0	0.0	0.0	57	10.126631
4	63	110.0	0	50000.0	50.0	Other	16.0	16.0	16.0	5	6	2.0	0.0	0.0	91	10.819778

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_height, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

               mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  ess_bulk  \
intercept     5.898  0.390     5.111      6.636      0.008    0.006    2370.0   
slope_height  0.057  0.006     0.046      0.069      0.000    0.000    2362.0   
sigma         0.878  0.016     0.847      0.908      0.000    0.000    3556.0   

              ess_tail  r_hat  
intercept       2520.0    1.0  
slope_height    2442.0    1.0  
sigma           2886.0    1.0  

Regression formula: log_earn = 5.90 + 0.06*height
Residual std dev (σ): 0.88 ± 0.02
Bayesian R²: 0.060 ± 0.011
LOO-ELPD: -2100.63 ± 38.77  (p_loo=3.9) <- should be lower than number of parameters (3) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -1.290 ± 0.024

Why natural log?

Natural logs (base e) are used because they are directly interpretable as proportional differences. Coefficient of 0.05 ~ 5% change. Sometimes written as ln. Log base 10 can be used, but it is not as interpretable, as a coefficient of 0.05 would correspond to a 12% change in the outcome variable, which is less intuitive than a 5% change.

Building a regression model on the log scale

Adding another predictor to the model: after adding male to the model, the coefficient for height has changed from 0.06 to 0.02, which means that the predicted percentage increase in earnings for a one inch increase in height is now 2% instead of 6%. The coefficient for male is 0.37 - exponentiating gives 45% increase in earnings for males compared to females.

earnings_model = fit_and_plot_bayes(earnings_clean, ['height', 'male'], 'log_earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_height, slope_male, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 9 seconds.

               mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  ess_bulk  \
intercept     7.972  0.508     6.957      8.927      0.011    0.007    2085.0   
slope_height  0.024  0.008     0.009      0.040      0.000    0.000    2060.0   
slope_male    0.372  0.062     0.250      0.493      0.001    0.001    2385.0   
sigma         0.869  0.015     0.838      0.896      0.000    0.000    4487.0   

              ess_tail  r_hat  
intercept       2762.0    1.0  
slope_height    2708.0    1.0  
slope_male      3244.0    1.0  
sigma           4411.0    1.0  

Regression formula: log_earn = 7.97 + 0.02*height + 0.37*male
Residual std dev (σ): 0.87 ± 0.01
Bayesian R²: 0.081 ± 0.012
LOO-ELPD: -2083.54 ± 38.68  (p_loo=4.8) <- should be lower than number of parameters (4) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -1.279 ± 0.024

Naming inputs: here male is a binary variable, so we can interpret the coefficient as the percentage difference in earnings between males and females, holding height constant. Male = 1 for men.

Residual deviation and R2: \sigma is the residual standard deviation on the log scale, 0.87, implying 68% of log earnings are within 0.87 of the predicted value. 70 inch woman has earnings 7.99 + 0.02 * 70 + 0.37 * 0 = 9.37 +- 0.87, which means that her earnings are predicted to be between e^{9.37 - 0.87} and e^{9.37 + 0.87}, which is between $5k and 28k. This is very wide, which means that the model is not very good at predicting earnings. R2 is 0.08, which means that the model explains only 8% of the variance in log earnings.

Including an interaction: between height and sex: log_earn = 8.56 + 0.01 * height - 0.97 * male + 0.02 * height_male. Intercept is when height and male are 0, which is not meaningful and has no real world interpretation. Height coefficient is 0.01, which is the predicted difference in log earnings for one inch difference in height when male is 0. Coefficient for male is predicted difference in log earnings between males and females when height is 0, which is not meaningful. Coefficient for height_male is the difference in slopes of lines predicting log earnings from height, comparing men to women. Thus a difference of 1 inch in height is associated with a 2% increase in earnings for men compared to women. The standard error is large compared to the coefficient, so we are not very confident in this estimate.

# interaction between height and male
earnings_clean['height_male'] = earnings_clean['height'] * earnings_clean['male']

earnings_model = fit_and_plot_bayes(earnings_clean, ['height', 'male', 'height_male'], 'log_earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_height, slope_male, slope_height_male, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 16 seconds.

                    mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept          8.562  0.672     7.268      9.886      0.015    0.008   
slope_height       0.015  0.010    -0.005      0.036      0.000    0.000   
slope_male        -0.967  1.047    -3.127      0.995      0.024    0.014   
slope_height_male  0.020  0.015    -0.009      0.052      0.000    0.000   
sigma              0.868  0.015     0.839      0.899      0.000    0.000   

                   ess_bulk  ess_tail  r_hat  
intercept            2071.0    3294.0    1.0  
slope_height         2073.0    3395.0    1.0  
slope_male           1959.0    3050.0    1.0  
slope_height_male    1922.0    2987.0    1.0  
sigma                4740.0    4720.0    1.0  

Regression formula: log_earn = 8.56 + 0.01*height + -0.97*male + 0.02*height_male
Residual std dev (σ): 0.87 ± 0.02
Bayesian R²: 0.082 ± 0.012
LOO-ELPD: -2083.70 ± 38.78  (p_loo=5.8) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -1.279 ± 0.024

Linear transformation to make coefficients more interpretable: we can rescale height by having a mean of 0 and a standard deviation of 1. This is done by standardising height: z_height = (height - mean(height)) / sd(height). The coefficient for z_height is now the predicted percentage increase in earnings for a one standard deviation increase in height, which is more interpretable than the coefficient for height in inches.

9.55 + 0.06* z_height + 0.36 * male + 0.08 * z_height_male

Intercept is where z_height and male are 0, which is the predicted log earnings for a person of average height and female = 9.55, which is $14k.

Coefficient for z_height is 0.06, which is the predicted difference in log earnings for a one standard deviation increase in height among females.

Coefficient for male is predicted difference in log earnings between males and females if z_height is 0, which is the predicted difference in log earnings between males and females of average height. Log difference of 0.36 corresponds to a 43% increase in earnings for males compared to females of average height.

Coefficient for z_height_male is the difference in slopes of lines predicting log earnings from z_height, comparing men to females. Thus a difference of 1 standard deviation in height is associated with an 8% increase in earnings for men compared to women.

# standardise height by subtracting mean and dividing by std
earnings_clean['z_height'] = (earnings_clean['height'] - earnings_clean['height'].mean()) / earnings_clean['height'].std()
earnings_clean['z_height_male'] = earnings_clean['z_height'] * earnings_clean['male']

earnings_model = fit_and_plot_bayes(earnings_clean, ['z_height', 'male', 'z_height_male'], 'log_earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_z_height, slope_male, slope_z_height_male, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.

                      mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept            9.549  0.037     9.471      9.617      0.001    0.000   
slope_z_height       0.058  0.041    -0.019      0.142      0.001    0.000   
slope_male           0.356  0.062     0.232      0.475      0.001    0.001   
slope_z_height_male  0.074  0.061    -0.047      0.189      0.001    0.001   
sigma                0.869  0.015     0.839      0.899      0.000    0.000   

                     ess_bulk  ess_tail  r_hat  
intercept              4017.0    5331.0    1.0  
slope_z_height         4278.0    5094.0    1.0  
slope_male             5352.0    5671.0    1.0  
slope_z_height_male    4812.0    5198.0    1.0  
sigma                  7316.0    4555.0    1.0  

Regression formula: log_earn = 9.55 + 0.06*z_height + 0.36*male + 0.07*z_height_male
Residual std dev (σ): 0.87 ± 0.02
Bayesian R²: 0.082 ± 0.013
LOO-ELPD: -2083.81 ± 38.73  (p_loo=5.9) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -1.279 ± 0.024

Further difficulties in interpretation

With the model log_earn = 8.00 + 0.02 * height + 0.37 * male

Interpreting the coefficient for height: comparing two adults of same sex, one who is 1 inch taller than the other, the taller adult is predicted to have 2% higher earnings.

Interpreting the coefficient for male: comparing two adults of the same height, but different sex, the male is predicted to have 45% higher earnings than the female. However, the same height means that a woman would be above average height and the male would be below average height, which is not a meaningful comparison.

Log-log model: transforming the input and outcome variables

If we transform both the input and outcome, the coefficient can be interpreted as the proportional difference in x.

log_earn = 2.97 + 1.57 * log_height + 0.37 * male

For each 1% increase in height, earnings are predicted to increase by 1.57%. More precisely, comparing two peoople who differ in height by factor of 1.01 (i.e., 1% increase), this corresponds to a difference of 0.01 in predicted log earnings, which corresponds to a 1.57% increase in earnings.

earnings_clean['log_height'] = np.log(earnings_clean['height'])

earnings_model = fit_and_plot_bayes(earnings_clean, ['log_height', 'male'], 'log_earn',
                       intercept_mu=0, intercept_sigma=50,
                       slope_mu=0, slope_sigma=50,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_height, slope_male, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 35 seconds.

                   mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept         2.830  2.182    -1.503      6.992      0.052    0.031   
slope_log_height  1.604  0.524     0.611      2.655      0.012    0.007   
slope_male        0.372  0.062     0.247      0.490      0.001    0.001   
sigma             0.868  0.015     0.839      0.897      0.000    0.000   

                  ess_bulk  ess_tail  r_hat  
intercept           1788.0    2422.0    1.0  
slope_log_height    1786.0    2399.0    1.0  
slope_male          2357.0    3420.0    1.0  
sigma               3460.0    3918.0    1.0  

Regression formula: log_earn = 2.83 + 1.60*log_height + 0.37*male
Residual std dev (σ): 0.87 ± 0.02
Bayesian R²: 0.081 ± 0.012
LOO-ELPD: -2083.48 ± 38.69  (p_loo=4.9) <- should be lower than number of parameters (4) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -1.279 ± 0.024

Taking logarithms even when not necessary

If variable has a narrow dynamic range (ratio between high and low values is close to 1), then log transformation may not be necessary.

However, log transformation can still be useful for interpretability, as it allows us to interpret coefficients as percentage changes rather than absolute changes.

12.5 Other transformations

Square root transformation

Linear (no transformation) treats all differences in the predictor as equally important, but this may not be the case. 0 vs 10k is very different to 80 vs 90k.

Log treats doubling the same, so 5k vs 10k is the same as 80k vs 160k, which may not be appropriate.

Square root transformation squashes the larger values more gently. sqrt(10k) = 100, sqrt(40k) = 200, sqrt(90k) = 300.

However, interpretation is challenging and negative values become positive.

Idiosyncratic transformations

Problem 1: Sometimes log doesnt work, for example zero earnings. We can first model if earnings are greater than 0, and then model log earnings for those with positive earnings. This is called a hurdle approach. These two parts do not combine simply additively, they interact non-linearly. Makes coefficients hard to interpret, so best way is to simulate predictions from the model and interpret those predictions.

Problem 2: negative values. We can create discrete buckets of the predictor, e.g., -10k to 0 (-1), 0-10k (0), 10-20k (1), etc. This is called binning. We can then use these buckets as categorical predictors in the regression model.

Using continuous rather than discrete predictors

Certain variables that might appear binary or discrete may actually be better represented as continuous variables. For example, handedness is often treated as a binary variable (left-handed vs right-handed), but it can be more accurately represented as a continuous variable that captures the degree of handedness. Binary predictors can be less informative than continuous predictors, as they do not capture the full range of variation in the data.

Using discrete rather than continuous predictors

Some variables dont have a clean mathematical relationship with the outcome variable, and it may be more appropriate to use discrete predictors. For example, age may have a non-linear relationship with the outcome variable, and it may be more appropriate to use age categories (e.g., 0-18, 19-35, 36-50, 51+) rather than treating age as a continuous variable. Discrete predictors can also be easier to interpret and communicate to a non-technical audience.

Index and indicator variables

Index variables divide a population into catagories: e.g., males = 1 and females = 0, age = 1 (ages 19-29), age = 2 (ages 30-39), etc. Indicator variables are binary variables that unpack the index variables into separate binary variables. e.g., sex_1 = 1 for females and sex_2 = 1 for males.

Indicator variables, identifiability, and the baseline condition]

'Dummy variable trap' is when we have perfect multicollinearity in our model, which occurs when one predictor variable can be perfectly predicted from the others. Predictors are collinear, meaning one predictor (or a combination of them) is just a redundant copy of others.

The issue for indicator variables: the intercept is a series of 1's. If we have a indicator, one of the indicators will be 1 and the other will be 0, and the sum will always be 1. So the indicators are perfectly collinear with the intercept, which means we cannot estimate the coefficients for the indicators.

The fix is to drop one of the indicators, which becomes the baseline condition. The coefficients for the remaining indicators are then interpreted as the difference in the outcome variable between that category and the baseline category.

The clean approach is to keep intercept and drop one indicator from each categorical predictor. This allows us to interpret the coefficients as differences from the baseline category, and it avoids issues with multicollinearity.

The dropped indicator is the baseline condition, and the coefficients for the remaining indicators are interpreted as the difference in the outcome variable between that category and the baseline category.

12.6 Building and comparing regression models for prediction

Start with simple model and slowly build in complexity, checking for problems along the way.

Key choices are how input variables should be combined or transformed and which predictors to include in the model. Critical as too many predictors can cause parameter estimates to be too variable as to be useless.

General principles

1. Include input variables (for substantive reasons, not just because they are available), might be expected to be important in predicting the outcome.
2. Not always include predictors as separate, they can be combined into an index or a composite variable
3. For inputs with large effects, consider interactions
4. Use standard errors to get an idea of uncertainties
5. Decide on which predictors to include based on context and purpose of the model:
   1. If coefficient of predictor is estimated precisely (small standard error) generally keep in model
   2. If standard error is large and no substantive reason to keep predictor, consider dropping it from the model
   3. If predictor is important for the problem, then keep it in the model even if standard error is large, but be cautious in interpreting the coefficient for that predictor
   4. If coefficient doesnt make sense, try to understand why. If SE is large, could be due to random variation. If SE is small, more effort to understand why coefficient is not as expected.

Dont discard important information.

Example: predicting the yields of mesquite bushes

Aim was to develop a method of estimating total production (biomass) of mesquite leaves using easily measured parameters before harvest. Two separate sets of measurements, one on a group of 26 bushes and another on a group of 20 bushes at a different time. All from same location (ranch). Outcome variable is total weight in grams. Inputs are:

diam1: diameter of canopy in meters along longer axis diam2: canopy diameter in meters along shorter axis canopy_height: height of canopy in meters total_height: height of canopy in meters density: plant unit density (number of primary stems per plant) group: which group (0 or 1)

# Import mesquite.dat from /ros_data
mesquite = pd.read_csv('../ros_data/mesquite.dat', 
                     sep=r'\s+')

# Display the first few rows of the dataset
display(mesquite.head(10))

# print value counts of group
print(mesquite['group'].value_counts())

# change group where == 'MCD' to 0 and 'ALS' to 1
mesquite['group'] = mesquite['group'].map({'MCD': 1, 'ALS': 0})

# print value counts of group
print(mesquite['group'].value_counts())

	obs	group	diam1	diam2	total_height	canopy_height	density	weight
0	1	MCD	1.8	1.15	1.30	1.00	1	401.3
1	2	MCD	1.7	1.35	1.35	1.33	1	513.7
2	3	MCD	2.8	2.55	2.16	0.60	1	1179.2
3	4	MCD	1.3	0.85	1.80	1.20	1	308.0
4	5	MCD	3.3	1.90	1.55	1.05	1	855.2
5	6	MCD	1.4	1.40	1.20	1.00	1	268.7
6	7	MCD	1.5	0.50	1.00	0.90	1	155.5
7	8	MCD	3.9	2.30	1.70	1.30	2	1253.2
8	9	MCD	1.8	1.35	0.80	0.60	1	328.0
9	10	MCD	2.1	1.60	1.20	0.80	1	614.6

group
MCD    26
ALS    20
Name: count, dtype: int64
group
1    26
0    20
Name: count, dtype: int64

fit_and_plot_bayes(mesquite, ['diam1', 'diam2', 'canopy_height', 'total_height', 'density', 'group'], 'weight',
                       intercept_mu=0, intercept_sigma=2000,
                       slope_mu=0, slope_sigma=1000,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=False,
                       show_residuals=False,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_diam1, slope_diam2, slope_canopy_height, slope_total_height, slope_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

                         mean       sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept           -1083.358  143.596 -1362.766   -796.707      1.954   
slope_diam1           192.910   90.487    15.344    366.192      1.390   
slope_diam2           365.649  100.533   169.825    566.287      1.580   
slope_canopy_height   338.854  167.986    16.334    673.660      2.438   
slope_total_height    -90.739  149.789  -388.020    201.866      2.134   
slope_density         131.413   28.284    77.103    187.893      0.351   
slope_group           357.527   81.796   198.507    519.278      1.255   
sigma                 222.364   17.584   188.343    256.413      0.241   

                     mcse_sd  ess_bulk  ess_tail  r_hat  
intercept              1.604    5391.0    5355.0    1.0  
slope_diam1            1.001    4241.0    4963.0    1.0  
slope_diam2            1.089    4054.0    5026.0    1.0  
slope_canopy_height    1.972    4743.0    5198.0    1.0  
slope_total_height     1.622    4951.0    5350.0    1.0  
slope_density          0.336    6527.0    5161.0    1.0  
slope_group            0.869    4233.0    5075.0    1.0  
sigma                  0.198    5345.0    5277.0    1.0  

Regression formula: weight = -1083.36 + 192.91*diam1 + 365.65*diam2 + 338.85*canopy_height + -90.74*total_height + 131.41*density + 357.53*group
Residual std dev (σ): 222.36 ± 17.58
Bayesian R²: 0.836 ± 0.017
LOO-ELPD: -338.28 ± 18.37  (p_loo=22.1) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -7.354 ± 0.399
  Warning: 3 observations with Pareto-k > 0.7 (unreliable LOO estimates)

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 528kB
  Dimensions:              (chain: 4, draw: 2000)
  Coordinates:
    * chain                (chain) int64 32B 0 1 2 3
    * draw                 (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables:
      intercept            (chain, draw) float64 64kB -1.295e+03 ... -1.144e+03
      slope_diam1          (chain, draw) float64 64kB 346.0 344.2 ... 81.92 -18.51
      slope_diam2          (chain, draw) float64 64kB 269.0 242.8 ... 420.7 449.3
      slope_canopy_height  (chain, draw) float64 64kB 428.2 436.2 ... 392.0 289.2
      slope_total_height   (chain, draw) float64 64kB -124.9 -206.1 ... 193.3
      slope_density        (chain, draw) float64 64kB 113.5 108.5 ... 137.5 114.7
      slope_group          (chain, draw) float64 64kB 428.8 317.9 ... 328.2 423.0
      sigma                (chain, draw) float64 64kB 215.9 210.9 ... 209.3 225.2
  Attributes:
      created_at:                 2026-06-03T22:04:31.427352+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              4.922168016433716
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (8)
    • intercept
      (chain, draw)
      float64
      -1.295e+03 ... -1.144e+03

      array([[-1295.21561748, -1059.34899189,  -990.48946479, ...,
              -1052.32002336, -1069.90834547, -1086.33389332],
             [-1174.09307292, -1074.92036022,  -989.83488902, ...,
               -837.12352896, -1009.66244229, -1042.39443939],
             [-1061.13103554, -1032.58079046, -1420.60163058, ...,
              -1178.863225  , -1165.97258001, -1027.95601856],
             [-1043.48453645, -1265.59264153, -1260.18053439, ...,
              -1037.39677371, -1171.71572199, -1144.15534484]], shape=(4, 2000))

    • slope_diam1
      (chain, draw)
      float64
      346.0 344.2 343.4 ... 81.92 -18.51

      array([[345.96462329, 344.21192574, 343.42356847, ..., 249.25416091,
              229.68669168, 173.63316973],
             [166.19850028, 134.25154823, 136.91204768, ..., 147.0554802 ,
              112.2211384 , 310.82588146],
             [117.11153149, 213.34651724, 242.8649101 , ..., 195.19594065,
               69.44045461, 237.91679635],
             [ 99.99741855, 272.6554863 , 276.44545054, ..., 218.53186827,
               81.92496148, -18.51160002]], shape=(4, 2000))

    • slope_diam2
      (chain, draw)
      float64
      269.0 242.8 268.0 ... 420.7 449.3

      array([[268.98896104, 242.84627432, 268.03081368, ..., 333.68580943,
              380.73460339, 453.11734596],
             [337.29366879, 327.78736561, 354.74716607, ..., 460.44177097,
              468.42539882, 315.71080921],
             [429.40146628, 303.09403052, 329.16221464, ..., 348.13677225,
              635.27321357, 295.29890905],
             [486.01887277, 198.90636114, 200.00827736, ..., 375.02003981,
              420.66052597, 449.28894057]], shape=(4, 2000))

    • slope_canopy_height
      (chain, draw)
      float64
      428.2 436.2 476.2 ... 392.0 289.2

      array([[428.22336366, 436.18511149, 476.19465791, ..., 338.38487714,
              345.00795124, 307.24902901],
             [474.34810058, 400.05774342,  87.6653382 , ..., 470.01057647,
              422.7458289 , -17.48382917],
             [193.07597891, 458.08942799, 613.0631927 , ..., 315.18144066,
              396.42637236, 544.03948232],
             [523.60712817, 234.31475704, 233.84283182, ..., 230.7413324 ,
              391.96228656, 289.211767  ]], shape=(4, 2000))

    • slope_total_height
      (chain, draw)
      float64
      -124.9 -206.1 ... 34.06 193.3

      array([[-124.908505  , -206.08148913, -358.94710361, ...,  -95.41182084,
              -113.56857002, -201.89021192],
             [ -65.94669363,  -24.71384515,  149.01565108, ..., -442.41479356,
              -259.92023181,   16.4002229 ],
             [ -16.00888746, -134.606076  ,  -77.73898713, ...,  -40.0648354 ,
              -216.10607867, -235.37276932],
             [-267.75114219,  215.69328271,  194.41100245, ...,  -67.45881856,
                34.06187642,  193.2558323 ]], shape=(4, 2000))

    • slope_density
      (chain, draw)
      float64
      113.5 108.5 127.9 ... 137.5 114.7

      array([[113.49868101, 108.45468272, 127.85253688, ...,  90.98657324,
               76.01866492, 167.44593053],
             [147.17780103, 149.71660215, 137.28818617, ..., 164.08832867,
              148.10578192, 122.09913382],
             [159.53517713,  89.56590558,  97.41497821, ..., 171.67750229,
              104.35372068, 143.74345176],
             [156.22775664, 106.89959834, 103.37572488, ..., 117.9454309 ,
              137.48616404, 114.71840742]], shape=(4, 2000))

    • slope_group
      (chain, draw)
      float64
      428.8 317.9 308.8 ... 328.2 423.0

      array([[428.7737763 , 317.94623909, 308.77214732, ..., 314.43171889,
              350.83849413, 499.33393121],
             [325.95906566, 337.16444758, 276.21564   , ..., 379.36871942,
              360.70883319, 313.77271219],
             [431.29541978, 317.73695639, 388.81210654, ..., 390.71392232,
              420.89750301, 410.59566698],
             [326.36615229, 289.11740768, 286.33845676, ..., 383.41827166,
              328.21200355, 422.97812753]], shape=(4, 2000))

    • sigma
      (chain, draw)
      float64
      215.9 210.9 226.7 ... 209.3 225.2

      array([[215.89843734, 210.90530474, 226.73112823, ..., 203.01947618,
              208.09924393, 240.8018473 ],
             [189.62147591, 183.93360936, 197.92268855, ..., 280.63074993,
              242.64677826, 210.16755384],
             [254.16740715, 196.40839042, 216.53328456, ..., 210.62738477,
              243.01605012, 191.38737348],
             [206.44175485, 227.53984922, 226.80840609, ..., 203.34033295,
              209.29217727, 225.17491297]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:31.427352+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    4.922168016433716

    tuning_steps :

    1000

  
  
• log_likelihood

  <xarray.Dataset> Size: 3MB
  Dimensions:  (chain: 4, draw: 2000, y_dim_0: 46)
  Coordinates:
    * chain    (chain) int64 32B 0 1 2 3
    * draw     (draw) int64 16kB 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (chain, draw, y_dim_0) float64 3MB -6.314 -6.372 ... -6.336 -8.423
  Attributes:
      created_at:                 2026-06-03T22:04:31.700965+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • y_dim_0: 46
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (chain, draw, y_dim_0)
      float64
      -6.314 -6.372 ... -6.336 -8.423

      array([[[ -6.31420602,  -6.37224014,  -7.19829227, ...,  -6.29795741,
                -6.36413791,  -9.37862755],
              [ -6.28250436,  -6.32273987,  -8.18323093, ...,  -6.27179333,
                -6.35559889,  -8.10101441],
              [ -6.3462902 ,  -6.34971906,  -9.58771201, ...,  -6.37422977,
                -6.78965208,  -7.95345033],
              ...,
              [ -6.23226617,  -6.24697554,  -7.18568475, ...,  -6.23316818,
                -6.25120224,  -8.19826272],
              [ -6.25716697,  -6.28969396,  -6.95884377, ...,  -6.25987649,
                -6.27475842,  -8.34765344],
              [ -6.43148407,  -6.50671917,  -6.77908212, ...,  -6.40319376,
                -6.71816382,  -8.18022692]],
      
             [[ -6.17380296,  -6.22286973,  -8.53055239, ...,  -6.25765082,
                -6.21383889,  -9.08885971],
              [ -6.13362938,  -6.20127905,  -8.14393542, ...,  -6.58047626,
                -6.25891618,  -8.6802889 ],
              [ -6.22916385,  -6.25369727,  -6.31744411, ...,  -6.39508369,
                -6.21907516,  -7.73461532],
      ...
              [ -6.26912822,  -6.28425035,  -7.02392598, ...,  -6.31627038,
                -6.41464167,  -8.71749154],
              [ -6.45471627,  -6.4393328 ,  -6.8775364 , ...,  -6.62558464,
                -6.45098066,  -8.56856196],
              [ -6.407427  ,  -6.84364668,  -8.53118769, ...,  -6.19808079,
                -6.2827252 ,  -8.13914044]],
      
             [[ -6.27569879,  -6.33666593,  -9.00781454, ...,  -6.25997141,
                -6.37550629,  -8.21569518],
              [ -6.3592269 ,  -6.36889899,  -6.63077932, ...,  -6.53636714,
                -6.35012899,  -8.60201358],
              [ -6.37526696,  -6.39167198,  -6.75296741, ...,  -6.60922032,
                -6.34450711,  -8.72817028],
              ...,
              [ -6.24495275,  -6.25708734,  -6.52176187, ...,  -6.23631143,
                -6.31181267,  -8.16283374],
              [ -6.28079683,  -6.29687348,  -7.1220945 , ...,  -6.29279628,
                -6.26270243,  -8.5605088 ],
              [ -6.33835533,  -6.43533455,  -6.43610133, ...,  -6.76281557,
                -6.33630846,  -8.42291038]]], shape=(4, 2000, 46))

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:31.700965+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  
• sample_stats

  <xarray.Dataset> Size: 1MB
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 32B 0 1 2 3
    * draw                   (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables: (12/18)
      energy                 (chain, draw) float64 64kB 389.8 388.6 ... 391.8
      divergences            (chain, draw) int64 64kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
      acceptance_rate        (chain, draw) float64 64kB 0.9805 0.9863 ... 0.9892
      step_size              (chain, draw) float64 64kB 0.1375 0.1375 ... 0.1186
      lp                     (chain, draw) float64 64kB -386.9 -386.0 ... -389.4
      step_size_bar          (chain, draw) float64 64kB 0.1124 0.1124 ... 0.09822
      ...                     ...
      tree_depth             (chain, draw) int64 64kB 6 6 5 4 5 6 ... 6 6 6 5 5 5
      energy_error           (chain, draw) float64 64kB -0.4797 ... -0.01394
      perf_counter_start     (chain, draw) float64 64kB 399.0 399.0 ... 402.0
      perf_counter_diff      (chain, draw) float64 64kB 0.001977 ... 0.0009177
      index_in_trajectory    (chain, draw) int64 64kB 32 11 -10 3 ... 15 15 -10
      reached_max_treedepth  (chain, draw) bool 8kB False False ... False False
  Attributes:
      created_at:                 2026-06-03T22:04:31.446747+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              4.922168016433716
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (18)
    • energy
      (chain, draw)
      float64
      389.8 388.6 388.0 ... 387.7 391.8

      array([[389.79329835, 388.63006698, 388.02736348, ..., 386.64254528,
              389.42416876, 392.7421252 ],
             [394.09103211, 392.37224628, 393.295322  , ..., 395.67122814,
              397.52650926, 393.78943968],
             [392.47679206, 390.97002639, 391.93966656, ..., 386.94083025,
              393.57211914, 401.04143615],
             [387.85037074, 391.26734183, 389.71777657, ..., 386.64912427,
              387.6917664 , 391.80218006]], shape=(4, 2000))

    • divergences
      (chain, draw)
      int64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0]], shape=(4, 2000))

    • acceptance_rate
      (chain, draw)
      float64
      0.9805 0.9863 ... 0.9998 0.9892

      array([[0.98045885, 0.98625688, 0.99684193, ..., 0.98440536, 0.95042112,
              0.60858401],
             [0.67791869, 0.98840242, 0.2945264 , ..., 0.99897388, 0.57648689,
              0.82036503],
             [0.89691948, 0.7923745 , 0.99773543, ..., 0.99990383, 0.94597983,
              0.5912883 ],
             [0.99624436, 0.99029901, 0.90592408, ..., 0.98965885, 0.9998021 ,
              0.98924038]], shape=(4, 2000))

    • step_size
      (chain, draw)
      float64
      0.1375 0.1375 ... 0.1186 0.1186

      array([[0.13746654, 0.13746654, 0.13746654, ..., 0.13746654, 0.13746654,
              0.13746654],
             [0.09750623, 0.09750623, 0.09750623, ..., 0.09750623, 0.09750623,
              0.09750623],
             [0.09603594, 0.09603594, 0.09603594, ..., 0.09603594, 0.09603594,
              0.09603594],
             [0.11855127, 0.11855127, 0.11855127, ..., 0.11855127, 0.11855127,
              0.11855127]], shape=(4, 2000))

    • lp
      (chain, draw)
      float64
      -386.9 -386.0 ... -386.1 -389.4

      array([[-386.85532949, -385.97021678, -387.09329624, ..., -385.18768286,
              -386.51185693, -387.66302022],
             [-387.21577441, -387.38738603, -386.72556987, ..., -392.17748809,
              -387.65910433, -389.54042488],
             [-387.33281619, -387.51996832, -388.06244638, ..., -385.88932388,
              -390.02117594, -393.39786035],
             [-386.74863953, -387.48890882, -387.75148423, ..., -385.04020221,
              -386.12634517, -389.4349833 ]], shape=(4, 2000))

    • step_size_bar
      (chain, draw)
      float64
      0.1124 0.1124 ... 0.09822 0.09822

      array([[0.11239506, 0.11239506, 0.11239506, ..., 0.11239506, 0.11239506,
              0.11239506],
             [0.10375925, 0.10375925, 0.10375925, ..., 0.10375925, 0.10375925,
              0.10375925],
             [0.10754418, 0.10754418, 0.10754418, ..., 0.10754418, 0.10754418,
              0.10754418],
             [0.0982165 , 0.0982165 , 0.0982165 , ..., 0.0982165 , 0.0982165 ,
              0.0982165 ]], shape=(4, 2000))

    • max_energy_error
      (chain, draw)
      float64
      -0.744 -0.194 ... -0.06909 -0.04329

      array([[-0.74395054, -0.19395078, -0.19257878, ...,  0.03919233,
               0.1023918 ,  1.03379261],
             [ 1.62057741,  0.10529208,  2.85421484, ..., -0.2330241 ,
               1.91610737, -1.14386763],
             [ 0.1947331 ,  0.96514871, -0.55233675, ..., -0.4030706 ,
               0.19442983,  2.76590645],
             [-0.02244086,  0.04026198,  0.20195285, ..., -0.10669825,
              -0.06908624, -0.04329308]], shape=(4, 2000))

    • n_steps
      (chain, draw)
      float64
      63.0 63.0 31.0 ... 31.0 31.0 31.0

      array([[63., 63., 31., ..., 31., 31., 31.],
             [63., 31., 63., ..., 31., 19., 63.],
             [31., 95., 63., ..., 63., 63., 63.],
             [63., 63., 31., ..., 31., 31., 31.]], shape=(4, 2000))

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

    • process_time_diff
      (chain, draw)
      float64
      0.001977 0.001979 ... 0.000918

      array([[0.001977, 0.001979, 0.001004, ..., 0.000967, 0.000963, 0.000937],
             [0.001923, 0.000968, 0.00192 , ..., 0.001064, 0.00062 , 0.001949],
             [0.000933, 0.002851, 0.001973, ..., 0.002022, 0.001977, 0.003307],
             [0.001864, 0.00184 , 0.000917, ..., 0.000919, 0.000924, 0.000918]],
            shape=(4, 2000))

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • tree_depth
      (chain, draw)
      int64
      6 6 5 4 5 6 6 5 ... 5 4 6 6 6 5 5 5

      array([[6, 6, 5, ..., 5, 5, 5],
             [6, 5, 6, ..., 5, 5, 6],
             [5, 7, 6, ..., 6, 6, 6],
             [6, 6, 5, ..., 5, 5, 5]], shape=(4, 2000))

    • energy_error
      (chain, draw)
      float64
      -0.4797 -0.02139 ... -0.01394

      array([[-4.79714014e-01, -2.13894024e-02, -1.23435973e-02, ...,
              -6.23185833e-04,  8.07863835e-02,  3.77990416e-01],
             [-6.91108296e-01, -6.37445984e-02,  3.03068665e-01, ...,
              -9.63285521e-02,  4.28612532e-02,  2.69701209e-01],
             [ 6.06859687e-02,  2.45267479e-01, -4.09964407e-01, ...,
              -1.86951543e-01,  1.75295713e-01,  2.07666746e+00],
             [ 8.37896221e-03, -9.82365408e-03,  2.01952847e-01, ...,
              -3.55615096e-02, -7.75491598e-03, -1.39398463e-02]],
            shape=(4, 2000))

    • perf_counter_start
      (chain, draw)
      float64
      399.0 399.0 399.0 ... 402.0 402.0

      array([[399.02834912, 399.03038904, 399.03245387, ..., 401.35806438,
              401.35907479, 401.36007913],
             [399.00591425, 399.00788975, 399.00890246, ..., 401.73499879,
              401.73615229, 401.73682579],
             [398.99750875, 398.99849958, 399.00142629, ..., 401.641765  ,
              401.64389417, 401.645925  ],
             [399.11485   , 399.11677221, 399.11865542, ..., 402.01963946,
              402.02059925, 402.02156375]], shape=(4, 2000))

    • perf_counter_diff
      (chain, draw)
      float64
      0.001977 0.002012 ... 0.0009177

      array([[0.00197738, 0.00201221, 0.001005  , ..., 0.00096742, 0.00096242,
              0.0009375 ],
             [0.00192275, 0.00096888, 0.00191908, ..., 0.00107025, 0.00062108,
              0.00194863],
             [0.00093308, 0.00286146, 0.00197267, ..., 0.00205592, 0.00197675,
              0.01387246],
             [0.00186408, 0.00184121, 0.0009165 , ..., 0.00091925, 0.00092292,
              0.00091771]], shape=(4, 2000))

    • index_in_trajectory
      (chain, draw)
      int64
      32 11 -10 3 15 ... -10 15 15 -10

      array([[ 32,  11, -10, ...,  12,  -3, -26],
             [-16,  -3,   8, ...,   8,  -8, -34],
             [-26,  32, -14, ..., -14,  33,  17],
             [ 20, -46,  -1, ...,  15,  15, -10]], shape=(4, 2000))

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:31.446747+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    4.922168016433716

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset> Size: 736B
  Dimensions:  (y_dim_0: 46)
  Coordinates:
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (y_dim_0) float64 368B 401.3 513.7 1.179e+03 ... 908.0 213.5
  Attributes:
      created_at:                 2026-06-03T22:04:31.450405+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • y_dim_0: 46
  • Coordinates: (1)
    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (y_dim_0)
      float64
      401.3 513.7 ... 908.0 213.5

      array([ 401.3,  513.7, 1179.2,  308. ,  855.2,  268.7,  155.5, 1253.2,
              328. ,  614.6,   60.2,  269.6,  448.4,  120.4,  378.7,  266.4,
              138.9, 1020.8,  635.7,  621.8,  579.8,  326.8,   66.7,   68. ,
              153.1,  256.4,  723. , 4052. ,  345. ,  330.9,  163.5, 1160. ,
              386.6,  693.5,  674.4,  217.5,  771.3,  341.7,  125.7,  462.5,
               64.5,  850.6,  226. , 1745.1,  908. ,  213.5])

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:31.450405+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  

LOO computation unstable and p_loo is greater than number of predictors, which suggests normal distribution is not good for this data.

Now try log scale

mesquite['log_weight'] = np.log(mesquite['weight'])
mesquite['log_diam1'] = np.log(mesquite['diam1'])
mesquite['log_diam2'] = np.log(mesquite['diam2'])
mesquite['log_canopy_height'] = np.log(mesquite['canopy_height'])
mesquite['log_total_height'] = np.log(mesquite['total_height'])
mesquite['log_density'] = np.log(mesquite['density'])

fit_and_plot_bayes(mesquite, ['log_diam1', 'log_diam2', 'log_canopy_height', 'log_total_height', 'log_density', 'group'], 'log_weight',
                       intercept_mu=0, intercept_sigma=2000,
                       slope_mu=0, slope_sigma=1000,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=False,
                       show_residuals=False,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_diam1, slope_log_diam2, slope_log_canopy_height, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.769  0.164     4.440      5.075      0.003   
slope_log_diam1          0.392  0.300    -0.195      0.980      0.005   
slope_log_diam2          1.154  0.219     0.713      1.570      0.004   
slope_log_canopy_height  0.373  0.291    -0.207      0.932      0.004   
slope_log_total_height   0.393  0.326    -0.267      0.996      0.005   
slope_log_density        0.108  0.128    -0.128      0.377      0.002   
slope_group              0.584  0.135     0.314      0.840      0.002   
sigma                    0.341  0.040     0.269      0.422      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    4151.0    4767.0    1.0  
slope_log_diam1            0.004    3582.0    4006.0    1.0  
slope_log_diam2            0.003    3763.0    4406.0    1.0  
slope_log_canopy_height    0.003    4415.0    5252.0    1.0  
slope_log_total_height     0.004    4443.0    5188.0    1.0  
slope_log_density          0.001    5025.0    5045.0    1.0  
slope_group                0.002    3905.0    4926.0    1.0  
sigma                      0.000    5287.0    3915.0    1.0  

Regression formula: log_weight = 4.77 + 0.39*log_diam1 + 1.15*log_diam2 + 0.37*log_canopy_height + 0.39*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.016
LOO-ELPD: -19.51 ± 5.28  (p_loo=7.7) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.424 ± 0.115
  Warning: 1 observations with Pareto-k > 0.7 (unreliable LOO estimates)

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 528kB
  Dimensions:                  (chain: 4, draw: 2000)
  Coordinates:
    * chain                    (chain) int64 32B 0 1 2 3
    * draw                     (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables:
      intercept                (chain, draw) float64 64kB 4.523 4.791 ... 4.823
      slope_log_diam1          (chain, draw) float64 64kB 0.3947 0.08292 ... 0.453
      slope_log_diam2          (chain, draw) float64 64kB 1.224 1.423 ... 1.043
      slope_log_canopy_height  (chain, draw) float64 64kB 0.702 0.2439 ... 0.5961
      slope_log_total_height   (chain, draw) float64 64kB 0.1876 0.3326 ... 0.1526
      slope_log_density        (chain, draw) float64 64kB 0.2322 0.1079 ... 0.1981
      slope_group              (chain, draw) float64 64kB 0.7947 0.7048 ... 0.6119
      sigma                    (chain, draw) float64 64kB 0.4057 0.3209 ... 0.3171
  Attributes:
      created_at:                 2026-06-03T22:04:35.866843+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              2.322632074356079
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (8)
    • intercept
      (chain, draw)
      float64
      4.523 4.791 4.621 ... 4.624 4.823

      array([[4.52306658, 4.79139983, 4.620577  , ..., 4.74944959, 4.79304058,
              4.94728519],
             [4.62725123, 4.63393667, 4.90065842, ..., 4.8686287 , 4.92899512,
              4.95479439],
             [4.99599201, 4.90005418, 4.84817342, ..., 4.9772962 , 4.85890828,
              4.83150664],
             [4.68701873, 4.85488944, 4.67505979, ..., 4.88645925, 4.62367648,
              4.8227456 ]], shape=(4, 2000))

    • slope_log_diam1
      (chain, draw)
      float64
      0.3947 0.08292 ... 0.4691 0.453

      array([[0.39469215, 0.08292362, 0.21499611, ..., 0.46316476, 0.56800653,
              0.41135803],
             [0.64263263, 0.18455961, 0.14028155, ..., 0.27063137, 0.40430396,
              0.09896914],
             [0.3615852 , 0.22785768, 0.39205514, ..., 0.11769716, 0.12476762,
              0.02123823],
             [0.45317446, 0.5887562 , 0.78941541, ..., 0.42025993, 0.46913459,
              0.45303876]], shape=(4, 2000))

    • slope_log_diam2
      (chain, draw)
      float64
      1.224 1.423 1.457 ... 1.296 1.043

      array([[1.22405091, 1.42344587, 1.45681969, ..., 1.03099537, 1.04344657,
              1.02214654],
             [1.15015904, 1.06651739, 1.23814702, ..., 1.1462218 , 1.02556702,
              1.39157791],
             [1.07505814, 1.2520462 , 1.37283842, ..., 1.12149654, 1.50392809,
              1.5280253 ],
             [1.19740875, 1.09541271, 0.97772579, ..., 0.98511589, 1.29579059,
              1.04343322]], shape=(4, 2000))

    • slope_log_canopy_height
      (chain, draw)
      float64
      0.702 0.2439 ... 0.5808 0.5961

      array([[ 0.70203413,  0.24387596,  0.39456364, ...,  0.32339621,
               0.28320237,  0.68441252],
             [ 0.54027734,  0.20972752, -0.11524851, ...,  0.64839818,
               0.57859937,  0.46059321],
             [ 0.8191074 ,  0.91001027,  1.00894192, ...,  0.70861079,
               0.00691244,  0.32013303],
             [ 0.74561921,  0.35069143,  0.40023505, ...,  0.50551953,
               0.58079947,  0.59606286]], shape=(4, 2000))

    • slope_log_total_height
      (chain, draw)
      float64
      0.1876 0.3326 ... 0.09447 0.1526

      array([[ 0.18763083,  0.33257278,  0.27562338, ...,  0.45261622,
               0.21106792,  0.3640739 ],
             [ 0.01480973,  0.95026184,  0.5586688 , ...,  0.13873504,
               0.11998894,  0.50802801],
             [-0.04274233, -0.09986234, -0.41512673, ...,  0.48969925,
               0.43447097,  0.60168422],
             [ 0.3948423 ,  0.11493473,  0.52677674, ...,  0.22271891,
               0.0944665 ,  0.15255079]], shape=(4, 2000))

    • slope_log_density
      (chain, draw)
      float64
      0.2322 0.1079 ... 0.08496 0.1981

      array([[ 0.23221212,  0.10788576,  0.13526167, ...,  0.15006375,
               0.13451956, -0.03849132],
             [ 0.1818731 ,  0.11402135,  0.22387892, ...,  0.15970507,
               0.24228288, -0.0908581 ],
             [ 0.03389838,  0.16529202,  0.05622013, ..., -0.02579837,
               0.10762364,  0.14201655],
             [-0.06860991, -0.04992482, -0.0189041 , ...,  0.09245911,
               0.08495961,  0.19808922]], shape=(4, 2000))

    • slope_group
      (chain, draw)
      float64
      0.7947 0.7048 ... 0.7717 0.6119

      array([[0.79470418, 0.70478619, 0.75531851, ..., 0.61680667, 0.42111718,
              0.60890985],
             [0.64769046, 0.72083993, 0.42749697, ..., 0.58336292, 0.59102818,
              0.4555297 ],
             [0.62853037, 0.57351789, 0.69872832, ..., 0.51717415, 0.64921587,
              0.60707226],
             [0.58006198, 0.54021222, 0.48896335, ..., 0.44424215, 0.77165785,
              0.61191529]], shape=(4, 2000))

    • sigma
      (chain, draw)
      float64
      0.4057 0.3209 ... 0.2835 0.3171

      array([[0.40565127, 0.32094027, 0.33070171, ..., 0.37964758, 0.38250827,
              0.32933283],
             [0.33976062, 0.38570749, 0.30710201, ..., 0.30846503, 0.29503802,
              0.3189579 ],
             [0.27736737, 0.26778256, 0.34381213, ..., 0.35580903, 0.35547183,
              0.36488723],
             [0.35850137, 0.3287538 , 0.31035539, ..., 0.33013534, 0.28350771,
              0.31709812]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:35.866843+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    2.322632074356079

    tuning_steps :

    1000

  
  
• log_likelihood

  <xarray.Dataset> Size: 3MB
  Dimensions:  (chain: 4, draw: 2000, y_dim_0: 46)
  Coordinates:
    * chain    (chain) int64 32B 0 1 2 3
    * draw     (draw) int64 16kB 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (chain, draw, y_dim_0) float64 3MB -0.17 -0.04161 ... -2.422
  Attributes:
      created_at:                 2026-06-03T22:04:36.119808+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • y_dim_0: 46
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (chain, draw, y_dim_0)
      float64
      -0.17 -0.04161 ... 0.02269 -2.422

      array([[[-1.70011627e-01, -4.16054764e-02, -5.44368858e-01, ...,
               -1.09543581e-01, -3.42816580e-01, -3.67644255e+00],
              [ 8.76649368e-02,  1.64131050e-01,  2.14020719e-01, ...,
               -8.99090969e-02, -1.33883959e-01, -2.13883842e+00],
              [-2.67355669e-02,  1.22611556e-01,  1.41048740e-01, ...,
                8.33670283e-02, -1.14684159e-01, -3.54259823e+00],
              ...,
              [ 1.93325837e-02,  2.01153407e-02,  2.68245278e-02, ...,
               -1.49171781e-02, -3.08131028e-02, -1.63000360e+00],
              [-1.63853999e-01, -2.04907806e-01, -2.23869741e-01, ...,
               -9.38494281e-03, -8.92490652e-02, -1.13870793e+00],
              [ 1.83749400e-01,  9.61219955e-02, -2.49539413e-03, ...,
                1.69728167e-01,  1.74746416e-01, -1.91217077e+00]],
      
             [[ 2.42602177e-02,  9.60964058e-02, -2.93857602e-01, ...,
                1.53744586e-01, -8.41633336e-02, -2.91233290e+00],
              [-2.57827536e-02, -1.78359917e-02,  3.25114089e-03, ...,
               -3.18368235e-01, -7.00896656e-02, -2.48142076e+00],
              [-1.09200017e-01, -3.25121583e-01,  2.49356207e-01, ...,
                1.11754356e-02, -3.83865212e-02, -6.77744868e-01],
      ...
               -6.83684216e-02,  1.13879734e-01, -2.03534987e+00],
              [ 8.39807238e-02,  8.77050306e-02, -2.47604028e-01, ...,
                1.07947223e-01,  7.54261274e-02, -9.31820866e-01],
              [-2.21633842e-02,  7.50891894e-02,  3.68864899e-02, ...,
                6.39041770e-02,  8.11257700e-02, -1.91342472e+00]],
      
             [[-3.40218297e-02,  9.94838329e-02, -2.31299863e-01, ...,
                1.06402018e-01,  9.58798392e-02, -3.46098013e+00],
              [ 1.70648966e-01,  1.70261223e-01,  1.07289476e-01, ...,
                1.16141463e-01,  1.33775269e-01, -1.18186927e+00],
              [ 2.07331833e-01,  2.06121072e-01,  2.48875145e-01, ...,
                1.43571119e-01,  1.34119934e-01, -2.55430762e+00],
              ...,
              [-3.44846070e-02,  3.85276436e-02, -8.40265183e-01, ...,
               -8.95761029e-02, -1.07893593e-02, -1.67193325e+00],
              [ 2.55345819e-01,  3.40266879e-01,  7.96824898e-02, ...,
                2.92077658e-01,  7.44071248e-02, -4.84778202e+00],
              [ 1.71710294e-01,  2.22546057e-01, -4.94973348e-01, ...,
                1.52061030e-01,  2.26852476e-02, -2.42222389e+00]]],
            shape=(4, 2000, 46))

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:36.119808+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  
• sample_stats

  <xarray.Dataset> Size: 1MB
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 32B 0 1 2 3
    * draw                   (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables: (12/18)
      energy                 (chain, draw) float64 64kB 83.64 80.83 ... 80.59
      divergences            (chain, draw) int64 64kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
      acceptance_rate        (chain, draw) float64 64kB 0.8935 0.9497 ... 0.5642
      step_size              (chain, draw) float64 64kB 0.2003 0.2003 ... 0.246
      lp                     (chain, draw) float64 64kB -77.77 -74.65 ... -72.6
      step_size_bar          (chain, draw) float64 64kB 0.2122 0.2122 ... 0.1969
      ...                     ...
      tree_depth             (chain, draw) int64 64kB 3 4 4 5 4 4 ... 5 5 4 4 4 4
      energy_error           (chain, draw) float64 64kB 0.1166 0.05265 ... 0.1894
      perf_counter_start     (chain, draw) float64 64kB 405.0 405.0 ... 406.5
      perf_counter_diff      (chain, draw) float64 64kB 0.0002503 ... 0.0004614
      index_in_trajectory    (chain, draw) int64 64kB 7 9 -4 -2 ... 11 -15 8 -8
      reached_max_treedepth  (chain, draw) bool 8kB False False ... False False
  Attributes:
      created_at:                 2026-06-03T22:04:35.880251+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              2.322632074356079
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (18)
    • energy
      (chain, draw)
      float64
      83.64 80.83 76.51 ... 76.95 80.59

      array([[83.63799958, 80.82827164, 76.50792295, ..., 78.01033089,
              76.35643135, 78.98765688],
             [77.0181942 , 77.89032031, 80.82179473, ..., 76.29400802,
              77.21953871, 79.65331092],
             [80.78773379, 82.63457298, 80.64374255, ..., 81.0187518 ,
              80.58818876, 79.73567475],
             [78.77309651, 77.03395012, 81.00661414, ..., 77.56248811,
              76.94525988, 80.59057889]], shape=(4, 2000))

    • divergences
      (chain, draw)
      int64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0]], shape=(4, 2000))

    • acceptance_rate
      (chain, draw)
      float64
      0.8935 0.9497 1.0 ... 0.9147 0.5642

      array([[0.89354298, 0.94972594, 1.        , ..., 0.99348139, 0.87057377,
              0.76870206],
             [0.90589783, 0.79875419, 0.85240993, ..., 0.78564464, 0.92854361,
              1.        ],
             [0.99211326, 0.2795574 , 1.        , ..., 0.95145373, 0.87302735,
              0.60914716],
             [0.95655269, 0.97588167, 0.35708116, ..., 0.6879551 , 0.91468283,
              0.56418084]], shape=(4, 2000))

    • step_size
      (chain, draw)
      float64
      0.2003 0.2003 ... 0.246 0.246

      array([[0.20025894, 0.20025894, 0.20025894, ..., 0.20025894, 0.20025894,
              0.20025894],
             [0.17559791, 0.17559791, 0.17559791, ..., 0.17559791, 0.17559791,
              0.17559791],
             [0.22333205, 0.22333205, 0.22333205, ..., 0.22333205, 0.22333205,
              0.22333205],
             [0.24595981, 0.24595981, 0.24595981, ..., 0.24595981, 0.24595981,
              0.24595981]], shape=(4, 2000))

    • lp
      (chain, draw)
      float64
      -77.77 -74.65 ... -74.53 -72.6

      array([[-77.77364963, -74.65071863, -75.34766158, ..., -73.52684462,
              -74.72060418, -76.91833306],
             [-73.53527327, -76.05599206, -76.00358717, ..., -72.80401556,
              -75.58558418, -75.09402169],
             [-77.5526898 , -77.82635925, -75.73811198, ..., -75.12053796,
              -75.79132733, -75.55099578],
             [-75.96095163, -74.889971  , -75.2044652 , ..., -73.94921293,
              -74.52652618, -72.59924992]], shape=(4, 2000))

    • step_size_bar
      (chain, draw)
      float64
      0.2122 0.2122 ... 0.1969 0.1969

      array([[0.21217098, 0.21217098, 0.21217098, ..., 0.21217098, 0.21217098,
              0.21217098],
             [0.20893847, 0.20893847, 0.20893847, ..., 0.20893847, 0.20893847,
              0.20893847],
             [0.20379356, 0.20379356, 0.20379356, ..., 0.20379356, 0.20379356,
              0.20379356],
             [0.19689628, 0.19689628, 0.19689628, ..., 0.19689628, 0.19689628,
              0.19689628]], shape=(4, 2000))

    • max_energy_error
      (chain, draw)
      float64
      0.2425 -0.3071 ... -0.8014 1.368

      array([[ 0.24245675, -0.30713343, -0.73607178, ..., -1.79312014,
               0.79038763,  1.0970607 ],
             [ 0.27311463,  0.70907091,  0.88565862, ...,  0.68936884,
              -0.45322284, -0.66603159],
             [-0.19828054,  3.29739995, -1.03783727, ..., -0.32589074,
               0.39561495,  1.40335607],
             [ 0.18803216, -0.12922279,  2.69106713, ...,  0.92248044,
              -0.80140382,  1.3675017 ]], shape=(4, 2000))

    • n_steps
      (chain, draw)
      float64
      7.0 15.0 15.0 ... 15.0 15.0 15.0

      array([[ 7., 15., 15., ..., 15., 31., 15.],
             [15., 15., 15., ..., 15.,  7., 15.],
             [15., 31., 15., ..., 15., 31.,  7.],
             [31.,  7., 15., ..., 15., 15., 15.]], shape=(4, 2000))

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

    • process_time_diff
      (chain, draw)
      float64
      0.000251 0.000489 ... 0.000461

      array([[0.000251, 0.000489, 0.000488, ..., 0.000484, 0.000959, 0.000485],
             [0.000489, 0.000486, 0.000486, ..., 0.000458, 0.00023 , 0.000461],
             [0.000486, 0.000929, 0.000472, ..., 0.000485, 0.000967, 0.000323],
             [0.001031, 0.00025 , 0.000487, ..., 0.00046 , 0.00046 , 0.000461]],
            shape=(4, 2000))

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • tree_depth
      (chain, draw)
      int64
      3 4 4 5 4 4 4 5 ... 5 5 5 5 4 4 4 4

      array([[3, 4, 4, ..., 4, 5, 4],
             [4, 4, 4, ..., 4, 3, 4],
             [4, 5, 4, ..., 4, 5, 3],
             [5, 3, 4, ..., 4, 4, 4]], shape=(4, 2000))

    • energy_error
      (chain, draw)
      float64
      0.1166 0.05265 ... -0.7715 0.1894

      array([[ 0.11656137,  0.05265274, -0.03670748, ..., -1.2423376 ,
               0.13631339,  0.50626299],
             [-0.06966241,  0.26884964, -0.14687077, ...,  0.0916119 ,
               0.27440223, -0.56854159],
             [ 0.05405064,  0.49460807, -0.5965159 , ..., -0.20141802,
               0.07999129,  0.28083047],
             [ 0.07503271, -0.03103203,  0.84909235, ...,  0.51601818,
              -0.77153113,  0.18943169]], shape=(4, 2000))

    • perf_counter_start
      (chain, draw)
      float64
      405.0 405.0 405.0 ... 406.5 406.5

      array([[404.97670017, 404.977001  , 404.977535  , ..., 406.33187213,
              406.33239792, 406.33339942],
             [404.98507438, 404.98561396, 404.98614404, ..., 406.37970183,
              406.38019829, 406.38046658],
             [404.97769567, 404.97822829, 404.97919988, ..., 406.462579  ,
              406.46310642, 406.46414054],
             [404.99879779, 404.99990488, 405.00019925, ..., 406.48817371,
              406.48867417, 406.48917467]], shape=(4, 2000))

    • perf_counter_diff
      (chain, draw)
      float64
      0.0002503 0.0004898 ... 0.0004614

      array([[0.00025029, 0.00048979, 0.00048883, ..., 0.00048375, 0.00095908,
              0.00048562],
             [0.00048879, 0.00048554, 0.00048492, ..., 0.00045779, 0.00022979,
              0.00046217],
             [0.00048692, 0.00092896, 0.00047254, ..., 0.00048525, 0.00096758,
              0.00036371],
             [0.00104254, 0.00025042, 0.00048775, ..., 0.00046063, 0.00046058,
              0.00046142]], shape=(4, 2000))

    • index_in_trajectory
      (chain, draw)
      int64
      7 9 -4 -2 -11 8 ... 12 11 -15 8 -8

      array([[  7,   9,  -4, ..., -11,  -7,   7],
             [ -8,  -8, -11, ...,  -9,   2,  -7],
             [-11,  -3,   7, ...,  13, -30,   3],
             [-10,  -7,  -6, ..., -15,   8,  -8]], shape=(4, 2000))

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:35.880251+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    2.322632074356079

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset> Size: 736B
  Dimensions:  (y_dim_0: 46)
  Coordinates:
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (y_dim_0) float64 368B 5.995 6.242 7.073 5.73 ... 7.465 6.811 5.364
  Attributes:
      created_at:                 2026-06-03T22:04:35.883937+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • y_dim_0: 46
  • Coordinates: (1)
    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (y_dim_0)
      float64
      5.995 6.242 7.073 ... 6.811 5.364

      array([5.99470928, 6.24163944, 7.07259152, 5.73009978, 6.75133536,
             5.59359552, 5.04664573, 7.13345556, 5.79301361, 6.42097165,
             4.09767235, 5.59693938, 6.10568569, 4.79081953, 5.93674433,
             5.58499894, 4.93375425, 6.92834191, 6.45472675, 6.4326185 ,
             6.36268322, 5.78934836, 4.20020495, 4.21950771, 5.0310913 ,
             5.54673873, 6.58340922, 8.30696587, 5.84354442, 5.80181621,
             5.09681299, 7.05617528, 5.95739057, 6.54175124, 6.51382341,
             5.38219885, 6.6480774 , 5.83393316, 4.83389812, 6.13664656,
             4.16666522, 6.74594198, 5.420535  , 7.46456714, 6.81124438,
             5.36363683])

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:35.883937+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  

ELPD for two models: -341 and -19 are not comparable as one was transformed. When comparing models with transformed continuous outcomes, we must account for the transformation.

# Re-fit both models (no plots) to get traces for comparison
trace_raw = fit_and_plot_bayes(
    mesquite, ['diam1', 'diam2', 'canopy_height', 'total_height', 'density', 'group'], 'weight',
    intercept_mu=0, intercept_sigma=2000, slope_mu=0, slope_sigma=1000, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

trace_log = fit_and_plot_bayes(
    mesquite, ['log_diam1', 'log_diam2', 'log_canopy_height', 'log_total_height', 'log_density', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=2000, slope_mu=0, slope_sigma=1000, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

# Pointwise LOO for both models
loo_raw = az.loo(trace_raw, pointwise=True)
loo_log = az.loo(trace_log, pointwise=True)

# Jacobian adjustment: outcome is log(weight), so log|Jacobian| = log(1/weight) = -log(weight)
# Subtract log(weight) from each pointwise ELPD to put the log model on the original scale
loo_log_jacobian = loo_log.copy()
loo_log_jacobian.loo_i.values -= np.log(mesquite['weight'].values)

n = len(loo_log_jacobian.loo_i)
loo_log_jacobian["elpd_loo"] = float(loo_log_jacobian.loo_i.sum())
loo_log_jacobian["se"] = float((n * loo_log_jacobian.loo_i.var()) ** 0.5)

print(f"LOO raw model:              {loo_raw.elpd_loo:.1f} ± {loo_raw.se:.1f}")
print(f"LOO log model (unadjusted): {loo_log.elpd_loo:.1f} ± {loo_log.se:.1f}")
print(f"LOO log model (adjusted):   {loo_log_jacobian['elpd_loo']:.1f} ± {loo_log_jacobian['se']:.1f}")
print()
az.compare({"raw": loo_raw, "log_jacobian": loo_log_jacobian})

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_diam1, slope_diam2, slope_canopy_height, slope_total_height, slope_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

                         mean       sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept           -1083.922  145.081 -1367.901   -803.578      1.924   
slope_diam1           191.674   93.663    15.446    384.091      1.449   
slope_diam2           367.219  101.872   166.836    563.528      1.505   
slope_canopy_height   339.323  169.671    -6.565    659.120      2.500   
slope_total_height    -90.681  150.924  -377.898    210.347      2.071   
slope_density         131.169   28.898    72.864    187.342      0.351   
slope_group           357.164   81.415   184.371    501.971      1.141   
sigma                 222.602   18.064   188.044    258.261      0.228   

                     mcse_sd  ess_bulk  ess_tail  r_hat  
intercept              1.579    5686.0    5086.0    1.0  
slope_diam1            1.119    4175.0    4620.0    1.0  
slope_diam2            1.164    4555.0    5456.0    1.0  
slope_canopy_height    1.821    4607.0    4799.0    1.0  
slope_total_height     1.643    5318.0    5180.0    1.0  
slope_density          0.344    6770.0    5646.0    1.0  
slope_group            0.835    5093.0    5458.0    1.0  
sigma                  0.215    6198.0    5051.0    1.0  

Regression formula: weight = -1083.92 + 191.67*diam1 + 367.22*diam2 + 339.32*canopy_height + -90.68*total_height + 131.17*density + 357.16*group
Residual std dev (σ): 222.60 ± 18.06
Bayesian R²: 0.836 ± 0.017
LOO-ELPD: -339.97 ± 19.69  (p_loo=24.0) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -7.391 ± 0.428
  Warning: 2 observations with Pareto-k > 0.7 (unreliable LOO estimates)

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_diam1, slope_log_diam2, slope_log_canopy_height, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.765  0.159     4.453      5.068      0.002   
slope_log_diam1          0.397  0.294    -0.189      0.966      0.005   
slope_log_diam2          1.147  0.218     0.738      1.599      0.003   
slope_log_canopy_height  0.370  0.291    -0.183      0.952      0.005   
slope_log_total_height   0.398  0.331    -0.231      1.070      0.005   
slope_log_density        0.110  0.128    -0.141      0.360      0.002   
slope_group              0.583  0.131     0.336      0.846      0.002   
sigma                    0.341  0.040     0.268      0.421      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    4308.0    5360.0    1.0  
slope_log_diam1            0.004    4011.0    4717.0    1.0  
slope_log_diam2            0.003    3985.0    4273.0    1.0  
slope_log_canopy_height    0.003    4119.0    4709.0    1.0  
slope_log_total_height     0.004    4517.0    4837.0    1.0  
slope_log_density          0.001    5694.0    4832.0    1.0  
slope_group                0.001    4625.0    4995.0    1.0  
sigma                      0.001    5092.0    4807.0    1.0  

Regression formula: log_weight = 4.77 + 0.40*log_diam1 + 1.15*log_diam2 + 0.37*log_canopy_height + 0.40*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.016
LOO-ELPD: -19.29 ± 5.22  (p_loo=7.4) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.419 ± 0.113
  Warning: 1 observations with Pareto-k > 0.7 (unreliable LOO estimates)
LOO raw model:              -340.0 ± 19.7
LOO log model (unadjusted): -19.3 ± 5.2
LOO log model (adjusted):   -291.6 ± 7.2

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
log_jacobian	0	-291.589546	7.429382	0.000000	1.000000e+00	7.207880	0.000000	True	log
raw	1	-339.967081	23.965311	48.377535	6.189493e-13	19.693817	17.334181	True	log

ELPD difference is 46, which means that the model with log transformation has a better fit to the data than the model without log transformation.

def posterior_predict(trace, data, predictors):
    post = trace.posterior
    intercept = post["intercept"].values.flatten()
    slopes = np.column_stack([post[f"slope_{p}"].values.flatten() for p in predictors])
    sigma = post["sigma"].values.flatten()
    X = data[predictors].values
    mu = intercept[:, None] + slopes @ X.T      # (n_draws, n_obs)
    return np.random.default_rng(42).normal(mu, sigma[:, None])

preds_1 = ['diam1', 'diam2', 'canopy_height', 'total_height', 'density', 'group']
preds_2 = ['log_diam1', 'log_diam2', 'log_canopy_height', 'log_total_height', 'log_density', 'group']

yrep_1 = posterior_predict(trace_raw, mesquite, preds_1)
yrep_2 = posterior_predict(trace_log, mesquite, preds_2)

n_sims = yrep_1.shape[0]
subset = np.random.default_rng(0).choice(n_sims, 100, replace=False)

# Model 1: raw weight — equivalent of ppc_dens_overlay(mesquite$weight, yrep_1[subset,])
idata_1 = az.from_dict(
    observed_data={"y": mesquite["weight"].values},
    posterior_predictive={"y": yrep_1[np.newaxis, subset, :]}
)
az.plot_ppc(idata_1, num_pp_samples=100, kind="kde")
plt.title("PPC: weight (raw scale)")
plt.show()

# Model 2: log weight — equivalent of ppc_dens_overlay(log(mesquite$weight), yrep_2[subset,])
idata_2 = az.from_dict(
    observed_data={"y": np.log(mesquite["weight"].values)},
    posterior_predictive={"y": yrep_2[np.newaxis, subset, :]}
)
az.plot_ppc(idata_2, num_pp_samples=100, kind="kde")
plt.title("PPC: log(weight)")
plt.show()

Constructing a simpler model

Thinking about a simpler model which makes sense and has good predictive performance by cross validation. Thinking geometrically, leaf weight can be predicted from volume of the leaf canopy:

canopy_volume = diam1 * diam2 * canopy_height

Simplification, as leaves are mostly on outer surface of canopy, not inside.

mesquite['canopy_volume'] = mesquite['diam1'] * mesquite['diam2'] * mesquite['canopy_height']
mesquite['log_canopy_volume'] = np.log(mesquite['canopy_volume'])

fit_and_plot_bayes(mesquite, ['log_canopy_volume'], 'log_weight',
                       intercept_mu=0, intercept_sigma=2000,
                       slope_mu=0, slope_sigma=1000,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       show_residuals=False,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_volume, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                5.169  0.085     5.000      5.334      0.001   
slope_log_canopy_volume  0.722  0.057     0.610      0.833      0.001   
sigma                    0.427  0.047     0.343      0.523      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.001    3722.0    4660.0    1.0  
slope_log_canopy_volume    0.001    3767.0    4540.0    1.0  
sigma                      0.001    5590.0    4961.0    1.0  

Regression formula: log_weight = 5.17 + 0.72*log_canopy_volume
Residual std dev (σ): 0.43 ± 0.05
Bayesian R²: 0.793 ± 0.029
LOO-ELPD: -26.64 ± 4.92  (p_loo=2.5) <- should be lower than number of parameters (3) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.579 ± 0.107

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 208kB
  Dimensions:                  (chain: 4, draw: 2000)
  Coordinates:
    * chain                    (chain) int64 32B 0 1 2 3
    * draw                     (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables:
      intercept                (chain, draw) float64 64kB 5.132 5.135 ... 5.072
      slope_log_canopy_volume  (chain, draw) float64 64kB 0.769 0.6727 ... 0.7685
      sigma                    (chain, draw) float64 64kB 0.4048 0.4599 ... 0.382
  Attributes:
      created_at:                 2026-06-03T22:04:47.975425+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              0.5926599502563477
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (3)
    • intercept
      (chain, draw)
      float64
      5.132 5.135 5.46 ... 5.271 5.072

      array([[5.13196901, 5.13477195, 5.4601399 , ..., 5.15192738, 5.25627606,
              5.17871698],
             [5.20016634, 5.20016634, 5.19505511, ..., 5.16479922, 5.23756642,
              5.2381136 ],
             [5.30844306, 5.13452348, 5.17148595, ..., 5.28808599, 5.27396707,
              5.31303783],
             [5.21335384, 5.25201892, 5.16977605, ..., 5.18963395, 5.27057399,
              5.07236289]], shape=(4, 2000))

    • slope_log_canopy_volume
      (chain, draw)
      float64
      0.769 0.6727 ... 0.6856 0.7685

      array([[0.76896903, 0.67272068, 0.60534824, ..., 0.68985943, 0.656134  ,
              0.72647592],
             [0.74007361, 0.74007361, 0.72684671, ..., 0.67979353, 0.71060495,
              0.72096118],
             [0.68220518, 0.67618305, 0.69762187, ..., 0.67781911, 0.66824784,
              0.6652302 ],
             [0.7293564 , 0.71464978, 0.72512556, ..., 0.6772866 , 0.68556468,
              0.76853123]], shape=(4, 2000))

    • sigma
      (chain, draw)
      float64
      0.4048 0.4599 ... 0.4138 0.382

      array([[0.40483814, 0.45988111, 0.54217081, ..., 0.4211508 , 0.47737689,
              0.3395552 ],
             [0.40258153, 0.40258153, 0.39856298, ..., 0.35326931, 0.38789485,
              0.41770934],
             [0.36929184, 0.38564579, 0.37865184, ..., 0.52313191, 0.47609324,
              0.39311363],
             [0.38204462, 0.51089834, 0.32384276, ..., 0.50828647, 0.41383921,
              0.38200125]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:47.975425+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    0.5926599502563477

    tuning_steps :

    1000

  
  
• log_likelihood

  <xarray.Dataset> Size: 3MB
  Dimensions:  (chain: 4, draw: 2000, y_dim_0: 46)
  Coordinates:
    * chain    (chain) int64 32B 0 1 2 3
    * draw     (draw) int64 16kB 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (chain, draw, y_dim_0) float64 3MB -0.2953 -0.2077 ... -0.7296
  Attributes:
      created_at:                 2026-06-03T22:04:48.112376+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • y_dim_0: 46
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (chain, draw, y_dim_0)
      float64
      -0.2953 -0.2077 ... 0.02744 -0.7296

      array([[[-2.95271290e-01, -2.07743260e-01, -2.07531079e+00, ...,
               -3.98072932e-02, -6.53275543e-02, -5.41267183e-01],
              [-4.66682169e-01, -4.42063974e-01, -2.31682836e+00, ...,
               -2.46160603e-01, -1.63332606e-01, -5.01110965e-01],
              [-3.21841880e-01, -3.25870838e-01, -1.21753548e+00, ...,
               -3.22670180e-01, -3.15654948e-01, -3.10713926e-01],
              ...,
              [-3.81731049e-01, -3.42640804e-01, -2.42462052e+00, ...,
               -1.08311510e-01, -5.80749744e-02, -4.53983457e-01],
              [-3.29023798e-01, -3.20123867e-01, -1.80868824e+00, ...,
               -2.22810282e-01, -1.79813506e-01, -3.32625572e-01],
              [-1.97131500e-01, -1.14734656e-01, -2.87642771e+00, ...,
                1.61123882e-01,  1.36335885e-01, -3.96359907e-01]],
      
             [[-2.11426152e-01, -1.52506010e-01, -1.96235306e+00, ...,
               -2.33069709e-02, -6.05089963e-02, -3.66472726e-01],
              [-2.11426152e-01, -1.52506010e-01, -1.96235306e+00, ...,
               -2.33069709e-02, -6.05089963e-02, -3.66472726e-01],
              [-2.29933299e-01, -1.73595568e-01, -2.11576977e+00, ...,
               -4.51610163e-04, -2.62136957e-02, -3.67873381e-01],
      ...
               -2.73993907e-01, -2.80136091e-01, -3.74127449e-01],
              [-2.98161491e-01, -2.85480000e-01, -1.68328543e+00, ...,
               -1.92602991e-01, -1.79350613e-01, -3.13997783e-01],
              [-1.11719918e-01, -9.75325944e-02, -2.01331838e+00, ...,
                4.89966013e-03,  6.24765850e-04, -1.27424137e-01]],
      
             [[-1.72046549e-01, -1.14170987e-01, -2.13876439e+00, ...,
                3.56089609e-02, -3.69449090e-03, -3.17692682e-01],
              [-3.42398657e-01, -3.18068834e-01, -1.41529971e+00, ...,
               -2.50363445e-01, -2.75485456e-01, -4.00156580e-01],
              [-2.13034355e-01, -1.20428951e-01, -3.21855902e+00, ...,
                2.08120009e-01,  1.89289595e-01, -4.34215591e-01],
              ...,
              [-4.31002591e-01, -4.12038065e-01, -1.80142610e+00, ...,
               -2.80430098e-01, -2.43862294e-01, -4.60571701e-01],
              [-1.84924861e-01, -1.60595173e-01, -1.92667437e+00, ...,
               -3.99309527e-02, -5.15231451e-02, -2.29378669e-01],
              [-4.08606939e-01, -2.89429349e-01, -2.62271959e+00, ...,
                4.03529426e-02,  2.74446298e-02, -7.29586855e-01]]],
            shape=(4, 2000, 46))

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:48.112376+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  
• sample_stats

  <xarray.Dataset> Size: 1MB
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 32B 0 1 2 3
    * draw                   (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables: (12/18)
      energy                 (chain, draw) float64 64kB 46.89 47.81 ... 47.77 47.7
      divergences            (chain, draw) int64 64kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
      acceptance_rate        (chain, draw) float64 64kB 1.0 0.7113 ... 0.9353
      step_size              (chain, draw) float64 64kB 0.7811 0.7811 ... 0.8513
      lp                     (chain, draw) float64 64kB -45.47 -46.78 ... -46.1
      step_size_bar          (chain, draw) float64 64kB 0.7117 0.7117 ... 0.775
      ...                     ...
      tree_depth             (chain, draw) int64 64kB 2 3 2 3 2 3 ... 2 3 3 2 2 3
      energy_error           (chain, draw) float64 64kB -0.2149 ... -0.03314
      perf_counter_start     (chain, draw) float64 64kB 418.2 418.2 ... 418.6
      perf_counter_diff      (chain, draw) float64 64kB 9.192e-05 ... 0.0001731
      index_in_trajectory    (chain, draw) int64 64kB 3 -6 3 2 1 ... -4 -2 -1 3 -2
      reached_max_treedepth  (chain, draw) bool 8kB False False ... False False
  Attributes:
      created_at:                 2026-06-03T22:04:47.989819+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              0.5926599502563477
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (18)
    • energy
      (chain, draw)
      float64
      46.89 47.81 52.53 ... 47.77 47.7

      array([[46.88876842, 47.80693614, 52.5267836 , ..., 47.3908753 ,
              48.20454385, 48.54033218],
             [45.73147293, 47.42159123, 45.44470686, ..., 47.52571551,
              47.07828485, 49.18241518],
             [48.7959182 , 49.6183681 , 46.3378797 , ..., 49.29755565,
              48.00439443, 49.47564081],
             [46.6719428 , 49.6218072 , 48.58171931, ..., 47.52107392,
              47.76858231, 47.70230234]], shape=(4, 2000))

    • divergences
      (chain, draw)
      int64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0]], shape=(4, 2000))

    • acceptance_rate
      (chain, draw)
      float64
      1.0 0.7113 0.8317 ... 0.9891 0.9353

      array([[1.        , 0.71125607, 0.83167097, ..., 0.570391  , 0.62969335,
              0.85104435],
             [0.98369522, 0.69232753, 1.        , ..., 0.86655822, 1.        ,
              0.753364  ],
             [0.83329148, 0.565372  , 1.        , ..., 0.91187003, 1.        ,
              0.59918188],
             [1.        , 0.52972329, 0.98582019, ..., 0.97312849, 0.98909189,
              0.93534622]], shape=(4, 2000))

    • step_size
      (chain, draw)
      float64
      0.7811 0.7811 ... 0.8513 0.8513

      array([[0.78108146, 0.78108146, 0.78108146, ..., 0.78108146, 0.78108146,
              0.78108146],
             [0.8306001 , 0.8306001 , 0.8306001 , ..., 0.8306001 , 0.8306001 ,
              0.8306001 ],
             [0.81460361, 0.81460361, 0.81460361, ..., 0.81460361, 0.81460361,
              0.81460361],
             [0.85131784, 0.85131784, 0.85131784, ..., 0.85131784, 0.85131784,
              0.85131784]], shape=(4, 2000))

    • lp
      (chain, draw)
      float64
      -45.47 -46.78 ... -45.83 -46.1

      array([[-45.47259069, -46.78105517, -51.61329079, ..., -45.62027538,
              -46.61633781, -46.90101471],
             [-45.47735983, -45.47735983, -45.23714303, ..., -47.01435631,
              -45.70644478, -45.67947855],
             [-47.51327944, -46.71778526, -45.57364365, ..., -48.02705732,
              -46.60525173, -46.80565534],
             [-45.71411919, -47.47719244, -47.98465548, ..., -47.19233586,
              -45.83056587, -46.10067902]], shape=(4, 2000))

    • step_size_bar
      (chain, draw)
      float64
      0.7117 0.7117 ... 0.775 0.775

      array([[0.71172417, 0.71172417, 0.71172417, ..., 0.71172417, 0.71172417,
              0.71172417],
             [0.77556943, 0.77556943, 0.77556943, ..., 0.77556943, 0.77556943,
              0.77556943],
             [0.72321916, 0.72321916, 0.72321916, ..., 0.72321916, 0.72321916,
              0.72321916],
             [0.77497042, 0.77497042, 0.77497042, ..., 0.77497042, 0.77497042,
              0.77497042]], shape=(4, 2000))

    • max_energy_error
      (chain, draw)
      float64
      -0.2149 0.8198 ... -0.2147 0.1733

      array([[-0.21490088,  0.81979238,  0.26346691, ...,  1.08691702,
               0.75610271,  0.49151977],
             [-0.08596809,  0.61640038, -0.13006729, ...,  0.43854523,
              -0.54020693,  1.34671984],
             [ 0.68261942,  0.65803623, -0.60562322, ..., -0.88988897,
              -0.32264209,  1.28643029],
             [-0.24391223,  1.73623026, -0.39280883, ..., -0.25816243,
              -0.21465838,  0.17331433]], shape=(4, 2000))

    • n_steps
      (chain, draw)
      float64
      3.0 7.0 3.0 7.0 ... 5.0 3.0 3.0 7.0

      array([[3., 7., 3., ..., 3., 3., 5.],
             [3., 3., 1., ..., 3., 3., 3.],
             [7., 3., 1., ..., 7., 3., 7.],
             [7., 5., 7., ..., 3., 3., 7.]], shape=(4, 2000))

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

    • process_time_diff
      (chain, draw)
      float64
      9.3e-05 0.000174 ... 0.000174

      array([[9.30e-05, 1.74e-04, 8.90e-05, ..., 8.30e-05, 8.50e-05, 1.25e-04],
             [8.70e-05, 8.50e-05, 4.00e-05, ..., 9.30e-05, 9.10e-05, 9.00e-05],
             [1.79e-04, 9.40e-05, 4.50e-05, ..., 1.76e-04, 9.00e-05, 1.80e-04],
             [1.81e-04, 1.36e-04, 1.82e-04, ..., 8.60e-05, 8.70e-05, 1.74e-04]],
            shape=(4, 2000))

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • tree_depth
      (chain, draw)
      int64
      2 3 2 3 2 3 2 2 ... 3 3 2 3 3 2 2 3

      array([[2, 3, 2, ..., 2, 2, 3],
             [2, 2, 1, ..., 2, 2, 2],
             [3, 2, 1, ..., 3, 2, 3],
             [3, 3, 3, ..., 2, 2, 3]], shape=(4, 2000))

    • energy_error
      (chain, draw)
      float64
      -0.2149 0.5078 ... -0.2147 -0.03314

      array([[-0.21490088,  0.50784313,  0.26346691, ...,  0.        ,
               0.09564552, -0.05380515],
             [ 0.03824226,  0.        , -0.13006729, ...,  0.43854523,
              -0.54020693, -0.36192205],
             [-0.29473945,  0.62300024, -0.60562322, ..., -0.0827999 ,
              -0.16517807, -0.32359728],
             [-0.06260109,  0.5066156 , -0.19588436, ...,  0.07930778,
              -0.21465838, -0.03313927]], shape=(4, 2000))

    • perf_counter_start
      (chain, draw)
      float64
      418.2 418.2 418.2 ... 418.6 418.6

      array([[418.23584329, 418.23597783, 418.23618787, ..., 418.55730046,
              418.5574145 , 418.55753221],
             [418.245715  , 418.24584787, 418.24596738, ..., 418.5668305 ,
              418.56695679, 418.567084  ],
             [418.24136267, 418.24159087, 418.24171925, ..., 418.57444242,
              418.57465308, 418.57477633],
             [418.25018262, 418.25041979, 418.25070417, ..., 418.60017329,
              418.60029108, 418.60040925]], shape=(4, 2000))

    • perf_counter_diff
      (chain, draw)
      float64
      9.192e-05 0.0001742 ... 0.0001731

      array([[9.19170e-05, 1.74167e-04, 8.84160e-05, ..., 8.22080e-05,
              8.49160e-05, 1.25958e-04],
             [8.75830e-05, 8.45000e-05, 4.02080e-05, ..., 9.19160e-05,
              9.15840e-05, 9.02080e-05],
             [1.79584e-04, 9.25410e-05, 4.47910e-05, ..., 1.75667e-04,
              8.95000e-05, 1.80042e-04],
             [1.81500e-04, 1.35125e-04, 1.82625e-04, ..., 8.51250e-05,
              8.59580e-05, 1.73083e-04]], shape=(4, 2000))

    • index_in_trajectory
      (chain, draw)
      int64
      3 -6 3 2 1 6 -1 ... 1 -4 -2 -1 3 -2

      array([[ 3, -6,  3, ...,  0, -2, -2],
             [ 2,  0,  1, ...,  3, -2,  3],
             [ 3, -2, -1, ...,  1, -1,  2],
             [ 2,  2, -3, ..., -1,  3, -2]], shape=(4, 2000))

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:47.989819+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    0.5926599502563477

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset> Size: 736B
  Dimensions:  (y_dim_0: 46)
  Coordinates:
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (y_dim_0) float64 368B 5.995 6.242 7.073 5.73 ... 7.465 6.811 5.364
  Attributes:
      created_at:                 2026-06-03T22:04:47.993748+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • y_dim_0: 46
  • Coordinates: (1)
    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (y_dim_0)
      float64
      5.995 6.242 7.073 ... 6.811 5.364

      array([5.99470928, 6.24163944, 7.07259152, 5.73009978, 6.75133536,
             5.59359552, 5.04664573, 7.13345556, 5.79301361, 6.42097165,
             4.09767235, 5.59693938, 6.10568569, 4.79081953, 5.93674433,
             5.58499894, 4.93375425, 6.92834191, 6.45472675, 6.4326185 ,
             6.36268322, 5.78934836, 4.20020495, 4.21950771, 5.0310913 ,
             5.54673873, 6.58340922, 8.30696587, 5.84354442, 5.80181621,
             5.09681299, 7.05617528, 5.95739057, 6.54175124, 6.51382341,
             5.38219885, 6.6480774 , 5.83393316, 4.83389812, 6.13664656,
             4.16666522, 6.74594198, 5.420535  , 7.46456714, 6.81124438,
             5.36363683])

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:47.993748+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  

Model has 0.7*log_canopy_volume, which suggests that this is attenuating the relationship between canopy volume and leaf weight, which could be due to the fact that leaves are mostly on the outer surface of the canopy, rather than being distributed throughout the entire volume.

# Re-fit both models (no plots) to get traces for comparison

trace_log1 = fit_and_plot_bayes(
    mesquite, ['log_diam1', 'log_diam2', 'log_canopy_height', 'log_total_height', 'log_density', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

trace_log2 = fit_and_plot_bayes(
    mesquite, ['log_canopy_volume'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

# Pointwise LOO for both models
loo_log1 = az.loo(trace_log1, pointwise=True)
loo_log2 = az.loo(trace_log2, pointwise=True)

print(f"LOO log1 model:              {loo_log1.elpd_loo:.1f} ± {loo_log1.se:.1f}")
print(f"LOO log2 model: {loo_log2.elpd_loo:.1f} ± {loo_log2.se:.1f}")

az.compare({"log1": loo_log1, "log2": loo_log2})

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_diam1, slope_log_diam2, slope_log_canopy_height, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.
Initializing NUTS using jitter+adapt_diag...

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.771  0.161     4.434      5.066      0.003   
slope_log_diam1          0.391  0.290    -0.170      0.947      0.005   
slope_log_diam2          1.153  0.218     0.737      1.588      0.003   
slope_log_canopy_height  0.373  0.293    -0.219      0.927      0.004   
slope_log_total_height   0.396  0.331    -0.229      1.078      0.005   
slope_log_density        0.107  0.128    -0.141      0.353      0.002   
slope_group              0.581  0.136     0.306      0.838      0.002   
sigma                    0.340  0.040     0.265      0.416      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    4045.0    4626.0    1.0  
slope_log_diam1            0.003    3420.0    4692.0    1.0  
slope_log_diam2            0.002    3923.0    5300.0    1.0  
slope_log_canopy_height    0.003    4665.0    5013.0    1.0  
slope_log_total_height     0.004    4780.0    5491.0    1.0  
slope_log_density          0.001    5748.0    5472.0    1.0  
slope_group                0.002    4717.0    4418.0    1.0  
sigma                      0.000    5644.0    5303.0    1.0  

Regression formula: log_weight = 4.77 + 0.39*log_diam1 + 1.15*log_diam2 + 0.37*log_canopy_height + 0.40*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.017
LOO-ELPD: -19.38 ± 5.24  (p_loo=7.5) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.421 ± 0.114

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_volume, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                5.169  0.087     4.996      5.339      0.001   
slope_log_canopy_volume  0.723  0.057     0.615      0.834      0.001   
sigma                    0.427  0.048     0.340      0.525      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.001    4469.0    4673.0    1.0  
slope_log_canopy_volume    0.001    3987.0    4588.0    1.0  
sigma                      0.001    5015.0    4170.0    1.0  

Regression formula: log_weight = 5.17 + 0.72*log_canopy_volume
Residual std dev (σ): 0.43 ± 0.05
Bayesian R²: 0.793 ± 0.028
LOO-ELPD: -26.72 ± 4.96  (p_loo=2.6) <- should be lower than number of parameters (3) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.581 ± 0.108
LOO log1 model:              -19.4 ± 5.2
LOO log2 model: -26.7 ± 5.0

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
log1	0	-19.378336	7.508495	0.000000	0.839799	5.236462	0.000000	False	log
log2	1	-26.724994	2.643242	7.346658	0.160201	4.959957	4.885636	False	log

The simpler model has worse elpd_loo

We now try canopy area and canopy shape

mesquite['log_canopy_area'] = np.log(mesquite['diam1'] * mesquite['diam2'])
mesquite['log_canopy_shape'] = np.log(mesquite['diam1'] / mesquite['diam2'])

fit_and_plot_bayes(mesquite, ['log_canopy_area', 'log_canopy_shape', 'log_canopy_volume', 'log_total_height', 'log_density', 'group'], 'log_weight',
                       intercept_mu=0, intercept_sigma=2000,
                       slope_mu=0, slope_sigma=1000,
                       sigma_sigma=50,
                       samples=2000, tune=1000, hdi_prob=0.95,
                       show_trace=True, show_forest=False,
                       show_posterior=False, show_regression=True,
                       show_residuals=False,
                       n_regression_lines=100)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_area, slope_log_canopy_shape, slope_log_canopy_volume, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.765  0.163     4.445      5.090      0.002   
slope_log_canopy_area    0.403  0.309    -0.220      0.990      0.006   
slope_log_canopy_shape  -0.375  0.241    -0.848      0.109      0.003   
slope_log_canopy_volume  0.370  0.295    -0.187      0.973      0.005   
slope_log_total_height   0.399  0.333    -0.272      1.044      0.005   
slope_log_density        0.108  0.126    -0.137      0.356      0.002   
slope_group              0.584  0.134     0.329      0.857      0.002   
sigma                    0.341  0.040     0.266      0.422      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    4296.0    5050.0    1.0  
slope_log_canopy_area      0.004    3123.0    4144.0    1.0  
slope_log_canopy_shape     0.003    5015.0    5200.0    1.0  
slope_log_canopy_volume    0.004    3111.0    3887.0    1.0  
slope_log_total_height     0.004    4158.0    4759.0    1.0  
slope_log_density          0.001    6285.0    6027.0    1.0  
slope_group                0.002    4936.0    4847.0    1.0  
sigma                      0.001    5010.0    4610.0    1.0  

Regression formula: log_weight = 4.76 + 0.40*log_canopy_area + -0.37*log_canopy_shape + 0.37*log_canopy_volume + 0.40*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.017
LOO-ELPD: -19.58 ± 5.26  (p_loo=7.7) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.426 ± 0.114
  Warning: 1 observations with Pareto-k > 0.7 (unreliable LOO estimates)

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 528kB
  Dimensions:                  (chain: 4, draw: 2000)
  Coordinates:
    * chain                    (chain) int64 32B 0 1 2 3
    * draw                     (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables:
      intercept                (chain, draw) float64 64kB 4.756 4.884 ... 4.55
      slope_log_canopy_area    (chain, draw) float64 64kB 0.3663 0.1098 ... 0.5732
      slope_log_canopy_shape   (chain, draw) float64 64kB -0.6298 ... -0.2926
      slope_log_canopy_volume  (chain, draw) float64 64kB 0.3932 0.7007 ... 0.1585
      slope_log_total_height   (chain, draw) float64 64kB 0.5352 ... 0.5985
      slope_log_density        (chain, draw) float64 64kB 0.05741 ... 0.3869
      slope_group              (chain, draw) float64 64kB 0.6184 0.7657 ... 0.7324
      sigma                    (chain, draw) float64 64kB 0.3303 0.3589 ... 0.299
  Attributes:
      created_at:                 2026-06-03T22:04:57.708457+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              3.8649630546569824
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (8)
    • intercept
      (chain, draw)
      float64
      4.756 4.884 4.895 ... 4.534 4.55

      array([[4.75623079, 4.88411945, 4.89498925, ..., 4.95637568, 4.9987286 ,
              4.45075825],
             [4.62509649, 4.56127749, 4.42068985, ..., 4.6962073 , 4.70353659,
              4.57694057],
             [4.87659435, 4.63164228, 5.04788887, ..., 4.59132456, 4.44660006,
              4.74161267],
             [4.93162378, 4.78015813, 4.81396142, ..., 4.91776187, 4.53404384,
              4.54969635]], shape=(4, 2000))

    • slope_log_canopy_area
      (chain, draw)
      float64
      0.3663 0.1098 ... 0.2007 0.5732

      array([[ 0.36628541,  0.10983165,  0.19326443, ..., -0.03698162,
               0.0754949 ,  1.00555461],
             [ 0.87756904,  0.60822438,  0.43208407, ...,  1.00857891,
               0.93739606,  0.79077022],
             [ 0.48748953,  0.22855157,  0.06716611, ...,  1.00535041,
               1.0313449 ,  0.41533934],
             [-0.48213413,  0.23007836,  0.15148185, ...,  0.04040802,
               0.20074646,  0.57315109]], shape=(4, 2000))

    • slope_log_canopy_shape
      (chain, draw)
      float64
      -0.6298 -0.8584 ... -0.4484 -0.2926

      array([[-0.62975835, -0.85838619, -0.83136246, ..., -0.64783813,
              -0.32354179, -0.40073782],
             [ 0.01437332, -0.35595706, -0.23842369, ..., -0.47784577,
              -0.51814467, -0.25308117],
             [-0.27107278, -0.59093294, -0.78600889, ...,  0.00547612,
              -0.118891  , -0.29819072],
             [-0.38990093, -0.25103373, -0.2870315 , ..., -0.58398707,
              -0.44842914, -0.29263084]], shape=(4, 2000))

    • slope_log_canopy_volume
      (chain, draw)
      float64
      0.3932 0.7007 ... 0.5411 0.1585

      array([[ 0.39322363,  0.70070647,  0.7030344 , ...,  0.72824935,
               0.66285157, -0.1544669 ],
             [-0.07585497,  0.20909888,  0.40063345, ..., -0.34824881,
              -0.24001527,  0.00249153],
             [ 0.31143858,  0.59137249,  0.59754661, ..., -0.14468465,
              -0.16850479,  0.3463051 ],
             [ 1.02614444,  0.74710441,  0.83236079, ...,  0.74166965,
               0.54107778,  0.15854792]], shape=(4, 2000))

    • slope_log_total_height
      (chain, draw)
      float64
      0.5352 0.09935 ... 0.3358 0.5985

      array([[ 0.5352357 ,  0.09935448,  0.02258294, ...,  0.29581566,
               0.03044768,  0.92280966],
             [ 0.83003371,  0.32420524,  0.57335816, ...,  0.95889208,
               1.03067915,  0.62444754],
             [ 0.17138871,  0.35989312,  0.0956792 , ...,  0.84774148,
               0.9552322 ,  0.3594651 ],
             [ 0.38637115, -0.70287744, -0.66795924, ..., -0.34179871,
               0.33577202,  0.59847631]], shape=(4, 2000))

    • slope_log_density
      (chain, draw)
      float64
      0.05741 -0.01346 ... 0.2451 0.3869

      array([[ 0.05741241, -0.01346399,  0.04157851, ...,  0.16416352,
              -0.0921318 ,  0.20848105],
             [-0.15512424,  0.22847919,  0.24374818, ...,  0.20201291,
               0.09571221,  0.11330718],
             [-0.0331365 ,  0.16243205,  0.38494341, ..., -0.09603659,
               0.03131858,  0.12967982],
             [ 0.02762345,  0.17328511,  0.21319595, ...,  0.31264277,
               0.24509673,  0.3869132 ]], shape=(4, 2000))

    • slope_group
      (chain, draw)
      float64
      0.6184 0.7657 ... 0.9041 0.7324

      array([[0.6184468 , 0.76567707, 0.74142698, ..., 0.64466397, 0.42251035,
              0.66080097],
             [0.51783091, 0.62131029, 0.74353186, ..., 0.64686404, 0.61481666,
              0.61467089],
             [0.48527913, 0.79851048, 0.59561764, ..., 0.49447664, 0.52130315,
              0.58581599],
             [0.65504267, 0.68258767, 0.64262906, ..., 0.72994692, 0.90408416,
              0.7324414 ]], shape=(4, 2000))

    • sigma
      (chain, draw)
      float64
      0.3303 0.3589 ... 0.3645 0.299

      array([[0.33028843, 0.35893895, 0.35545728, ..., 0.27514757, 0.36931403,
              0.29257235],
             [0.42429954, 0.36149648, 0.30414333, ..., 0.28920007, 0.27215218,
              0.29140074],
             [0.2874683 , 0.35925727, 0.39265696, ..., 0.38611846, 0.3387074 ,
              0.32514886],
             [0.36467607, 0.4401686 , 0.42991332, ..., 0.36076744, 0.36453829,
              0.29900791]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:57.708457+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    3.8649630546569824

    tuning_steps :

    1000

  
  
• log_likelihood

  <xarray.Dataset> Size: 3MB
  Dimensions:  (chain: 4, draw: 2000, y_dim_0: 46)
  Coordinates:
    * chain    (chain) int64 32B 0 1 2 3
    * draw     (draw) int64 16kB 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (chain, draw, y_dim_0) float64 3MB -0.01169 0.135 ... -3.608
  Attributes:
      created_at:                 2026-06-03T22:04:58.089608+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • y_dim_0: 46
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (chain, draw, y_dim_0)
      float64
      -0.01169 0.135 ... -0.3383 -3.608

      array([[[-1.16856751e-02,  1.34972454e-01,  1.76210150e-01, ...,
                1.42690403e-01,  1.68277028e-01, -2.81084622e+00],
              [ 5.54778552e-02,  5.58788235e-02, -3.62345741e-02, ...,
                8.30723013e-02,  7.78085661e-02, -2.59499763e+00],
              [ 9.44744886e-02,  2.03619149e-02,  8.05939584e-02, ...,
                1.66157432e-02,  1.10137807e-01, -2.36110821e+00],
              ...,
              [ 3.00931397e-01,  3.28418618e-01, -2.94289525e-01, ...,
                3.58080291e-01,  3.51802340e-01, -3.54098832e+00],
              [-3.28687745e-02,  5.16423850e-02, -1.01324232e+00, ...,
               -1.28610118e-01,  2.46824603e-02, -1.17810939e+00],
              [ 7.32427863e-02, -1.56207984e-01, -9.47763828e-01, ...,
                3.08904047e-01,  2.68104417e-01, -3.50664540e+00]],
      
             [[-6.70697046e-02, -1.74862872e-01, -3.56032202e-01, ...,
               -1.72522445e-01, -6.16245401e-02, -9.84444558e-01],
              [-2.27437744e-01, -2.57841276e-01, -1.70260934e-02, ...,
               -8.34421868e-02, -3.95831988e-01, -2.46833406e+00],
              [ 9.41665979e-02,  1.42636572e-01,  2.53823126e-01, ...,
                2.59608157e-01,  9.09546227e-02, -5.57750300e+00],
      ...
                2.52417184e-02,  1.84058511e-02, -1.04479929e+00],
              [-1.36212156e-02, -3.57293183e-01, -5.40633221e-01, ...,
                1.55936857e-01,  1.62792413e-01, -2.21210605e+00],
              [ 9.46418505e-02,  1.07020245e-01,  5.81161992e-02, ...,
                6.02214884e-02,  2.00649645e-02, -2.15791094e+00]],
      
             [[ 6.22815205e-02,  3.95131283e-02, -1.64692973e+00, ...,
               -4.13003217e-01,  1.12931494e-02, -3.07810413e+00],
              [-1.34210206e-01, -9.96900746e-02, -1.14022544e+00, ...,
               -1.16685996e-01, -5.30649586e-01, -1.39359008e+00],
              [-1.17832346e-01, -7.47913425e-02, -1.21241892e+00, ...,
               -8.19361450e-02, -3.25609653e-01, -1.49800541e+00],
              ...,
              [ 3.63397086e-02,  9.71078915e-02, -1.21073361e+00, ...,
               -1.60068112e-03, -5.34480558e-01, -1.91629410e+00],
              [ 2.69096893e-02,  8.54282834e-02, -1.19384349e-01, ...,
               -2.61070870e-01, -5.12930034e-01, -4.02742620e+00],
              [ 1.55177467e-01,  7.72304832e-02,  2.88345862e-01, ...,
                1.06012483e-01, -3.38266829e-01, -3.60797043e+00]]],
            shape=(4, 2000, 46))

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:58.089608+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  
• sample_stats

  <xarray.Dataset> Size: 1MB
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 32B 0 1 2 3
    * draw                   (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
  Data variables: (12/18)
      energy                 (chain, draw) float64 64kB 79.19 80.82 ... 82.83
      divergences            (chain, draw) int64 64kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
      acceptance_rate        (chain, draw) float64 64kB 0.9858 0.9807 ... 0.6856
      step_size              (chain, draw) float64 64kB 0.121 0.121 ... 0.08538
      lp                     (chain, draw) float64 64kB -73.5 -76.51 ... -75.92
      step_size_bar          (chain, draw) float64 64kB 0.1227 0.1227 ... 0.1034
      ...                     ...
      tree_depth             (chain, draw) int64 64kB 5 6 3 6 5 5 ... 5 5 6 5 6 6
      energy_error           (chain, draw) float64 64kB -0.06507 ... -1.183
      perf_counter_start     (chain, draw) float64 64kB 425.7 425.7 ... 428.3
      perf_counter_diff      (chain, draw) float64 64kB 0.0007213 ... 0.00198
      index_in_trajectory    (chain, draw) int64 64kB -8 -11 -5 -12 ... 24 12 -9
      reached_max_treedepth  (chain, draw) bool 8kB False False ... False False
  Attributes:
      created_at:                 2026-06-03T22:04:57.722565+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2
      sampling_time:              3.8649630546569824
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

  • Data variables: (18)
    • energy
      (chain, draw)
      float64
      79.19 80.82 79.6 ... 80.37 82.83

      array([[79.19405946, 80.81598537, 79.59863109, ..., 77.97586877,
              80.25588914, 78.99859439],
             [81.58792298, 83.73190616, 80.73722115, ..., 82.67168214,
              84.78822542, 83.063827  ],
             [79.42337077, 78.72467581, 79.21654077, ..., 82.10781806,
              78.36852877, 76.68796065],
             [80.64908224, 82.74469336, 82.5917984 , ..., 80.67180674,
              80.37014435, 82.832048  ]], shape=(4, 2000))

    • divergences
      (chain, draw)
      int64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 1, 1, 1],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0],
             [0, 0, 0, ..., 0, 0, 0]], shape=(4, 2000))

    • acceptance_rate
      (chain, draw)
      float64
      0.9858 0.9807 ... 0.844 0.6856

      array([[0.98580872, 0.98065277, 0.60088295, ..., 0.25301951, 0.7473732 ,
              0.95734364],
             [0.76415309, 0.71863369, 0.74051748, ..., 0.90387882, 0.62724543,
              0.58738558],
             [0.32475876, 0.99304462, 0.99263525, ..., 0.68275616, 0.99361473,
              0.89724562],
             [0.99842207, 0.99985222, 0.96866929, ..., 0.89490533, 0.84396056,
              0.68563606]], shape=(4, 2000))

    • step_size
      (chain, draw)
      float64
      0.121 0.121 ... 0.08538 0.08538

      array([[0.12099286, 0.12099286, 0.12099286, ..., 0.12099286, 0.12099286,
              0.12099286],
             [0.07970563, 0.07970563, 0.07970563, ..., 0.07970563, 0.07970563,
              0.07970563],
             [0.12778862, 0.12778862, 0.12778862, ..., 0.12778862, 0.12778862,
              0.12778862],
             [0.08537906, 0.08537906, 0.08537906, ..., 0.08537906, 0.08537906,
              0.08537906]], shape=(4, 2000))

    • lp
      (chain, draw)
      float64
      -73.5 -76.51 ... -76.35 -75.92

      array([[-73.50192705, -76.5082761 , -77.39719858, ..., -75.28780034,
              -76.53733815, -77.57266624],
             [-78.92253988, -76.07483281, -76.62305826, ..., -78.92297641,
              -77.4587571 , -73.87711683],
             [-75.66539865, -74.83895643, -76.94378367, ..., -76.76170331,
              -75.01350691, -71.71818084],
             [-79.0559429 , -80.60634885, -79.50449406, ..., -75.74291427,
              -76.35004725, -75.91646661]], shape=(4, 2000))

    • step_size_bar
      (chain, draw)
      float64
      0.1227 0.1227 ... 0.1034 0.1034

      array([[0.12274484, 0.12274484, 0.12274484, ..., 0.12274484, 0.12274484,
              0.12274484],
             [0.11388639, 0.11388639, 0.11388639, ..., 0.11388639, 0.11388639,
              0.11388639],
             [0.10594491, 0.10594491, 0.10594491, ..., 0.10594491, 0.10594491,
              0.10594491],
             [0.10335405, 0.10335405, 0.10335405, ..., 0.10335405, 0.10335405,
              0.10335405]], shape=(4, 2000))

    • max_energy_error
      (chain, draw)
      float64
      -0.1507 -0.08187 ... 0.5013 2.545

      array([[-0.15073819, -0.0818662 ,  1.02805655, ...,  5.05013953,
               1.75255003, -1.17180954],
             [ 0.63240012,  1.99735988,  1.43449799, ..., -0.46489129,
               1.57477477,  1.30745466],
             [ 2.54296196, -0.77416644, -0.12476714, ...,  0.91387707,
              -0.11709142,  0.26239563],
             [-0.9993915 , -0.33133847, -0.20581177, ...,  0.28351153,
               0.50125641,  2.54479043]], shape=(4, 2000))

    • n_steps
      (chain, draw)
      float64
      23.0 63.0 7.0 ... 31.0 63.0 63.0

      array([[23., 63.,  7., ..., 31., 31., 63.],
             [63., 63., 15., ..., 15., 31., 15.],
             [31., 31., 31., ..., 31., 31., 63.],
             [15., 31., 31., ..., 31., 63., 63.]], shape=(4, 2000))

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

    • process_time_diff
      (chain, draw)
      float64
      0.000721 0.001907 ... 0.001949

      array([[0.000721, 0.001907, 0.000248, ..., 0.001029, 0.000979, 0.001938],
             [0.001911, 0.0019  , 0.000486, ..., 0.000466, 0.000933, 0.000479],
             [0.000933, 0.000938, 0.001002, ..., 0.000929, 0.00094 , 0.001834],
             [0.000524, 0.001008, 0.001017, ..., 0.00109 , 0.001962, 0.001949]],
            shape=(4, 2000))

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

    • tree_depth
      (chain, draw)
      int64
      5 6 3 6 5 5 5 5 ... 5 5 5 5 6 5 6 6

      array([[5, 6, 3, ..., 5, 5, 6],
             [6, 6, 4, ..., 4, 5, 4],
             [5, 5, 5, ..., 5, 5, 6],
             [4, 5, 5, ..., 5, 6, 6]], shape=(4, 2000))

    • energy_error
      (chain, draw)
      float64
      -0.06507 -0.01839 ... 0.3537 -1.183

      array([[-0.06507017, -0.01839356,  0.76569094, ..., -0.12351911,
               0.57778413, -0.69099196],
             [ 0.27962637,  0.28570437, -1.24701137, ...,  0.32595661,
              -0.8750935 ,  0.65444622],
             [ 0.76285807, -0.38612548, -0.00316673, ...,  0.14239423,
              -0.07677456,  0.0575101 ],
             [-0.57299179, -0.18230886, -0.06542098, ..., -0.05504579,
               0.35372269, -1.1828889 ]], shape=(4, 2000))

    • perf_counter_start
      (chain, draw)
      float64
      425.7 425.7 425.7 ... 428.3 428.3

      array([[425.74015721, 425.74093887, 425.742893  , ..., 427.91761321,
              427.91869242, 427.91971792],
             [425.74730979, 425.74927221, 425.75122008, ..., 427.91333104,
              427.91383712, 427.914813  ],
             [425.80491092, 425.80589617, 425.80687767, ..., 428.18476233,
              428.18573537, 428.18671817],
             [425.81819133, 425.81877579, 425.81983325, ..., 428.32673042,
              428.32792012, 428.32994529]], shape=(4, 2000))

    • perf_counter_diff
      (chain, draw)
      float64
      0.0007213 0.001906 ... 0.00198

      array([[0.00072129, 0.00190633, 0.0002475 , ..., 0.00102871, 0.00097983,
              0.00193746],
             [0.00191067, 0.00190083, 0.00048558, ..., 0.00046575, 0.00093408,
              0.00047804],
             [0.00093287, 0.00093808, 0.00102704, ..., 0.00093008, 0.00094062,
              0.00183413],
             [0.00052392, 0.00100758, 0.00103   , ..., 0.00110583, 0.00196271,
              0.00198029]], shape=(4, 2000))

    • index_in_trajectory
      (chain, draw)
      int64
      -8 -11 -5 -12 11 ... 27 14 24 12 -9

      array([[ -8, -11,  -5, ...,  -9,  13, -25],
             [ 23, -12,  14, ...,   9,   3,  11],
             [  2,  23,  16, ...,  10,   8, -20],
             [  9,  21,  -3, ...,  24,  12,  -9]], shape=(4, 2000))

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • Attributes: (6)

    created_at :

    2026-06-03T22:04:57.722565+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

    sampling_time :

    3.8649630546569824

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset> Size: 736B
  Dimensions:  (y_dim_0: 46)
  Coordinates:
    * y_dim_0  (y_dim_0) int64 368B 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
  Data variables:
      y        (y_dim_0) float64 368B 5.995 6.242 7.073 5.73 ... 7.465 6.811 5.364
  Attributes:
      created_at:                 2026-06-03T22:04:57.726634+00:00
      arviz_version:              0.23.4
      inference_library:          pymc
      inference_library_version:  5.28.2

  xarray.Dataset

  • Dimensions:
    • y_dim_0: 46
  • Coordinates: (1)
    • y_dim_0
      (y_dim_0)
      int64
      0 1 2 3 4 5 6 ... 40 41 42 43 44 45

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • Data variables: (1)
    • y
      (y_dim_0)
      float64
      5.995 6.242 7.073 ... 6.811 5.364

      array([5.99470928, 6.24163944, 7.07259152, 5.73009978, 6.75133536,
             5.59359552, 5.04664573, 7.13345556, 5.79301361, 6.42097165,
             4.09767235, 5.59693938, 6.10568569, 4.79081953, 5.93674433,
             5.58499894, 4.93375425, 6.92834191, 6.45472675, 6.4326185 ,
             6.36268322, 5.78934836, 4.20020495, 4.21950771, 5.0310913 ,
             5.54673873, 6.58340922, 8.30696587, 5.84354442, 5.80181621,
             5.09681299, 7.05617528, 5.95739057, 6.54175124, 6.51382341,
             5.38219885, 6.6480774 , 5.83393316, 4.83389812, 6.13664656,
             4.16666522, 6.74594198, 5.420535  , 7.46456714, 6.81124438,
             5.36363683])

  • Attributes: (4)

    created_at :

    2026-06-03T22:04:57.726634+00:00

    arviz_version :

    0.23.4

    inference_library :

    pymc

    inference_library_version :

    5.28.2

  
  

# Re-fit both models (no plots) to get traces for comparison

trace_log4 = fit_and_plot_bayes(
    mesquite, ['log_diam1', 'log_diam2', 'log_canopy_height', 'log_total_height', 'log_density', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

trace_log3 = fit_and_plot_bayes(
    mesquite, ['log_canopy_area', 'log_canopy_shape', 'log_canopy_volume', 'log_total_height', 'log_density', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)


# Pointwise LOO for both models
loo_log4 = az.loo(trace_log4, pointwise=True)
loo_log3 = az.loo(trace_log3, pointwise=True)

print(f"LOO log4 model: {loo_log4.elpd_loo:.1f} ± {loo_log4.se:.1f}")
print(f"LOO log3 model: {loo_log3.elpd_loo:.1f} ± {loo_log3.se:.1f}")

az.compare({"log4": loo_log4, "log3": loo_log3})

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_diam1, slope_log_diam2, slope_log_canopy_height, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.770  0.161     4.453      5.086      0.002   
slope_log_diam1          0.393  0.292    -0.154      0.981      0.005   
slope_log_diam2          1.151  0.220     0.714      1.579      0.003   
slope_log_canopy_height  0.375  0.291    -0.208      0.931      0.004   
slope_log_total_height   0.392  0.324    -0.221      1.043      0.005   
slope_log_density        0.109  0.126    -0.145      0.353      0.002   
slope_group              0.581  0.135     0.313      0.837      0.002   
sigma                    0.340  0.040     0.263      0.418      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    4653.0    4929.0    1.0  
slope_log_diam1            0.003    3795.0    4639.0    1.0  
slope_log_diam2            0.002    3986.0    5022.0    1.0  
slope_log_canopy_height    0.003    5026.0    5314.0    1.0  
slope_log_total_height     0.004    5143.0    4905.0    1.0  
slope_log_density          0.002    5318.0    5429.0    1.0  
slope_group                0.002    5056.0    4661.0    1.0  
sigma                      0.001    4733.0    4727.0    1.0  

Regression formula: log_weight = 4.77 + 0.39*log_diam1 + 1.15*log_diam2 + 0.37*log_canopy_height + 0.39*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.017
LOO-ELPD: -19.42 ± 5.25  (p_loo=7.6) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.422 ± 0.114
  Warning: 1 observations with Pareto-k > 0.7 (unreliable LOO estimates)

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_area, slope_log_canopy_shape, slope_log_canopy_volume, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.765  0.163     4.444      5.086      0.003   
slope_log_canopy_area    0.399  0.302    -0.185      0.999      0.006   
slope_log_canopy_shape  -0.377  0.240    -0.840      0.102      0.004   
slope_log_canopy_volume  0.373  0.289    -0.170      0.941      0.006   
slope_log_total_height   0.398  0.319    -0.247      1.006      0.005   
slope_log_density        0.109  0.129    -0.163      0.342      0.002   
slope_group              0.586  0.133     0.325      0.843      0.002   
sigma                    0.340  0.040     0.268      0.421      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    3753.0    4616.0    1.0  
slope_log_canopy_area      0.004    2607.0    4014.0    1.0  
slope_log_canopy_shape     0.003    3911.0    5151.0    1.0  
slope_log_canopy_volume    0.003    2695.0    3899.0    1.0  
slope_log_total_height     0.003    3793.0    5072.0    1.0  
slope_log_density          0.002    5193.0    4765.0    1.0  
slope_group                0.002    4951.0    5426.0    1.0  
sigma                      0.000    5260.0    4840.0    1.0  

Regression formula: log_weight = 4.76 + 0.40*log_canopy_area + -0.38*log_canopy_shape + 0.37*log_canopy_volume + 0.40*log_total_height + 0.11*log_density + 0.59*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.017
LOO-ELPD: -19.43 ± 5.29  (p_loo=7.6) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.422 ± 0.115
LOO log4 model: -19.4 ± 5.3
LOO log3 model: -19.4 ± 5.3

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
log4	0	-19.424221	7.569362	0.000000	0.918638	5.250184	0.000000	True	log
log3	1	-19.431993	7.613154	0.007771	0.081362	5.286802	0.136447	False	log

Regression formula: log_weight = 4.77 + 0.40log_canopy_area + -0.37log_canopy_shape + 0.37log_canopy_volume + 0.40log_total_height + 0.11log_density + 0.58group

Interpretation:

1. canopy volume and area are positively associated with leaf weight - makes sense - larger canopy means more leaves, so more weight. However large SEs, likely due to collinearity between canopy area and volume, as they are both derived from diam1, diam2, and canopy_height.
2. canopy shape is negatively associated with leaf weight - more spherical canopies have more leaves. Again large SE, likely due to collinearity.
3. Total height is positively associated with leaf weight, which makes sense as taller bushes are likely to have more leaves. However, the SE is large, again, collinearity.
4. Not clear how to interpret density, as no prior information to suggest important. Inclined to exclude.
5. Group coefficient is large and SE small.

# Re-fit both models (no plots) to get traces for comparison

trace_log5 = fit_and_plot_bayes(
    mesquite, ['log_canopy_shape', 'log_canopy_volume', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)

trace_log3 = fit_and_plot_bayes(
    mesquite, ['log_canopy_area', 'log_canopy_shape', 'log_canopy_volume', 'log_total_height', 'log_density', 'group'], 'log_weight',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=False, show_posterior=False, show_regression=False, show_residuals=False
)


# Pointwise LOO for both models
loo_log5 = az.loo(trace_log5, pointwise=True)
loo_log3 = az.loo(trace_log3, pointwise=True)

print(f"LOO log5 model: {loo_log5.elpd_loo:.1f} ± {loo_log5.se:.1f}")
print(f"LOO log3 model: {loo_log3.elpd_loo:.1f} ± {loo_log3.se:.1f}")

az.compare({"log5": loo_log5, "log3": loo_log3})

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_shape, slope_log_canopy_volume, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 1 seconds.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.882  0.139     4.611      5.161      0.002   
slope_log_canopy_shape  -0.380  0.222    -0.812      0.062      0.003   
slope_log_canopy_volume  0.806  0.054     0.697      0.911      0.001   
slope_group              0.585  0.122     0.339      0.817      0.002   
sigma                    0.341  0.039     0.267      0.416      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    3286.0    4017.0    1.0  
slope_log_canopy_shape     0.003    4667.0    4790.0    1.0  
slope_log_canopy_volume    0.001    3734.0    4773.0    1.0  
slope_group                0.001    4087.0    5030.0    1.0  
sigma                      0.001    5656.0    4966.0    1.0  

Regression formula: log_weight = 4.88 + -0.38*log_canopy_shape + 0.81*log_canopy_volume + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.870 ± 0.016
LOO-ELPD: -18.06 ± 5.29  (p_loo=5.4) <- should be lower than number of parameters (5) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.393 ± 0.115

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_log_canopy_area, slope_log_canopy_shape, slope_log_canopy_volume, slope_log_total_height, slope_log_density, slope_group, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
There were 3 divergences after tuning. Increase `target_accept` or reparameterize.

                          mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  \
intercept                4.768  0.160     4.439      5.072      0.003   
slope_log_canopy_area    0.393  0.305    -0.220      0.971      0.006   
slope_log_canopy_shape  -0.376  0.235    -0.823      0.102      0.004   
slope_log_canopy_volume  0.380  0.289    -0.164      0.959      0.005   
slope_log_total_height   0.391  0.321    -0.204      1.051      0.005   
slope_log_density        0.107  0.127    -0.137      0.361      0.002   
slope_group              0.584  0.132     0.339      0.852      0.002   
sigma                    0.340  0.040     0.265      0.421      0.001   

                         mcse_sd  ess_bulk  ess_tail  r_hat  
intercept                  0.002    3940.0    5013.0    1.0  
slope_log_canopy_area      0.004    2925.0    4045.0    1.0  
slope_log_canopy_shape     0.003    4297.0    5047.0    1.0  
slope_log_canopy_volume    0.004    2885.0    4012.0    1.0  
slope_log_total_height     0.004    3798.0    4684.0    1.0  
slope_log_density          0.002    5958.0    4782.0    1.0  
slope_group                0.002    5135.0    5219.0    1.0  
sigma                      0.001    4921.0    3657.0    1.0  

Regression formula: log_weight = 4.77 + 0.39*log_canopy_area + -0.38*log_canopy_shape + 0.38*log_canopy_volume + 0.39*log_total_height + 0.11*log_density + 0.58*group
Residual std dev (σ): 0.34 ± 0.04
Bayesian R²: 0.872 ± 0.016
LOO-ELPD: -19.29 ± 5.26  (p_loo=7.5) <- should be lower than number of parameters (8) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -0.419 ± 0.114
LOO log5 model: -18.1 ± 5.3
LOO log3 model: -19.3 ± 5.3

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
log5	0	-18.064652	5.399723	0.00000	1.000000e+00	5.294748	0.00000	False	log
log3	1	-19.290272	7.464691	1.22562	5.551115e-17	5.260358	1.46315	False	log

Model is better and has less predictors. Dropped collinear predictors, marginal uncertainty reduced. Leaf weight is almost proportional to canopy volume at 0.8 vs 0.4 in previous model.

A simpler model with weak priors means data drives the inference. When we have highly correlated predictors, the coefficients are poorly identified and estimates are noisy. Excluding or combining can get more stable estimates and easier to interpret. Alternatively keep all predictors but use stronger priors.

Finally interactions would be good, but there is not enough data to estimate them reliably.

To conclude: it was better to transform predictors and outcome to log scale, and to drop some predictors to get more stable estimates and better predictive performance. Cross validation is useful for comparing models. Dont chase too many models, or it could result in overfitting.

12.7 Models for regression coefficients

Demonstrate standardisation of predictors and models for regression coefficients.

Example of high school studens in Portugal.

Predict final year maths grade for 343 stidents using a large number of predictors. First not standardised, then standardised.

# import student-merged-all.csv

students = pd.read_csv('../ros_data/student-merged-all.csv', skiprows=0)
display(students.head())


predictors = [
    "school", "sex", "age", "address", "famsize", "Pstatus", "Medu", "Fedu",
    "traveltime", "studytime", "failures", "schoolsup", "famsup", "paid", "activities",
    "nursery", "higher", "internet", "romantic", "famrel", "freetime", "goout", "Dalc",
    "Walc", "health", "absences"
]

data_G3mat = students.loc[students["G3mat"] > 0, ["G3mat"] + predictors]

student_fit_1 = fit_and_plot_bayes(
    data_G3mat, predictors, 'G3mat',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=True, show_posterior=False, show_regression=False, show_residuals=False
)

	G1mat	G2mat	G3mat	G1por	G2por	G3por	school	sex	age	address	...	higher	internet	romantic	famrel	freetime	goout	Dalc	Walc	health	absences
0	7.0	10.0	10.0	13.0	13.0	13.0	0	0	15	0	...	1.0	1	0.0	3.0	1.0	2.0	1.0	1.0	1.0	2.0
1	NaN	NaN	NaN	8.0	9.0	9.0	0	0	15	0	...	NaN	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	8.0	6.0	5.0	13.0	11.0	11.0	0	0	15	0	...	1.0	1	1.0	3.0	3.0	4.0	2.0	4.0	5.0	2.0
3	14.0	13.0	13.0	14.0	13.0	12.0	0	0	15	0	...	1.0	0	0.0	4.0	3.0	1.0	1.0	1.0	2.0	8.0
4	10.0	9.0	8.0	10.0	11.0	10.0	0	0	15	0	...	1.0	1	0.0	4.0	3.0	2.0	1.0	1.0	5.0	2.0

5 rows × 32 columns

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_school, slope_sex, slope_age, slope_address, slope_famsize, slope_Pstatus, slope_Medu, slope_Fedu, slope_traveltime, slope_studytime, slope_failures, slope_schoolsup, slope_famsup, slope_paid, slope_activities, slope_nursery, slope_higher, slope_internet, slope_romantic, slope_famrel, slope_freetime, slope_goout, slope_Dalc, slope_Walc, slope_health, slope_absences, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 23 seconds.

                    mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept         13.667  2.906     8.097     19.459      0.050    0.033   
slope_school      -0.683  0.555    -1.756      0.419      0.007    0.006   
slope_sex          0.561  0.357    -0.135      1.248      0.004    0.004   
slope_age         -0.108  0.147    -0.405      0.174      0.002    0.002   
slope_address      0.515  0.408    -0.285      1.304      0.005    0.004   
slope_famsize      0.315  0.341    -0.339      0.991      0.004    0.004   
slope_Pstatus     -0.395  0.505    -1.401      0.574      0.005    0.005   
slope_Medu         0.246  0.185    -0.110      0.610      0.002    0.002   
slope_Fedu         0.234  0.182    -0.114      0.603      0.002    0.002   
slope_traveltime   0.044  0.246    -0.430      0.534      0.003    0.002   
slope_studytime    0.438  0.204     0.044      0.835      0.002    0.002   
slope_failures    -0.741  0.253    -1.222     -0.232      0.003    0.003   
slope_schoolsup   -2.273  0.463    -3.145     -1.328      0.005    0.005   
slope_famsup      -0.452  0.337    -1.135      0.189      0.004    0.004   
slope_paid        -0.556  0.330    -1.212      0.074      0.004    0.004   
slope_activities   0.116  0.318    -0.488      0.747      0.004    0.003   
slope_nursery     -0.321  0.403    -1.114      0.472      0.004    0.004   
slope_higher       0.233  0.835    -1.423      1.877      0.010    0.009   
slope_internet     0.650  0.441    -0.175      1.531      0.005    0.005   
slope_romantic    -0.114  0.344    -0.791      0.541      0.004    0.004   
slope_famrel       0.051  0.176    -0.306      0.385      0.002    0.002   
slope_freetime     0.067  0.163    -0.260      0.378      0.002    0.002   
slope_goout       -0.397  0.164    -0.725     -0.089      0.002    0.002   
slope_Dalc         0.067  0.223    -0.354      0.514      0.003    0.002   
slope_Walc        -0.252  0.175    -0.587      0.097      0.002    0.002   
slope_health      -0.183  0.114    -0.400      0.044      0.001    0.001   
slope_absences    -0.068  0.020    -0.108     -0.030      0.000    0.000   
sigma              2.848  0.112     2.643      3.071      0.001    0.001   

                  ess_bulk  ess_tail  r_hat  
intercept           3436.0    4739.0    1.0  
slope_school        6612.0    5395.0    1.0  
slope_sex           7706.0    6315.0    1.0  
slope_age           3646.0    4529.0    1.0  
slope_address       7066.0    6399.0    1.0  
slope_famsize       8993.0    6498.0    1.0  
slope_Pstatus       8946.0    6365.0    1.0  
slope_Medu          5782.0    6077.0    1.0  
slope_Fedu          5912.0    5786.0    1.0  
slope_traveltime    7243.0    6248.0    1.0  
slope_studytime     7853.0    5978.0    1.0  
slope_failures      6467.0    5793.0    1.0  
slope_schoolsup     7244.0    6513.0    1.0  
slope_famsup        7706.0    6329.0    1.0  
slope_paid          8069.0    6493.0    1.0  
slope_activities    8204.0    6301.0    1.0  
slope_nursery       8370.0    6047.0    1.0  
slope_higher        6411.0    5157.0    1.0  
slope_internet      8201.0    6368.0    1.0  
slope_romantic      7054.0    5887.0    1.0  
slope_famrel        8279.0    6308.0    1.0  
slope_freetime      7729.0    6273.0    1.0  
slope_goout         7816.0    6070.0    1.0  
slope_Dalc          5452.0    6106.0    1.0  
slope_Walc          5177.0    5687.0    1.0  
slope_health        7125.0    5630.0    1.0  
slope_absences      6972.0    5897.0    1.0  
sigma               8989.0    6239.0    1.0  

Regression formula: G3mat = 13.67 + -0.68*school + 0.56*sex + -0.11*age + 0.52*address + 0.32*famsize + -0.40*Pstatus + 0.25*Medu + 0.23*Fedu + 0.04*traveltime + 0.44*studytime + -0.74*failures + -2.27*schoolsup + -0.45*famsup + -0.56*paid + 0.12*activities + -0.32*nursery + 0.23*higher + 0.65*internet + -0.11*romantic + 0.05*famrel + 0.07*freetime + -0.40*goout + 0.07*Dalc + -0.25*Walc + -0.18*health + -0.07*absences
Residual std dev (σ): 2.85 ± 0.11
Bayesian R²: 0.302 ± 0.031
LOO-ELPD: -913.11 ± 12.91  (p_loo=25.8) <- should be lower than number of parameters (28) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -2.502 ± 0.035

Effect of absences looks small and certain - below we standardise predictors.

from sklearn.preprocessing import scale

datasd_G3mat = data_G3mat.copy()

datasd_G3mat[predictors] = scale(datasd_G3mat[predictors])

student_fit_2 = fit_and_plot_bayes(
    datasd_G3mat, predictors, 'G3mat',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=50, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=True, show_posterior=False, show_regression=False, show_residuals=False
)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_school, slope_sex, slope_age, slope_address, slope_famsize, slope_Pstatus, slope_Medu, slope_Fedu, slope_traveltime, slope_studytime, slope_failures, slope_schoolsup, slope_famsup, slope_paid, slope_activities, slope_nursery, slope_higher, slope_internet, slope_romantic, slope_famrel, slope_freetime, slope_goout, slope_Dalc, slope_Walc, slope_health, slope_absences, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.

                    mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept         11.534  0.149    11.240     11.824      0.001    0.002   
slope_school      -0.214  0.179    -0.553      0.154      0.002    0.002   
slope_sex          0.277  0.172    -0.057      0.615      0.002    0.002   
slope_age         -0.142  0.183    -0.501      0.210      0.002    0.002   
slope_address      0.210  0.174    -0.136      0.551      0.002    0.002   
slope_famsize      0.145  0.157    -0.149      0.468      0.001    0.002   
slope_Pstatus     -0.123  0.158    -0.431      0.187      0.001    0.002   
slope_Medu         0.273  0.201    -0.142      0.651      0.002    0.002   
slope_Fedu         0.252  0.203    -0.127      0.665      0.002    0.002   
slope_traveltime   0.029  0.167    -0.299      0.351      0.002    0.002   
slope_studytime    0.364  0.173     0.030      0.709      0.002    0.002   
slope_failures    -0.499  0.176    -0.855     -0.165      0.002    0.002   
slope_schoolsup   -0.795  0.163    -1.109     -0.474      0.002    0.002   
slope_famsup      -0.223  0.169    -0.568      0.097      0.002    0.002   
slope_paid        -0.275  0.161    -0.574      0.054      0.002    0.002   
slope_activities   0.058  0.161    -0.252      0.375      0.001    0.002   
slope_nursery     -0.128  0.159    -0.447      0.176      0.001    0.002   
slope_higher       0.043  0.168    -0.292      0.368      0.002    0.002   
slope_internet     0.239  0.161    -0.072      0.559      0.002    0.002   
slope_romantic    -0.053  0.157    -0.361      0.251      0.001    0.002   
slope_famrel       0.047  0.158    -0.266      0.350      0.002    0.002   
slope_freetime     0.069  0.170    -0.253      0.412      0.002    0.002   
slope_goout       -0.432  0.186    -0.792     -0.071      0.002    0.002   
slope_Dalc         0.057  0.202    -0.334      0.471      0.002    0.002   
slope_Walc        -0.324  0.226    -0.757      0.124      0.003    0.002   
slope_health      -0.253  0.155    -0.565      0.043      0.001    0.002   
slope_absences    -0.555  0.168    -0.892     -0.239      0.002    0.002   
sigma              2.850  0.108     2.634      3.058      0.001    0.001   

                  ess_bulk  ess_tail  r_hat  
intercept          13416.0    6018.0    1.0  
slope_school        8483.0    6209.0    1.0  
slope_sex          10199.0    6838.0    1.0  
slope_age           8455.0    7152.0    1.0  
slope_address      10536.0    6433.0    1.0  
slope_famsize      11021.0    6526.0    1.0  
slope_Pstatus      11199.0    6753.0    1.0  
slope_Medu          8389.0    6723.0    1.0  
slope_Fedu          8457.0    6500.0    1.0  
slope_traveltime   11080.0    5811.0    1.0  
slope_studytime     9953.0    6077.0    1.0  
slope_failures     10594.0    6587.0    1.0  
slope_schoolsup    10613.0    6122.0    1.0  
slope_famsup       10964.0    6391.0    1.0  
slope_paid         11181.0    6994.0    1.0  
slope_activities   12128.0    6331.0    1.0  
slope_nursery      11298.0    5672.0    1.0  
slope_higher        9754.0    6120.0    1.0  
slope_internet     10578.0    6630.0    1.0  
slope_romantic     12062.0    6188.0    1.0  
slope_famrel       10819.0    5910.0    1.0  
slope_freetime      9502.0    6430.0    1.0  
slope_goout         8855.0    6662.0    1.0  
slope_Dalc          8190.0    5763.0    1.0  
slope_Walc          6309.0    6355.0    1.0  
slope_health       13244.0    6257.0    1.0  
slope_absences     10461.0    6644.0    1.0  
sigma              10476.0    6086.0    1.0  

Regression formula: G3mat = 11.53 + -0.21*school + 0.28*sex + -0.14*age + 0.21*address + 0.15*famsize + -0.12*Pstatus + 0.27*Medu + 0.25*Fedu + 0.03*traveltime + 0.36*studytime + -0.50*failures + -0.79*schoolsup + -0.22*famsup + -0.27*paid + 0.06*activities + -0.13*nursery + 0.04*higher + 0.24*internet + -0.05*romantic + 0.05*famrel + 0.07*freetime + -0.43*goout + 0.06*Dalc + -0.32*Walc + -0.25*health + -0.55*absences
Residual std dev (σ): 2.85 ± 0.11
Bayesian R²: 0.302 ± 0.032
LOO-ELPD: -913.39 ± 12.90  (p_loo=26.0) <- should be lower than number of parameters (28) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -2.502 ± 0.035
LOO R²: 0.166
  *** WARNING: LOO R² (0.166) < Bayesian R² (0.302) — model may be overfitting ***

After standardising predictors, the coefficient for absences is larger and less certain, which suggests that the effect of absences on final year maths grade is more substantial than it appeared in the unstandardised model. This is because standardising the predictors allows us to compare the coefficients on a common scale, which can reveal relationships that may not be apparent in the unstandardised model.

LOO R² (0.168) < Bayesian R² (0.302) — model may be overfitting

We have 26 predictors each has SD of 1. We use gentle prior: normal mean 0 and SD 2.5 - modest for single predictor.

The issue: Model is predicting linear combination, all 26 predictors and their coefficients add together. Because priors are independent, variances add up. Variance of predicted means is 26 * 2.5^2 = 162.5, and SD is 12.75. Meanwhile prior for noise \sigma SD is 3.3.

Setup: R2 is explained variance / total variance = var(predicted means) / (var(predicted means) + var(noise)^2). This is saying signal has SD 12.7 and noise 3.3 - R2 is being pushed towards 1.

Lesson: weakly informative priors per coefficient can be very informative for things we care about (R2) when we have many predictors.

The fix: decide on a budget. For example, guess that predictors together explain 30% of variance:

Explained part = 0.3 * (var(y), split across 26 predictors, so \sqrt(0.3 / 26) * sd(y) Residual part = 0.7 * (var(y) - \sigma gets prior mean \sqrt(0.7) * sd(y)

Now total signal variance plus noise variance = actual data variance.

prior_scale = datasd_G3mat['G3mat'].std() * np.sqrt(0.3 / 26)

student_fit_2 = fit_and_plot_bayes(
    datasd_G3mat, predictors, 'G3mat',
    intercept_mu=0, intercept_sigma=50, slope_mu=0, slope_sigma=prior_scale, sigma_sigma=50,
    samples=2000, tune=1000,
    show_trace=False, show_forest=True, show_posterior=False, show_regression=False, show_residuals=False
)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slope_school, slope_sex, slope_age, slope_address, slope_famsize, slope_Pstatus, slope_Medu, slope_Fedu, slope_traveltime, slope_studytime, slope_failures, slope_schoolsup, slope_famsup, slope_paid, slope_activities, slope_nursery, slope_higher, slope_internet, slope_romantic, slope_famrel, slope_freetime, slope_goout, slope_Dalc, slope_Walc, slope_health, slope_absences, sigma]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.

                    mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  \
intercept         11.535  0.148    11.248     11.820      0.001    0.002   
slope_school      -0.159  0.154    -0.457      0.136      0.001    0.002   
slope_sex          0.259  0.151    -0.026      0.559      0.001    0.002   
slope_age         -0.143  0.157    -0.451      0.157      0.001    0.002   
slope_address      0.184  0.147    -0.111      0.466      0.001    0.002   
slope_famsize      0.123  0.143    -0.146      0.411      0.001    0.002   
slope_Pstatus     -0.100  0.142    -0.384      0.176      0.001    0.002   
slope_Medu         0.256  0.166    -0.063      0.586      0.002    0.002   
slope_Fedu         0.210  0.166    -0.092      0.553      0.002    0.002   
slope_traveltime  -0.007  0.145    -0.306      0.262      0.001    0.002   
slope_studytime    0.299  0.150     0.007      0.594      0.001    0.002   
slope_failures    -0.454  0.151    -0.740     -0.142      0.001    0.002   
slope_schoolsup   -0.665  0.144    -0.939     -0.370      0.001    0.002   
slope_famsup      -0.189  0.144    -0.480      0.086      0.001    0.002   
slope_paid        -0.218  0.147    -0.523      0.051      0.001    0.002   
slope_activities   0.057  0.140    -0.213      0.332      0.001    0.002   
slope_nursery     -0.095  0.143    -0.378      0.176      0.001    0.002   
slope_higher       0.059  0.150    -0.240      0.350      0.001    0.002   
slope_internet     0.207  0.143    -0.085      0.479      0.001    0.002   
slope_romantic    -0.051  0.142    -0.334      0.220      0.001    0.002   
slope_famrel       0.046  0.142    -0.222      0.337      0.001    0.002   
slope_freetime     0.043  0.150    -0.246      0.337      0.001    0.002   
slope_goout       -0.364  0.153    -0.657     -0.055      0.001    0.002   
slope_Dalc        -0.013  0.168    -0.349      0.314      0.002    0.002   
slope_Walc        -0.275  0.177    -0.623      0.072      0.002    0.002   
slope_health      -0.207  0.143    -0.483      0.068      0.001    0.002   
slope_absences    -0.470  0.148    -0.762     -0.184      0.001    0.002   
sigma              2.838  0.108     2.624      3.051      0.001    0.001   

                  ess_bulk  ess_tail  r_hat  
intercept          18064.0    5667.0    1.0  
slope_school       12386.0    6078.0    1.0  
slope_sex          12162.0    6410.0    1.0  
slope_age          13259.0    6762.0    1.0  
slope_address      13427.0    6171.0    1.0  
slope_famsize      16883.0    5566.0    1.0  
slope_Pstatus      14885.0    6388.0    1.0  
slope_Medu         12093.0    6495.0    1.0  
slope_Fedu         10966.0    5892.0    1.0  
slope_traveltime   14582.0    6586.0    1.0  
slope_studytime    12889.0    6312.0    1.0  
slope_failures     14834.0    6380.0    1.0  
slope_schoolsup    15160.0    6783.0    1.0  
slope_famsup       12776.0    6826.0    1.0  
slope_paid         14289.0    6341.0    1.0  
slope_activities   15789.0    6661.0    1.0  
slope_nursery      14115.0    6267.0    1.0  
slope_higher       14519.0    5911.0    1.0  
slope_internet     12767.0    5826.0    1.0  
slope_romantic     15212.0    6059.0    1.0  
slope_famrel       16214.0    5835.0    1.0  
slope_freetime     13235.0    6442.0    1.0  
slope_goout        12318.0    6371.0    1.0  
slope_Dalc         11314.0    6398.0    1.0  
slope_Walc          8993.0    6210.0    1.0  
slope_health       15532.0    5819.0    1.0  
slope_absences     12646.0    6688.0    1.0  
sigma              14365.0    5837.0    1.0  

Regression formula: G3mat = 11.54 + -0.16*school + 0.26*sex + -0.14*age + 0.18*address + 0.12*famsize + -0.10*Pstatus + 0.26*Medu + 0.21*Fedu + -0.01*traveltime + 0.30*studytime + -0.45*failures + -0.67*schoolsup + -0.19*famsup + -0.22*paid + 0.06*activities + -0.10*nursery + 0.06*higher + 0.21*internet + -0.05*romantic + 0.05*famrel + 0.04*freetime + -0.36*goout + -0.01*Dalc + -0.27*Walc + -0.21*health + -0.47*absences
Residual std dev (σ): 2.84 ± 0.11
Bayesian R²: 0.254 ± 0.031
LOO-ELPD: -908.96 ± 12.70  (p_loo=21.0) <- should be lower than number of parameters (28) for reliable estimates
When comparing models, better performing models have higher (less negative) LOO-ELPD - a meaningful difference between models is 2x SE.
LOO log score (per obs): -2.490 ± 0.035
LOO R²: 0.186
  *** WARNING: LOO R² (0.186) < Bayesian R² (0.254) — model may be overfitting ***

However, if we think that only a few predictors are important, we can use regularised horseshoe priors, which are designed to handle situations where we have many predictors but only a few of them are expected to have a substantial effect on the outcome variable. The regularised horseshoe prior allows for some coefficients to be large while shrinking the others towards zero.

Instead of dividing the budget across 26 predictors, we assume say 6 are important and 20 not. But we dont know which 6 are important, so we use a regularised horseshoe prior, which allows for some coefficients to be large while shrinking the others towards zero.

How: Each coefficient gets a normal prior with mean 0 and SD τ·λ_j

τ global scale - how much all get shrunk towards zero λ_j local scale - how much each coefficient gets shrunk towards zero

λ has a clever prior: a half-Cauchy, which is mostly small but big tail. Most get small λ, but some can draw a large λ, which allows for some coefficients to be large while shrinking the others towards zero. Predictors are encouraged to be large or small, but not in between.

Why horseshoe? Because of the either large or small nature of the coefficients, the prior distribution has a shape that resembles a horseshoe when plotted - values bunched near 1 and 0. Prior is essentially saying predictor is either important (large coefficient) or not important (small coefficient), but not in between.

Spike and slab - spike is the sharp peak at zero. Slab is is wide flat distribution, to allow coefficients to roam. Plain horseshoe can allow escaping coefficients to be very large. The regularised horseshoe adds a regularisation term to the slab, which prevents coefficients from becoming too large, while still allowing for some coefficients to be large if the data supports it.

If we expect 6 out of 26 (p0 = 6). Slab scale is same as above, but for 6 instead of 26.

The payoff: prior assumes that most predictors are not important, so variance does not pile up like before. R2 is more cautious. Nudges towards simpler models.

Summary: even budget above, says everyone gets a small slice, horseshoe says most get no slice, but some get a big slice and that is decided by the data.

# --- Horseshoe prior parameters ---
# p0: our prior belief about how many predictors actually matter
# global_scale: shrinks most coefficients toward zero globally
# slab_scale: allows the ~p0 important coefficients to escape shrinkage

p  = len(predictors)
n  = len(datasd_G3mat)
p0 = 6  # expected number of relevant predictors

slab_scale   = np.sqrt(0.3 / p0) * datasd_G3mat['G3mat'].std()
global_scale = (p0 / (p - p0)) / np.sqrt(n)
slab_df      = 4  # rstanarm default; controls how heavy the slab tails are

y = datasd_G3mat['G3mat'].values
X = datasd_G3mat[predictors].values

with pm.Model() as hs_model:
    sigma = pm.HalfNormal("sigma", sigma=50)

    # Global shrinkage (scaled by sigma to match rstanarm's internal scaling)
    tau = pm.HalfCauchy("tau", beta=global_scale * sigma)

    # Local shrinkage: each predictor gets its own scale, allowing a few to be large
    lam = pm.HalfCauchy("lam", beta=1, shape=p)

    # Slab: prevents important coefficients from being over-shrunk
    c2 = pm.InverseGamma("c2", alpha=slab_df / 2, beta=slab_df / 2 * slab_scale**2)

    # Regularised local scale combining global, local, and slab (Piironen & Vehtari 2017)
    lam_tilde = lam * pm.math.sqrt(c2 / (c2 + tau**2 * lam**2))

    intercept = pm.Normal("intercept", mu=0, sigma=50)
    beta      = pm.Normal("beta", mu=0, sigma=tau * lam_tilde, shape=p)

    mu = intercept + pm.math.dot(X, beta)
    pm.Normal("y", mu=mu, sigma=sigma, observed=y)

    student_fit_3 = pm.sample(2000, tune=1000, idata_kwargs={"log_likelihood": True})

# Print posterior summary for intercept and slopes only (excluding sampler internals)
print(pm.summary(student_fit_3, var_names=["intercept", "beta"]))

# Label beta coefficients by predictor name for readability
coord_map = {f"beta[{i}]": name for i, name in enumerate(predictors)}
az.plot_forest(
    student_fit_3,
    var_names=["intercept", "beta"],
    combined=True,
    figsize=(8, 10)
)
# Replace generic "beta[i]" y-tick labels with predictor names
ax = plt.gca()
ax.set_yticklabels([coord_map.get(t.get_text(), t.get_text()) for t in ax.get_yticklabels()])
plt.title("Horseshoe prior: posterior coefficients")
plt.tight_layout()
The plt.show()

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, tau, lam, c2, intercept, beta]

/opt/anaconda3/envs/ros_pymc/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" for 
Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 9 seconds.
There were 1578 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

             mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \
intercept  11.539  0.148  11.259   11.813      0.004    0.002    1605.0   
beta[0]    -0.119  0.139  -0.408    0.107      0.004    0.004    1409.0   
beta[1]     0.174  0.161  -0.061    0.483      0.007    0.002     425.0   
beta[2]    -0.109  0.139  -0.405    0.103      0.003    0.003    1906.0   
beta[3]     0.145  0.150  -0.075    0.433      0.008    0.003     456.0   
beta[4]     0.061  0.114  -0.149    0.292      0.004    0.004     804.0   
beta[5]    -0.067  0.116  -0.352    0.100      0.005    0.005     611.0   
beta[6]     0.233  0.184  -0.064    0.566      0.007    0.002     975.0   
beta[7]     0.141  0.161  -0.088    0.477      0.005    0.004     863.0   
beta[8]    -0.027  0.102  -0.242    0.171      0.002    0.002    1602.0   
beta[9]     0.190  0.176  -0.074    0.515      0.009    0.005     363.0   
beta[10]   -0.503  0.177  -0.827   -0.173      0.007    0.003     675.0   
beta[11]   -0.697  0.156  -0.981   -0.382      0.003    0.002    2883.0   
beta[12]   -0.111  0.134  -0.374    0.096      0.004    0.002    1364.0   
beta[13]   -0.149  0.146  -0.438    0.076      0.003    0.002    1892.0   
beta[14]    0.042  0.104  -0.167    0.250      0.003    0.003    1472.0   
beta[15]   -0.022  0.108  -0.250    0.180      0.007    0.003     271.0   
beta[16]    0.030  0.102  -0.151    0.253      0.002    0.002    1488.0   
beta[17]    0.142  0.145  -0.076    0.420      0.005    0.002    1088.0   
beta[18]   -0.029  0.101  -0.250    0.150      0.002    0.002    1336.0   
beta[19]    0.033  0.102  -0.164    0.234      0.005    0.002     416.0   
beta[20]    0.012  0.105  -0.189    0.230      0.003    0.002     865.0   
beta[21]   -0.315  0.182  -0.626    0.022      0.005    0.002    1315.0   
beta[22]   -0.022  0.115  -0.272    0.205      0.003    0.003    1816.0   
beta[23]   -0.225  0.189  -0.594    0.060      0.006    0.002    1232.0   
beta[24]   -0.122  0.132  -0.390    0.077      0.003    0.002    1924.0   
beta[25]   -0.482  0.167  -0.801   -0.173      0.003    0.003    2308.0   

           ess_tail  r_hat  
intercept    4705.0   1.01  
beta[0]       592.0   1.01  
beta[1]       670.0   1.01  
beta[2]      3014.0   1.00  
beta[3]      1146.0   1.01  
beta[4]       630.0   1.01  
beta[5]       160.0   1.01  
beta[6]      3333.0   1.00  
beta[7]       635.0   1.01  
beta[8]      2198.0   1.00  
beta[9]       649.0   1.01  
beta[10]     1186.0   1.01  
beta[11]     2314.0   1.01  
beta[12]     2180.0   1.00  
beta[13]     3655.0   1.02  
beta[14]     2109.0   1.01  
beta[15]      245.0   1.02  
beta[16]     3283.0   1.01  
beta[17]     1664.0   1.01  
beta[18]     3132.0   1.01  
beta[19]     3146.0   1.01  
beta[20]      948.0   1.00  
beta[21]     2581.0   1.00  
beta[22]     1727.0   1.01  
beta[23]     3300.0   1.00  
beta[24]     1611.0   1.01  
beta[25]     3135.0   1.01  

The forest plot shows most coefficients are small and close to zero, but a few are larger.

Using this, we have a better idea of which predictors are important, and we can build a model with just those predictors.

However, we dont know how many predictors are important.

Other models for regression coefficients

We can use different regularisation priors, not just horseshoe. E.g., ridge, lasso, horseshoe etc.

We can work back from R2. For some prior choices, we can write a formula for implied prior on R2. I.e. I want R2 prior to look like this and then pick the prior on the coefficients that gives me that R2 prior.

Horseshoe is better than lasso when we have many predictors but only a few are expected to be important, as it allows for some coefficients to be large while shrinking the others towards zero, whereas lasso shrinks all coefficients towards zero, which can lead to biased estimates for the important predictors.

Using to choose predictors and warning. We can fit horseshoe and then build a model with fewer predictors. But only if predictors are perfectly independent. Two correlated predictors can be wrongly discarded, when the pair is important but the individual predictors are not.

Big takeaway: gentle priors can be strong on things like R2, fix is to think about sensible budget and structure of the problem (every parameter matters a small amount or few parameters matter a lot).

12.8 Bibliographic note

12.9 Exercises