Research & Monitoring / Cash Flow & Risk Modelling / Methodology

Bayesian Hierarchical Modelling with PyMC

Bayesian hierarchical modelling, also called multilevel modelling, is a statistical framework for complex ecological data. It allows TT-OS to model variation at several levels at once: individual plants, plots, farm blocks and the wider site.

At the Elegguá Project, this approach can support soil organic carbon monitoring, carbon sequestration estimates, plant performance analysis, biodiversity indicators and future additionality work. The purpose is not to pretend the project has perfect data. The purpose is to model uncertainty honestly and improve estimates as new field evidence arrives.

Why Hierarchical Bayesian Models?

Ecological systems are noisy. One plot may show high soil carbon gains while another plot nearby may show slower change because of shade, species mix, soil texture, slope, water availability or management history. A simple average can hide this structure. A hierarchical model keeps the structure visible.

• Partial pooling: individual plot estimates are pulled towards the group mean when data is limited. This reduces overfitting on noisy small samples while still respecting local differences.
• Prior knowledge: regional West African agroforestry studies and early TT-OS observations can inform prior assumptions.
• Uncertainty propagation: the model naturally produces distributions that can feed into carbon, NPV, risk and confidence-score calculations.
• Dynamic updating: posterior results from one monitoring season can become priors for the next season.
• Spatial hierarchy: plot-level data can be understood within block-level and farm-level patterns.

Applications at the Elegguá Project

PyMC Code Example: Hierarchical Model for Soil Carbon Sequestration

The example below models observed annual soil organic carbon change across forty plots grouped into five farm blocks. It is a practical starting point and should be adapted with real soil testing data, field identifiers and verified measurement units.

import pymc as pm
import numpy as np
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt

# Example data structure.
# Replace this simulated data with real TT-OS soil testing data.
np.random.seed(42)

n_plots = 40

data = pd.DataFrame({
    "plot_id": np.arange(n_plots),
    "block_id": np.repeat(np.arange(5), 8),
    "companion_species": np.random.choice(
        ["Plantain", "Breadfruit", "Mixed"],
        n_plots
    ),
    "soil_depth_cm": np.random.normal(30, 5, n_plots),
    "initial_soc": np.random.normal(25, 8, n_plots),
    "years_since_planting": np.random.uniform(1, 10, n_plots),
    "observed_soc_change": np.random.normal(2.5, 1.5, n_plots)
})

plot_idx = data["plot_id"].values
block_idx = data["block_id"].values
n_blocks = data["block_id"].nunique()

with pm.Model() as soil_carbon_model:
    # Farm-level hyperpriors
    mu_global = pm.Normal("mu_global", mu=2.0, sigma=1.5)
    sigma_global = pm.HalfNormal("sigma_global", sigma=2.0)

    # Block-level variation
    mu_block = pm.Normal(
        "mu_block",
        mu=mu_global,
        sigma=sigma_global,
        shape=n_blocks
    )

    sigma_block = pm.HalfNormal("sigma_block", sigma=1.0)

    # Plot-level parameters with partial pooling
    mu_plot = pm.Normal(
        "mu_plot",
        mu=mu_block[block_idx],
        sigma=sigma_block,
        shape=n_plots
    )

    # Observation noise
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=2.0)

    observed_change = pm.Normal(
        "observed_change",
        mu=mu_plot[plot_idx],
        sigma=sigma_obs,
        observed=data["observed_soc_change"]
    )

    trace = pm.sample(
        2000,
        tune=1000,
        target_accept=0.9,
        return_inferencedata=True
    )

print(az.summary(trace, var_names=["mu_global", "mu_block"]))

az.plot_trace(trace, var_names=["mu_global"])
plt.show()

with soil_carbon_model:
    posterior_predictive = pm.sample_posterior_predictive(trace)

Key Features of This Model

• Three levels: global farm-level estimate, block-level variation and plot-level estimates.
• Partial pooling: plot estimates are stabilised by information from the wider block and farm.
• Extensible covariates: companion species, soil depth, texture, shade, slope, water condition and years since planting can be added.
• Diagnostics: ArviZ can be used for trace plots, posterior summaries and model comparison.
• Posterior predictive checks: the model can be checked against observed patterns rather than accepted blindly.

Extension: Add Companion Species as a Covariate

A more advanced model can estimate whether different companion species groups are associated with different soil carbon change patterns. This should be interpreted carefully because species mix may be confounded with plot age, management, shade, soil texture or water condition.

# Example extension concept.
# Convert companion_species to an index before modelling.

species_codes, species_index = np.unique(
    data["companion_species"],
    return_inverse=True
)

with pm.Model() as soil_carbon_species_model:
    mu_global = pm.Normal("mu_global", mu=2.0, sigma=1.5)
    sigma_global = pm.HalfNormal("sigma_global", sigma=2.0)

    mu_block = pm.Normal(
        "mu_block",
        mu=mu_global,
        sigma=sigma_global,
        shape=n_blocks
    )

    beta_species = pm.Normal(
        "beta_species",
        mu=0.0,
        sigma=0.8,
        shape=len(species_codes)
    )

    beta_years = pm.Normal("beta_years", mu=0.0, sigma=0.3)

    expected_change = (
        mu_block[block_idx]
        + beta_species[species_index]
        + beta_years * data["years_since_planting"].values
    )

    sigma_obs = pm.HalfNormal("sigma_obs", sigma=2.0)

    observed_change = pm.Normal(
        "observed_change",
        mu=expected_change,
        sigma=sigma_obs,
        observed=data["observed_soc_change"]
    )

    trace_species = pm.sample(
        2000,
        tune=1000,
        target_accept=0.9,
        return_inferencedata=True
    )

Integration with Monte Carlo Cash Flow Modelling

The Bayesian model does not need to stand alone. Posterior distributions can feed directly into the TT-OS cash-flow and risk modelling layer. Instead of using a single fixed sequestration assumption, the Monte Carlo model can sample sequestration values from the Bayesian posterior.

# Conceptual integration.
# posterior_mu_global would be extracted from the PyMC trace.

posterior_mu_global = trace.posterior["mu_global"].values.flatten()

def sample_sequestration_rate(rng):
    return rng.choice(posterior_mu_global)

def carbon_revenue_from_posterior(rng, hectares, carbon_price):
    sequestration_rate = sample_sequestration_rate(rng)
    return sequestration_rate * hectares * carbon_price

This helps TT-OS communicate uncertainty more honestly. A model can report a range, probability and credible interval rather than one over-confident point estimate.

Implementation Workflow

1. Collect structured field data: plot ID, block ID, date, soil depth, SOC, bulk density, companion species, shade context and management notes.
2. Clean and validate the CSV: check units, missing values, duplicate plot IDs and measurement dates.
3. Run a baseline model: start with global, block and plot levels before adding many covariates.
4. Check diagnostics: inspect trace plots, effective sample size, R-hat and posterior predictive checks.
5. Add covariates carefully: companion species, years since planting, soil depth, shade and moisture context.
6. Connect to TT-OS: feed credible intervals and posterior summaries into confidence scoring, monitoring reports and scenario models.
7. Update annually: use new soil testing data to update posteriors and refine assumptions.

Data Columns Recommended for TT-OS Soil Carbon Modelling

Website-Ready Summary

TT-OS can apply PyMC to build hierarchical models of soil carbon dynamics and agroforestry performance. This approach accounts for variation across spatial scales, provides uncertainty quantification, supports adaptive management and strengthens future carbon feasibility work.

The model is open to research collaboration. Partners experienced in probabilistic programming, soil science, agroforestry, carbon monitoring or ecological statistics can help refine the model using real TT-OS field data.

Next Steps

• Create a clean soil testing CSV template.
• Add unit rules for SOC, bulk density and depth.
• Build a notebook-style version for research partners.
• Connect posterior estimates into the Monte Carlo cash-flow model.
• Publish periodic monitoring summaries with credible intervals and clear limitations.

Open Cash Flow & Risk Modelling | Open Confidence Scoring | Open Research & Monitoring