From 4b31a9a641e369e6a8489cf3812c153b27f4e1b1 Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 17:11:15 -0400 Subject: [PATCH 1/3] docs: migrate centered parameterization tutorial (#1273) --- .../centered_vs_noncentered_basic_logic.ipynb | 2083 ++++++++++------- .../centered_vs_noncentered_basic_logic.py | 759 ++++++ 2 files changed, 1941 insertions(+), 901 deletions(-) create mode 100644 docs/tutorials/centered_vs_noncentered_basic_logic.py diff --git a/docs/tutorials/centered_vs_noncentered_basic_logic.ipynb b/docs/tutorials/centered_vs_noncentered_basic_logic.ipynb index e33b088ae..be4c1b247 100644 --- a/docs/tutorials/centered_vs_noncentered_basic_logic.ipynb +++ b/docs/tutorials/centered_vs_noncentered_basic_logic.ipynb @@ -1,248 +1,388 @@ { "cells": [ { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, + "id": "Hbol", + "metadata": {}, + "outputs": [], + "source": [ + "import logging\n", + "import os\n", + "import warnings\n", + "from tempfile import gettempdir\n", + "\n", + "# This notebook only constructs models. Molab can expose a CUDA plugin even\n", + "# when no usable GPU is present, so select CPU before importing HSSM/JAX.\n", + "os.environ[\"JAX_PLATFORMS\"] = \"cpu\"\n", + "os.environ[\"JAX_SKIP_CUDA_CONSTRAINTS_CHECK\"] = \"1\"\n", + "os.environ[\"XLA_PYTHON_CLIENT_PREALLOCATE\"] = \"false\"\n", + "os.environ.setdefault(\n", + " \"MPLCONFIGDIR\", f\"{gettempdir()}/hssm-parameterization-matplotlib\"\n", + ")\n", + "\n", + "warnings.filterwarnings(\"ignore\")\n", + "logging.getLogger(\"jax._src.xla_bridge\").setLevel(logging.CRITICAL)\n", + "logging.getLogger(\"matplotlib\").setLevel(logging.ERROR)\n", + "\n", + "import bambi as bmb\n", + "import marimo as mo\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pymc as pm\n", + "from pytensor.graph.traversal import ancestors\n", + "\n", + "import hssm\n", + "\n", + "logging.getLogger(\"hssm\").setLevel(logging.WARNING)\n", + "hssm.set_floatX(\"float64\")\n", + "pd.set_option(\"display.max_colwidth\", 100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "MJUe", "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "# Centered vs. non-centered parameterizations\n", + "\n", + "Centering is a choice about **how a hierarchical effect is represented for\n", + "computation**. It should not change the statistical model. That distinction\n", + "becomes important when a group-specific prior has its own population mean.\n", + "\n", + "This tutorial uses **HSSM 0.4.0** and **Bambi 0.20.0**.\n", + "It constructs models and inspects their PyMC graphs; no MCMC is run.\n", + "\n", + "By the end, you should be able to:\n", + "\n", + "1. distinguish the mathematical non-centered transformation from Bambi's\n", + " current shortcut;\n", + "2. recognize a valid zero-mean group deviation in either parameterization;\n", + "3. understand why HSSM rejects an explicit non-centered free group mean; and\n", + "4. choose exactly one owner for each population location." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ + "mo.md(f\"\"\"\n", "# Centered vs. non-centered parameterizations\n", "\n", - "When you fit a hierarchical model in HSSM you implicitly choose between two parameterizations of the group-specific (random) effects. They are statistically equivalent, but **the node names in the underlying PyMC graph differ**, and that difference interacts with the way you specify priors. Misalignment between the two leads to a subtle footgun: a prior you supplied can be silently dropped, leaving a *disconnected* free RV in the graph.\n", + "Centering is a choice about **how a hierarchical effect is represented for\n", + "computation**. It should not change the statistical model. That distinction\n", + "becomes important when a group-specific prior has its own population mean.\n", + "\n", + "This tutorial uses **HSSM {hssm.__version__}** and **Bambi {bmb.__version__}**.\n", + "It constructs models and inspects their PyMC graphs; no MCMC is run.\n", "\n", - "This tutorial walks through:\n", + "By the end, you should be able to:\n", "\n", - "1. The two parameterizations and the equations behind them.\n", - "2. How HSSM/bambi name the PyMC nodes in each case.\n", - "3. The disconnected-node footgun, the warning HSSM now emits, and two ways to fix it." + "1. distinguish the mathematical non-centered transformation from Bambi's\n", + " current shortcut;\n", + "2. recognize a valid zero-mean group deviation in either parameterization;\n", + "3. understand why HSSM rejects an explicit non-centered free group mean; and\n", + "4. choose exactly one owner for each population location.\n", + "\"\"\")" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "vblA", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "## 1. The two parameterizations\n", + "## 1. Same distribution, different coordinates\n", "\n", - "For a group-specific intercept $u_g$ with group-level mean $\\mu$ and group-level scale $\\sigma$:\n", + "Let $u_g$ be a coefficient for group $g$, with population location\n", + "$\\mu$ and group scale $\\sigma$.\n", "\n", - "**Centered (`noncentered=False`)** — sample the group effect directly:\n", + "In the **centered** parameterization we sample the coefficient directly:\n", "\n", - "$$u_g \\sim \\mathcal{N}(\\mu, \\sigma)$$\n", + "$$\n", + "u_g \\sim \\mathcal{N}(\\mu, \\sigma).\n", + "$$\n", "\n", - "**Non-centered (`noncentered=True`, the bambi default)** — sample a standard-normal offset and rescale:\n", + "A mathematically equivalent **non-centered** parameterization samples a\n", + "standard-normal coordinate and transforms it:\n", "\n", - "$$z_g \\sim \\mathcal{N}(0, 1), \\qquad u_g = z_g \\cdot \\sigma$$\n", + "$$\n", + "z_g \\sim \\mathcal{N}(0, 1),\n", + "\\qquad\n", + "u_g = \\mu + \\sigma z_g.\n", + "$$\n", "\n", - "Both produce the same prior distribution on $u_g$ when $\\mu = 0$. The non-centered form usually samples better because it decouples $u_g$ from $\\sigma$ in the geometry of the posterior. **But notice that the non-centered form does not use $\\mu$ at all.** That is the source of the footgun we will see below: a `mu` hyperprior you supply on the group term is silently ignored under non-centered." + "These equations define the same prior distribution for $u_g$. Which\n", + "coordinates sample better depends on the amount of information in the data\n", + "and the posterior geometry; neither form is universally superior." ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "bkHC", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "## 2. Setup" + "> **The current Bambi boundary**\n", + "\n", + "For a non-centered group term, Bambi currently constructs\n", + "\n", + "$$\n", + "u_g = \\sigma z_g,\n", + "$$\n", + "\n", + "rather than $\\mu + \\sigma z_g$. This shortcut is faithful when the group\n", + "term is a zero-mean deviation: `mu` is absent or fixed entirely to zero. A\n", + "free or nonzero `mu` would be dropped.\n", + "\n", + "HSSM validates explicit group priors before asking Bambi to build the PyMC\n", + "model. If that shortcut would discard part of the requested prior, HSSM\n", + "raises a `ValueError` instead of silently changing the model." ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, + "id": "lEQa", "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
participant_idstimrtresponsethetadbsconf
00LL1.211.00.6562751HC
10WL1.631.0-0.3278891LC
20WW1.031.0-0.4802851HC
30WL2.771.01.9274271LC
40WW1.14-1.0-0.2132361HC
\n", - "
" - ], - "text/plain": [ - " participant_id stim rt response theta dbs conf\n", - "0 0 LL 1.21 1.0 0.656275 1 HC\n", - "1 0 WL 1.63 1.0 -0.327889 1 LC\n", - "2 0 WW 1.03 1.0 -0.480285 1 HC\n", - "3 0 WL 2.77 1.0 1.927427 1 LC\n", - "4 0 WW 1.14 -1.0 -0.213236 1 HC" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "import logging\n", - "import os\n", - "import shutil\n", - "import warnings\n", + "tutorial_data = pd.DataFrame(\n", + " {\n", + " \"rt\": 0.38 + 0.015 * np.arange(16),\n", + " \"response\": np.where(np.arange(16) % 2, 1, -1),\n", + " \"theta\": np.linspace(-1.0, 1.0, 16),\n", + " \"participant_id\": np.repeat(np.arange(4), 4),\n", + " }\n", + ")\n", + "raw_bambi_data = tutorial_data.loc[:, [\"participant_id\"]].assign(\n", + " y=np.array(\n", + " [\n", + " -0.20,\n", + " 0.05,\n", + " 0.10,\n", + " -0.05,\n", + " 0.30,\n", + " 0.35,\n", + " 0.20,\n", + " 0.40,\n", + " -0.35,\n", + " -0.20,\n", + " -0.10,\n", + " -0.25,\n", + " 0.15,\n", + " 0.25,\n", + " 0.05,\n", + " 0.30,\n", + " ]\n", + " )\n", + ")\n", "\n", - "# NOTE: The block below is a niche workaround that most users will not need.\n", - "# It only matters on macOS when the notebook kernel was launched from a\n", - "# non-terminal context (some IDE integrations, GUI launchers) so that\n", - "# `/opt/homebrew/bin` is missing from PATH. In that situation\n", - "# `pm.model_to_graphviz` raises `ExecutableNotFound: PosixPath('dot')` even\n", - "# though Graphviz is installed. We prepend the most common locations only if\n", - "# `dot` is not already on PATH, so this is a no-op for everyone else.\n", - "if shutil.which(\"dot\") is None:\n", - " for _candidate in (\"/opt/homebrew/bin\", \"/usr/local/bin\", \"/opt/conda/bin\"):\n", - " if os.path.isfile(os.path.join(_candidate, \"dot\")):\n", - " os.environ[\"PATH\"] = _candidate + os.pathsep + os.environ[\"PATH\"]\n", - " break\n", "\n", - "import pymc as pm\n", + "def zero_mean_group_prior(*, noncentered=None):\n", + " \"\"\"Return a fresh zero-mean hierarchical Normal group prior.\"\"\"\n", + " return hssm.Prior(\n", + " \"Normal\",\n", + " mu=0.0,\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + " noncentered=noncentered,\n", + " )\n", "\n", - "import hssm\n", "\n", - "warnings.filterwarnings(\"ignore\")\n", + "def free_location_group_prior(*, noncentered=None):\n", + " \"\"\"Return a fresh hierarchical Normal with a free population mean.\"\"\"\n", + " return hssm.Prior(\n", + " \"Normal\",\n", + " mu=hssm.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + " noncentered=noncentered,\n", + " )\n", "\n", - "# Show HSSM warnings inline so we can see the footgun message later.\n", - "logging.basicConfig(level=logging.INFO, format=\"%(levelname)s %(name)s: %(message)s\")\n", - "logging.getLogger(\"hssm\").setLevel(logging.WARNING)\n", "\n", - "cav_data = hssm.load_data(\"cavanagh_theta\")\n", - "cav_data.head()" + "def matched_include(group_prior):\n", + " \"\"\"Place a common intercept and matching group intercept in one spec.\"\"\"\n", + " return [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", + " \"prior\": {\n", + " \"Intercept\": hssm.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " \"1|participant_id\": group_prior,\n", + " },\n", + " }\n", + " ]\n", + "\n", + "\n", + "def hssm_model(*, formula, priors, noncentered=True):\n", + " \"\"\"Build a tiny analytical DDM without sampling or init-value work.\"\"\"\n", + " return hssm.HSSM(\n", + " data=tutorial_data,\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=[{\"name\": \"v\", \"formula\": formula, \"prior\": priors}],\n", + " p_outlier=0.0,\n", + " prior_settings=None,\n", + " noncentered=noncentered,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + " )\n", + "\n", + "\n", + "def disconnected_free_rvs(pymc_model):\n", + " \"\"\"Return free-RV names that are not ancestors of observed variables.\"\"\"\n", + " connected = {\n", + " id(variable)\n", + " for observed_rv in pymc_model.observed_RVs\n", + " for variable in ancestors([observed_rv])\n", + " }\n", + " return sorted(rv.name for rv in pymc_model.free_RVs if id(rv) not in connected)\n", + "\n", + "\n", + "def capture_hssm_build(builder):\n", + " \"\"\"Build a model while collecting HSSM warning messages.\"\"\"\n", + " logger = logging.getLogger(\"hssm\")\n", + " messages = []\n", + "\n", + " class _MessageHandler(logging.Handler):\n", + " def emit(self, record):\n", + " messages.append(record.getMessage())\n", + "\n", + " handler = _MessageHandler(level=logging.WARNING)\n", + " previous_handlers = list(logger.handlers)\n", + " previous_level = logger.level\n", + " previous_propagate = logger.propagate\n", + " logger.handlers = [handler]\n", + " logger.setLevel(logging.WARNING)\n", + " logger.propagate = False\n", + " try:\n", + " result = builder()\n", + " finally:\n", + " logger.handlers = previous_handlers\n", + " logger.setLevel(previous_level)\n", + " logger.propagate = previous_propagate\n", + " return result, tuple(messages)\n", + "\n", + "\n", + "# Assert the factories return fresh prior trees. HSSM/Bambi attach names\n", + "# while preparing priors, so tutorial cases must not share mutable objects.\n", + "assert zero_mean_group_prior() is not zero_mean_group_prior()\n", + "assert free_location_group_prior() is not free_location_group_prior()\n", + "assert isinstance(free_location_group_prior().args[\"mu\"], bmb.Prior)" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "PKri", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "## 3. A centered model\n", + "## 2. A valid zero-mean group deviation\n", + "\n", + "Consider\n", + "\n", + "```python\n", + "v ~ 1 + (1 | participant_id)\n", + "```\n", "\n", - "We fit a DDM where the drift rate `v` has a participant-level intercept. We pass `noncentered=False` so bambi builds the term as `1|participant_id ~ Normal(mu, sigma)`." + "The common `Intercept` owns the population location. The participant term\n", + "is a deviation around that location, so its prior has `mu=0` and a\n", + "hierarchical `sigma`. This statistical model can be represented faithfully\n", + "in either parameterization." ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, + "id": "Xref", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "['t',\n", - " 'z',\n", - " 'a',\n", - " 'v_Intercept',\n", - " 'v_1|participant_id_sigma',\n", - " 'v_1|participant_id']" + "text/html": [ + "
effective parameterizationsampled group coordinategroup coefficient in graphdisconnected free RVs
0centeredv_1|participant_idfree RVnone
1non-centeredv_1|participant_id_offsetdeterministicnone
" ] }, - "execution_count": 2, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_centered = hssm.HSSM(\n", - " data=cav_data,\n", - " model=\"ddm\",\n", - " include=[\n", - " {\n", - " \"name\": \"v\",\n", - " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", - " \"prior\": {\n", - " \"Intercept\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 1.5},\n", - " \"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": 0.0,\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - " },\n", - " },\n", - " }\n", - " ],\n", - " p_outlier=0.0,\n", - " noncentered=False,\n", + "matched_centered_model, matched_centered_messages = capture_hssm_build(\n", + " lambda: hssm.HSSM(\n", + " data=tutorial_data,\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=matched_include(zero_mean_group_prior()),\n", + " p_outlier=0.0,\n", + " prior_settings=None,\n", + " noncentered=False,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + " )\n", + ")\n", + "matched_noncentered_model, matched_noncentered_messages = capture_hssm_build(\n", + " lambda: hssm.HSSM(\n", + " data=tutorial_data,\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=matched_include(zero_mean_group_prior()),\n", + " p_outlier=0.0,\n", + " prior_settings=None,\n", + " noncentered=True,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + " )\n", ")\n", "\n", - "[rv.name for rv in model_centered.pymc_model.free_RVs]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Under the centered parameterization the group-specific term shows up as a single free RV named `v_1|participant_id`, plus its hyperprior `v_1|participant_id_sigma`. We can confirm this visually with the PyMC graph:" + "_centered_names = {rv.name for rv in matched_centered_model.pymc_model.free_RVs}\n", + "_noncentered_names = {rv.name for rv in matched_noncentered_model.pymc_model.free_RVs}\n", + "assert matched_centered_messages == ()\n", + "assert matched_noncentered_messages == ()\n", + "assert \"v_1|participant_id\" in _centered_names\n", + "assert \"v_1|participant_id_offset\" not in _centered_names\n", + "assert \"v_1|participant_id_offset\" in _noncentered_names\n", + "assert \"v_1|participant_id\" not in _noncentered_names\n", + "assert \"v_1|participant_id_mu\" not in _centered_names | _noncentered_names\n", + "assert disconnected_free_rvs(matched_centered_model.pymc_model) == []\n", + "assert disconnected_free_rvs(matched_noncentered_model.pymc_model) == []\n", + "\n", + "matched_parameterization_table = pd.DataFrame(\n", + " [\n", + " {\n", + " \"effective parameterization\": \"centered\",\n", + " \"sampled group coordinate\": \"v_1|participant_id\",\n", + " \"group coefficient in graph\": \"free RV\",\n", + " \"disconnected free RVs\": \"none\",\n", + " },\n", + " {\n", + " \"effective parameterization\": \"non-centered\",\n", + " \"sampled group coordinate\": \"v_1|participant_id_offset\",\n", + " \"group coefficient in graph\": \"deterministic\",\n", + " \"disconnected free RVs\": \"none\",\n", + " },\n", + " ]\n", + ")\n", + "matched_parameterization_table" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, + "id": "SFPL", "metadata": {}, "outputs": [ { @@ -251,215 +391,165 @@ "\n", "\n", - "\n", "\n", - "\n", + "\n", "\n", - "\n", + "\n", "\n", - "clusterv_1|participant_id__factor_dim (14)\n", - "\n", - "v_1|participant_id__factor_dim (14)\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", "\n", "\n", - "cluster__obs__ (3988)\n", - "\n", - "__obs__ (3988)\n", + "cluster__obs__ (16)\n", + "\n", + "__obs__ (16)\n", "\n", "\n", - "cluster__obs__ (3988) x rt,response_extra_dim_0 (2)\n", - "\n", - "__obs__ (3988) x rt,response_extra_dim_0 (2)\n", + "cluster__obs__ (16) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (16) x rt,response_extra_dim_0 (2)\n", "\n", - "\n", + "\n", "\n", + "z\n", + "\n", + "z\n", + "~\n", + "Uniform\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "z->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", "v_Intercept\n", - "\n", - "v_Intercept\n", - "~\n", - "Normal\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v\n", - "\n", - "v\n", - "~\n", - "Deterministic\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", "\n", "\n", "\n", "v_Intercept->v\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "a\n", - "\n", - "a\n", - "~\n", - "Halfnormal\n", + "\n", + "\n", "\n", - "\n", - "\n", - "rt,response\n", - "\n", - "rt,response\n", - "~\n", - "Ddm_RV\n", + "\n", + "\n", + "t\n", + "\n", + "t\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "a->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "t->rt,response\n", + "\n", + "\n", "\n", "\n", - "\n", + "\n", "v_1|participant_id_sigma\n", - "\n", - "v_1|participant_id_sigma\n", - "~\n", - "Halfnormal\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", "\n", "\n", "\n", "v_1|participant_id\n", - "\n", - "v_1|participant_id\n", - "~\n", - "Normal\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v_1|participant_id_sigma->v_1|participant_id\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "t\n", - "\n", - "t\n", - "~\n", - "Halfnormal\n", - "\n", - "\n", - "\n", - "t->rt,response\n", - "\n", - "\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "z\n", - "\n", - "z\n", - "~\n", - "Uniform\n", + "a\n", + "\n", + "a\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "z->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "a->rt,response\n", + "\n", + "\n", "\n", "\n", "\n", "v_1|participant_id->v\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n", "v->rt,response\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n" - ], - "text/plain": [ - "" ] }, - "execution_count": 3, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_centered.graph()" + "matched_centered_graph = pm.model_to_graphviz(matched_centered_model.pymc_model)\n", + "matched_centered_graph" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. A non-centered model\n", - "\n", - "Same prior dictionary, but `noncentered=True` (which is also the default if you do not pass the kwarg)." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "['t',\n", - " 'z',\n", - " 'a',\n", - " 'v_Intercept',\n", - " 'v_1|participant_id_sigma',\n", - " 'v_1|participant_id_offset']" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" + "id": "BYtC", + "metadata": { + "marimo": { + "md_prefix": "" } - ], - "source": [ - "model_noncentered = hssm.HSSM(\n", - " data=cav_data,\n", - " model=\"ddm\",\n", - " include=[\n", - " {\n", - " \"name\": \"v\",\n", - " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", - " \"prior\": {\n", - " \"Intercept\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 1.5},\n", - " \"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": 0.0,\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - " },\n", - " },\n", - " }\n", - " ],\n", - " p_outlier=0.0,\n", - " noncentered=True,\n", - ")\n", - "\n", - "[rv.name for rv in model_noncentered.pymc_model.free_RVs]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + }, "source": [ - "Now there is a new free RV called `v_1|participant_id_offset` (the standardized $z_g$), and `v_1|participant_id` is a *deterministic* equal to `offset * sigma`. The naming convention has changed." + "In the centered graph, `v_1|participant_id` is sampled directly. Its scale\n", + "hyperprior points to the group coefficients, and those coefficients point\n", + "to the trial-wise drift rate." ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, + "id": "RGSE", "metadata": {}, "outputs": [ { @@ -468,265 +558,254 @@ "\n", "\n", - "\n", "\n", - "\n", + "\n", "\n", - "\n", + "\n", "\n", - "clusterv_1|participant_id__factor_dim (14)\n", - "\n", - "v_1|participant_id__factor_dim (14)\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", "\n", "\n", - "cluster__obs__ (3988)\n", - "\n", - "__obs__ (3988)\n", + "cluster__obs__ (16)\n", + "\n", + "__obs__ (16)\n", "\n", "\n", - "cluster__obs__ (3988) x rt,response_extra_dim_0 (2)\n", - "\n", - "__obs__ (3988) x rt,response_extra_dim_0 (2)\n", + "cluster__obs__ (16) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (16) x rt,response_extra_dim_0 (2)\n", "\n", - "\n", + "\n", "\n", + "z\n", + "\n", + "z\n", + "~\n", + "Uniform\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "z->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", "v_Intercept\n", - "\n", - "v_Intercept\n", - "~\n", - "Normal\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v\n", - "\n", - "v\n", - "~\n", - "Deterministic\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", "\n", "\n", "\n", "v_Intercept->v\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "a\n", - "\n", - "a\n", - "~\n", - "Halfnormal\n", + "\n", + "\n", "\n", - "\n", - "\n", - "rt,response\n", - "\n", - "rt,response\n", - "~\n", - "Ddm_RV\n", + "\n", + "\n", + "t\n", + "\n", + "t\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "a->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "t->rt,response\n", + "\n", + "\n", "\n", "\n", - "\n", + "\n", "v_1|participant_id_sigma\n", - "\n", - "v_1|participant_id_sigma\n", - "~\n", - "Halfnormal\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", "\n", "\n", "\n", "v_1|participant_id\n", - "\n", - "v_1|participant_id\n", - "~\n", - "Deterministic\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Deterministic\n", "\n", "\n", "\n", "v_1|participant_id_sigma->v_1|participant_id\n", - "\n", - "\n", + "\n", + "\n", "\n", - "\n", - "\n", - "t\n", - "\n", - "t\n", - "~\n", - "Halfnormal\n", - "\n", - "\n", - "\n", - "t->rt,response\n", - "\n", - "\n", - "\n", - "\n", + "\n", "\n", - "z\n", - "\n", - "z\n", - "~\n", - "Uniform\n", + "a\n", + "\n", + "a\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "z->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "a->rt,response\n", + "\n", + "\n", "\n", "\n", "\n", "v_1|participant_id_offset\n", - "\n", - "v_1|participant_id_offset\n", - "~\n", - "Normal\n", + "\n", + "v_1|participant_id_offset\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v_1|participant_id_offset->v_1|participant_id\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n", "v_1|participant_id->v\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n", "v->rt,response\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n" - ], - "text/plain": [ - "" ] }, - "execution_count": 5, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "pm.model_to_graphviz(model_noncentered.pymc_model)" + "matched_noncentered_graph = pm.model_to_graphviz(matched_noncentered_model.pymc_model)\n", + "matched_noncentered_graph" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "Kclp", + "metadata": { + "marimo": { + "md_prefix": "" + } + }, + "source": [ + "In the non-centered graph, Bambi samples a standard-normal `offset` and\n", + "combines it with `sigma`. The resulting `v_1|participant_id` is a\n", + "deterministic node. Because the intended group mean is exactly zero, this\n", + "graph represents the same prior as the centered graph above." + ] + }, + { + "cell_type": "markdown", + "id": "emfo", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "## 5. The footgun: a `mu` hyperprior under non-centered\n", + "## 3. Why HSSM needs a pre-build guard\n", "\n", - "Suppose you write the group prior with a *hyperprior* on the group mean — perfectly natural if you are coming from a centered, fully Bayesian mindset:\n", + "To isolate the underlying mechanism, the next example uses Bambi directly\n", + "with a tiny Gaussian response. The requested group prior has a free\n", + "population location:\n", "\n", "```python\n", - "\"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 0.5}, # <-- hyperprior on mu\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - "}\n", + "Normal(\n", + " mu=Normal(0, 0.5),\n", + " sigma=HalfNormal(0.5),\n", + ")\n", "```\n", "\n", - "Under `noncentered=True` bambi reparameterizes the term as `offset * sigma` and **never uses `mu`**. The `mu` hyperprior is still created in the PyMC graph (as `v_1|participant_id_mu`), but it is a *disconnected* free RV: it has no path to the observed data. HSSM now warns about this twice — once with a targeted, actionable message, and once with a general \"these RVs are not wired in\" report." + "This is a legitimate centered hierarchy. Under Bambi's current\n", + "non-centered shortcut, however, the `mu` prior is created and then omitted\n", + "from `offset * sigma`." ] }, { "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "User prior for '1|participant_id' on parameter 'v' supplies a hyperprior on `mu`, but the effective parameterization is non-centered. bambi will reparameterize this term as `offset * sigma` and drop the `mu` hyperprior, leaving it as a disconnected node in the PyMC graph. Either pass `noncentered=False` to `HSSM(...)` so that `mu` is used in the centered Normal, or move the location prior to the common `Intercept` (e.g. use a formula like 'v ~ 1 + (1|participant_id)' and attach the `mu` prior to 'Intercept'). To silence this warning without changing the model, set the `mu` argument to a scalar (e.g. `mu=0`).\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING hssm: User prior for '1|participant_id' on parameter 'v' supplies a hyperprior on `mu`, but the effective parameterization is non-centered. bambi will reparameterize this term as `offset * sigma` and drop the `mu` hyperprior, leaving it as a disconnected node in the PyMC graph. Either pass `noncentered=False` to `HSSM(...)` so that `mu` is used in the centered Normal, or move the location prior to the common `Intercept` (e.g. use a formula like 'v ~ 1 + (1|participant_id)' and attach the `mu` prior to 'Intercept'). To silence this warning without changing the model, set the `mu` argument to a scalar (e.g. `mu=0`).\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The PyMC graph contains free random variables that do not influence the likelihood: 'v_1|participant_id_mu'. This typically happens when a hyperprior is supplied for a parameter that the chosen parameterization does not use (e.g. `mu` under `noncentered=True`). These nodes will be sampled but will not affect inference; consider switching the parameterization or adjusting the prior specification.\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING hssm: The PyMC graph contains free random variables that do not influence the likelihood: 'v_1|participant_id_mu'. This typically happens when a hyperprior is supplied for a parameter that the chosen parameterization does not use (e.g. `mu` under `noncentered=True`). These nodes will be sampled but will not affect inference; consider switching the parameterization or adjusting the prior specification.\n" - ] - } - ], - "source": [ - "model_footgun = hssm.HSSM(\n", - " data=cav_data,\n", - " model=\"ddm\",\n", - " include=[\n", - " {\n", - " \"name\": \"v\",\n", - " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", - " \"prior\": {\n", - " \"Intercept\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 1.5},\n", - " \"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 0.5},\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - " },\n", - " },\n", - " }\n", - " ],\n", - " p_outlier=0.0,\n", - " noncentered=True,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 7, + "execution_count": null, + "id": "Hstk", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "['v_1|participant_id_mu']" + "text/html": [ + "
quantityvalue
0sampled group coordinatessigma, Intercept, 1|participant_id_mu, 1|participant_id_sigma, 1|participant_id_offset
1disconnected free RVs1|participant_id_mu
" ] }, - "execution_count": 7, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "from hssm.param.parameterization_check import find_disconnected_free_rvs\n", + "raw_bambi_group_prior = bmb.Prior(\n", + " \"Normal\",\n", + " mu=bmb.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " sigma=bmb.Prior(\"HalfNormal\", sigma=0.5),\n", + ")\n", + "raw_bambi_model = bmb.Model(\n", + " \"y ~ 1 + (1|participant_id)\",\n", + " raw_bambi_data,\n", + " family=\"gaussian\",\n", + " priors={\"1|participant_id\": raw_bambi_group_prior},\n", + " noncentered=True,\n", + ")\n", + "raw_bambi_model.build()\n", + "raw_bambi_pymc_model = raw_bambi_model.backend.model\n", + "raw_bambi_orphans = disconnected_free_rvs(raw_bambi_pymc_model)\n", + "\n", + "_free_names = [rv.name for rv in raw_bambi_pymc_model.free_RVs]\n", + "assert \"1|participant_id_mu\" in _free_names\n", + "assert \"1|participant_id_offset\" in _free_names\n", + "assert raw_bambi_orphans == [\"1|participant_id_mu\"]\n", "\n", - "find_disconnected_free_rvs(model_footgun.pymc_model)" + "raw_bambi_summary = pd.DataFrame(\n", + " {\n", + " \"quantity\": [\"sampled group coordinates\", \"disconnected free RVs\"],\n", + " \"value\": [\", \".join(_free_names), \", \".join(raw_bambi_orphans)],\n", + " }\n", + ")\n", + "raw_bambi_summary" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, + "id": "nWHF", "metadata": {}, "outputs": [ { @@ -735,279 +814,296 @@ "\n", "\n", - "\n", "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", "\n", - "clusterv_1|participant_id__factor_dim (14)\n", - "\n", - "v_1|participant_id__factor_dim (14)\n", + "clusterparticipant_id__factor_dim (4)\n", + "\n", + "participant_id__factor_dim (4)\n", "\n", "\n", - "cluster__obs__ (3988)\n", - "\n", - "__obs__ (3988)\n", + "cluster__obs__ (16)\n", + "\n", + "__obs__ (16)\n", "\n", - "\n", - "cluster__obs__ (3988) x rt,response_extra_dim_0 (2)\n", - "\n", - "__obs__ (3988) x rt,response_extra_dim_0 (2)\n", - "\n", - "\n", + "\n", "\n", - "v_Intercept\n", - "\n", - "v_Intercept\n", - "~\n", - "Normal\n", + "Intercept\n", + "\n", + "Intercept\n", + "~\n", + "Normal\n", "\n", - "\n", - "\n", - "v\n", - "\n", - "v\n", - "~\n", - "Deterministic\n", + "\n", + "\n", + "mu\n", + "\n", + "mu\n", + "~\n", + "Deterministic\n", "\n", - "\n", + "\n", "\n", - "v_Intercept->v\n", - "\n", - "\n", + "Intercept->mu\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "a\n", - "\n", - "a\n", - "~\n", - "Halfnormal\n", + "1|participant_id_sigma\n", + "\n", + "1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "rt,response\n", - "\n", - "rt,response\n", - "~\n", - "Ddm_RV\n", + "\n", + "\n", + "1|participant_id\n", + "\n", + "1|participant_id\n", + "~\n", + "Deterministic\n", "\n", - "\n", - "\n", - "a->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "1|participant_id_sigma->1|participant_id\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "v_1|participant_id_sigma\n", - "\n", - "v_1|participant_id_sigma\n", - "~\n", - "Halfnormal\n", + "1|participant_id_mu\n", + "\n", + "1|participant_id_mu\n", + "~\n", + "Normal\n", "\n", - "\n", - "\n", - "v_1|participant_id\n", - "\n", - "v_1|participant_id\n", - "~\n", - "Deterministic\n", - "\n", - "\n", - "\n", - "v_1|participant_id_sigma->v_1|participant_id\n", - "\n", - "\n", - "\n", - "\n", + "\n", "\n", - "t\n", - "\n", - "t\n", - "~\n", - "Halfnormal\n", - "\n", - "\n", - "\n", - "t->rt,response\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "v_1|participant_id_mu\n", - "\n", - "v_1|participant_id_mu\n", - "~\n", - "Normal\n", + "sigma\n", + "\n", + "sigma\n", + "~\n", + "Halfstudentt\n", "\n", - "\n", - "\n", - "z\n", - "\n", - "z\n", - "~\n", - "Uniform\n", + "\n", + "\n", + "y\n", + "\n", + "y\n", + "~\n", + "Normal\n", "\n", - "\n", + "\n", "\n", - "z->rt,response\n", - "\n", - "\n", + "sigma->y\n", + "\n", + "\n", "\n", - "\n", - "\n", - "v_1|participant_id_offset\n", - "\n", - "v_1|participant_id_offset\n", - "~\n", - "Normal\n", + "\n", + "\n", + "1|participant_id_offset\n", + "\n", + "1|participant_id_offset\n", + "~\n", + "Normal\n", "\n", - "\n", + "\n", "\n", - "v_1|participant_id_offset->v_1|participant_id\n", - "\n", - "\n", + "1|participant_id_offset->1|participant_id\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "v_1|participant_id->v\n", - "\n", - "\n", + "1|participant_id->mu\n", + "\n", + "\n", "\n", - "\n", - "\n", - "v->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "mu->y\n", + "\n", + "\n", "\n", "\n", "\n" - ], - "text/plain": [ - "" ] }, - "execution_count": 8, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "pm.model_to_graphviz(model_footgun.pymc_model)" + "raw_bambi_orphan_graph = pm.model_to_graphviz(raw_bambi_pymc_model)\n", + "raw_bambi_orphan_graph" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "Notice the floating `v_1|participant_id_mu` node — it has no arrow pointing toward the response. That is the orphan." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "id": "iLit", + "metadata": { + "marimo": { + "md_prefix": "" + } + }, "source": [ - "### Fix 1 — switch to the centered parameterization\n", - "\n", - "If you really want a hyperprior on the group mean, use `noncentered=False` so that `mu` is wired in." + "> The floating `1|participant_id_mu` node has no path to the observed\n", + "> response. Sampling it would consume computation without changing the\n", + "> likelihood. This graph documents the **raw Bambi behavior that HSSM\n", + "> prevents**; it is not a graph that current HSSM will construct from the\n", + "> same explicit prior." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, + "id": "ZHCJ", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "[]" + "text/markdown": [ + "**Current HSSM result**\n", + "\n", + "```text\n", + "Explicit group-specific prior specification(s) cannot be represented faithfully by bambi:\n", + "- User prior for group term '1|participant_id' on parameter 'v' is incompatible with the effective parameterization: outer prior supplies a `mu` hyperprior that bambi creates and then omits from `offset * sigma`, leaving a disconnected node. Continuing would either fail in bambi or change the requested prior. Keep the common formula term '1' and use a plain built-in Normal with hierarchical `sigma`, absent or fixed-all-zero `mu`, and no additional arguments for its zero-mean group deviation. To retain the explicit prior instead, set `noncentered=False` on this prior and on any nested hierarchical hyperpriors (or remove their overrides and make the effective component setting for 'v' centered); remove the common effect as well if a free group mean should own the population location.\n", + "```" ] }, - "execution_count": 9, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_fix_centered = hssm.HSSM(\n", - " data=cav_data,\n", - " model=\"ddm\",\n", - " include=[\n", - " {\n", - " \"name\": \"v\",\n", - " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", - " \"prior\": {\n", - " \"Intercept\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 1.5},\n", - " \"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 0.5},\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - " },\n", - " },\n", - " }\n", - " ],\n", - " p_outlier=0.0,\n", - " noncentered=False,\n", - ")\n", + "try:\n", + " hssm.HSSM(\n", + " data=tutorial_data,\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=matched_include(free_location_group_prior()),\n", + " p_outlier=0.0,\n", + " prior_settings=None,\n", + " noncentered=True,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + " )\n", + "except ValueError as exc:\n", + " hssm_preflight_error = str(exc)\n", + "else:\n", + " raise AssertionError(\"HSSM did not reject the incompatible group prior\")\n", + "\n", + "assert \"cannot be represented faithfully by bambi\" in hssm_preflight_error\n", + "assert \"1|participant_id\" in hssm_preflight_error\n", + "assert \"omits from `offset * sigma`\" in hssm_preflight_error\n", + "assert \"noncentered=False\" in hssm_preflight_error\n", + "mo.md(f\"\"\"\n", + "**Current HSSM result**\n", + "\n", + "```text\n", + "{hssm_preflight_error}\n", + "```\n", + "\"\"\")" + ] + }, + { + "cell_type": "markdown", + "id": "ROlb", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, + "source": [ + "## 4. Centering retains the mean—but location ownership still matters\n", + "\n", + "If we center the explicit hierarchy, Bambi uses its `mu`. With both a\n", + "common intercept and a free group mean, the predictor for participant $g$\n", + "contains\n", "\n", - "find_disconnected_free_rvs(model_fix_centered.pymc_model)" + "$$\n", + "\\eta_g = \\beta_0 + u_g,\n", + "\\qquad\n", + "u_g \\sim \\mathcal{N}(\\mu_u, \\sigma_u).\n", + "$$\n", + "\n", + "The likelihood sees $\\beta_0 + \\mu_u$, not the two locations separately.\n", + "Shifting one up and the other down leaves the predictor unchanged. The PyMC\n", + "graph is fully connected, but the likelihood has a ridge along that shift\n", + "direction." ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, + "id": "qnkX", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "Hierarchical Sequential Sampling Model\n", - "Model: ddm\n", + "text/markdown": [ + "**HSSM's centered location warning**\n", "\n", - "Response variable: rt,response\n", - "Likelihood: analytical\n", - "Observations: 3988\n", + "```text\n", + "User prior for '1|participant_id' on parameter 'v' has a free `mu`, and its Formulae expression 'Intercept' also occurs as a common effect under the effective centered parameterization. The data only constrains their sum; the common and group locations are non-identifiable individually and the posterior will have a ridge along the anti-diagonal. Keep the common 'Intercept' effect and set `mu=0` on the matching group term, or remove that common effect if the group-level mean should own the location.\n", + "```\n", "\n", - "Parameters:\n", - "\n", - "v:\n", - " Formula: v ~ 1 + (1|participant_id)\n", - " Priors:\n", - " v_Intercept ~ Normal(mu: 0.0, sigma: 1.5)\n", - " v_1|participant_id ~ Normal(mu: Normal(mu: 0.0, sigma: 0.5), sigma: HalfNormal(sigma: 0.5))\n", - " Link: identity\n", - " Explicit bounds: (-inf, inf)\n", - "\n", - "a:\n", - " Prior: HalfNormal(sigma: 2.0)\n", - " Explicit bounds: (0.0, inf)\n", - "\n", - "z:\n", - " Prior: Uniform(lower: 0.0, upper: 1.0)\n", - " Explicit bounds: (0.0, 1.0)\n", - "\n", - "t:\n", - " Prior: HalfNormal(sigma: 2.0)\n", - " Explicit bounds: (0.0, inf)" + "The disconnected-node check returns `[]`: this is an identifiability\n", + "problem in the likelihood, not an orphan-node problem." ] }, - "execution_count": 13, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_fix_centered" + "centered_ridge_model, centered_ridge_messages = capture_hssm_build(\n", + " lambda: hssm.HSSM(\n", + " data=tutorial_data,\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=matched_include(free_location_group_prior()),\n", + " p_outlier=0.0,\n", + " prior_settings=None,\n", + " noncentered=False,\n", + " process_initvals=False,\n", + " initval_jitter=0.0,\n", + " )\n", + ")\n", + "centered_ridge_warning = next(\n", + " message\n", + " for message in centered_ridge_messages\n", + " if \"non-identifiable individually\" in message\n", + ")\n", + "\n", + "assert \"posterior will have a ridge\" in centered_ridge_warning\n", + "assert \"common 'Intercept' effect\" in centered_ridge_warning\n", + "assert disconnected_free_rvs(centered_ridge_model.pymc_model) == []\n", + "assert \"v_1|participant_id_mu\" in {\n", + " rv.name for rv in centered_ridge_model.pymc_model.free_RVs\n", + "}\n", + "mo.md(f\"\"\"\n", + "**HSSM's centered location warning**\n", + "\n", + "```text\n", + "{centered_ridge_warning}\n", + "```\n", + "\n", + "The disconnected-node check returns `[]`: this is an identifiability\n", + "problem in the likelihood, not an orphan-node problem.\n", + "\"\"\")" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, + "id": "TqIu", "metadata": {}, "outputs": [ { @@ -1016,320 +1112,499 @@ "\n", "\n", - "\n", "\n", - "\n", + "\n", "\n", - "\n", + "\n", "\n", - "clusterv_1|participant_id__factor_dim (14)\n", - "\n", - "v_1|participant_id__factor_dim (14)\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", "\n", "\n", - "cluster__obs__ (3988)\n", - "\n", - "__obs__ (3988)\n", + "cluster__obs__ (16)\n", + "\n", + "__obs__ (16)\n", "\n", "\n", - "cluster__obs__ (3988) x rt,response_extra_dim_0 (2)\n", - "\n", - "__obs__ (3988) x rt,response_extra_dim_0 (2)\n", + "cluster__obs__ (16) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (16) x rt,response_extra_dim_0 (2)\n", "\n", - "\n", + "\n", "\n", + "z\n", + "\n", + "z\n", + "~\n", + "Uniform\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "z->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", "v_Intercept\n", - "\n", - "v_Intercept\n", - "~\n", - "Normal\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v\n", - "\n", - "v\n", - "~\n", - "Deterministic\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", "\n", "\n", "\n", "v_Intercept->v\n", - "\n", - "\n", + "\n", + "\n", "\n", - "\n", - "\n", - "a\n", - "\n", - "a\n", - "~\n", - "Halfnormal\n", - "\n", - "\n", - "\n", - "rt,response\n", - "\n", - "rt,response\n", - "~\n", - "Ddm_RV\n", - "\n", - "\n", - "\n", - "a->rt,response\n", - "\n", - "\n", - "\n", - "\n", + "\n", "\n", - "v_1|participant_id_sigma\n", - "\n", - "v_1|participant_id_sigma\n", - "~\n", - "Halfnormal\n", + "v_1|participant_id_mu\n", + "\n", + "v_1|participant_id_mu\n", + "~\n", + "Normal\n", "\n", "\n", "\n", "v_1|participant_id\n", - "\n", - "v_1|participant_id\n", - "~\n", - "Normal\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Normal\n", "\n", - "\n", - "\n", - "v_1|participant_id_sigma->v_1|participant_id\n", - "\n", - "\n", + "\n", + "\n", + "v_1|participant_id_mu->v_1|participant_id\n", + "\n", + "\n", "\n", "\n", "\n", "t\n", - "\n", - "t\n", - "~\n", - "Halfnormal\n", + "\n", + "t\n", + "~\n", + "Halfnormal\n", "\n", "\n", - "\n", + "\n", "t->rt,response\n", - "\n", - "\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "v_1|participant_id_mu\n", - "\n", - "v_1|participant_id_mu\n", - "~\n", - "Normal\n", + "v_1|participant_id_sigma\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "v_1|participant_id_mu->v_1|participant_id\n", - "\n", - "\n", + "\n", + "\n", + "v_1|participant_id_sigma->v_1|participant_id\n", + "\n", + "\n", "\n", - "\n", + "\n", "\n", - "z\n", - "\n", - "z\n", - "~\n", - "Uniform\n", + "a\n", + "\n", + "a\n", + "~\n", + "Halfnormal\n", "\n", - "\n", - "\n", - "z->rt,response\n", - "\n", - "\n", + "\n", + "\n", + "a->rt,response\n", + "\n", + "\n", "\n", "\n", "\n", "v_1|participant_id->v\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n", "v->rt,response\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n" - ], - "text/plain": [ - "" ] }, - "execution_count": 11, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_fix_centered.graph()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "> **Note — Fix 1 trades one problem for another.** Switching to `noncentered=False` removes the disconnected `_mu` node and the graph now looks healthy. But look at the linear predictor for participant `g`:\n", - ">\n", - "> ```\n", - "> v_g = Intercept + u_g\n", - "> = (Intercept + mu_u) + eps_g * sigma_u\n", - "> ```\n", - ">\n", - "> The data only sees the sum `Intercept + mu_u`, so the likelihood is invariant under shifting mass between the two. Their joint posterior has a ridge along the anti-diagonal: samples of `v_Intercept` and `v_1|participant_id_mu` become anticorrelated, the ESS of both drops, and only their sum is well-identified. This is the textbook reason hierarchical models conventionally use **mean-zero random effects**.\n", - ">\n", - "> The practical rule that follows: do **not** estimate the group-mean `mu` from a hyperprior *simultaneously with* fixed effects — the two trade off along exactly this ridge. Keeping a hyperprior on the group **sigma** is fine, and is precisely what lets you recover individual-participant parameters.\n", - ">\n", - "> HSSM detects this layout automatically and emits a second warning: \"User prior for '1|participant_id' on parameter 'v' has a non-trivial `mu`, and the formula also includes a common `Intercept`...\" Fix 2 below is the clean fix — it sets `mu=0` on the group term so the common `Intercept` owns the location and nothing is redundant." + "centered_ridge_graph = pm.model_to_graphviz(centered_ridge_model.pymc_model)\n", + "centered_ridge_graph" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "Vxnm", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "### Fix 2 — move the location prior to the common intercept\n", + "## 5. A valid group-owned location\n", + "\n", + "A free group mean is appropriate when the group term is the unique owner of\n", + "that population location. Remove the matching common intercept:\n", + "\n", + "```python\n", + "v ~ 0 + (1 | participant_id)\n", + "```\n", "\n", - "If you want to keep the non-centered parameterization (usually a good idea for sampling), express the location of the group mean through the *common* `Intercept` term, which is shared across the model regardless of parameterization." + "The outer group prior below carries `noncentered=False`. This per-prior\n", + "setting overrides the model-level `noncentered=True`, retaining the\n", + "requested $\\mu_u$ without changing the parameterization of unrelated\n", + "components." ] }, { "cell_type": "code", "execution_count": null, + "id": "DnEU", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "[]" + "text/html": [ + "
model defaultgroup-prior overridepopulation-location ownerdisconnected free RVs
0non-centeredcentered1|participant_id munone
" ] }, - "execution_count": 10, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_fix_intercept = hssm.HSSM(\n", - " data=cav_data,\n", - " model=\"ddm\",\n", - " include=[\n", + "unique_owner_model, unique_owner_messages = capture_hssm_build(\n", + " lambda: hssm_model(\n", + " formula=\"v ~ 0 + (1|participant_id)\",\n", + " priors={\"1|participant_id\": free_location_group_prior(noncentered=False)},\n", + " noncentered=True,\n", + " )\n", + ")\n", + "_free_names = {rv.name for rv in unique_owner_model.pymc_model.free_RVs}\n", + "assert unique_owner_messages == ()\n", + "assert \"v_Intercept\" not in unique_owner_model.pymc_model.named_vars\n", + "assert \"v_1|participant_id_mu\" in _free_names\n", + "assert \"v_1|participant_id\" in _free_names\n", + "assert \"v_1|participant_id_offset\" not in _free_names\n", + "assert disconnected_free_rvs(unique_owner_model.pymc_model) == []\n", + "\n", + "unique_owner_summary = pd.DataFrame(\n", + " [\n", " {\n", - " \"name\": \"v\",\n", - " \"formula\": \"v ~ 1 + (1|participant_id)\",\n", - " \"prior\": {\n", - " # Location goes here, where both parameterizations use it:\n", - " \"Intercept\": {\"name\": \"Normal\", \"mu\": 0.0, \"sigma\": 0.5},\n", - " \"1|participant_id\": {\n", - " \"name\": \"Normal\",\n", - " \"mu\": 0.0, # scalar, not a hyperprior\n", - " \"sigma\": {\"name\": \"HalfNormal\", \"sigma\": 0.5},\n", - " },\n", - " },\n", + " \"model default\": \"non-centered\",\n", + " \"group-prior override\": \"centered\",\n", + " \"population-location owner\": \"1|participant_id mu\",\n", + " \"disconnected free RVs\": \"none\",\n", " }\n", - " ],\n", - " p_outlier=0.0,\n", - " noncentered=True,\n", + " ]\n", ")\n", - "\n", - "find_disconnected_free_rvs(model_fix_intercept.pymc_model)" + "unique_owner_summary" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, + "id": "ulZA", "metadata": {}, "outputs": [ { "data": { - "text/plain": [ - "Hierarchical Sequential Sampling Model\n", - "Model: ddm\n", - "\n", - "Response variable: rt,response\n", - "Likelihood: analytical\n", - "Observations: 3988\n", - "\n", - "Parameters:\n", - "\n", - "v:\n", - " Formula: v ~ 1 + (1|participant_id)\n", - " Priors:\n", - " v_Intercept ~ Normal(mu: 0.0, sigma: 0.5)\n", - " v_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 0.5))\n", - " Link: identity\n", - " Explicit bounds: (-inf, inf)\n", - "\n", - "a:\n", - " Prior: HalfNormal(sigma: 2.0)\n", - " Explicit bounds: (0.0, inf)\n", - "\n", - "z:\n", - " Prior: Uniform(lower: 0.0, upper: 1.0)\n", - " Explicit bounds: (0.0, 1.0)\n", - "\n", - "t:\n", - " Prior: HalfNormal(sigma: 2.0)\n", - " Explicit bounds: (0.0, inf)" + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (16)\n", + "\n", + "__obs__ (16)\n", + "\n", + "\n", + "cluster__obs__ (16) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (16) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "z\n", + "\n", + "z\n", + "~\n", + "Uniform\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "z->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id_mu\n", + "\n", + "v_1|participant_id_mu\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_1|participant_id\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_1|participant_id_mu->v_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "t\n", + "\n", + "t\n", + "~\n", + "Halfnormal\n", + "\n", + "\n", + "\n", + "t->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma->v_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "a\n", + "\n", + "a\n", + "~\n", + "Halfnormal\n", + "\n", + "\n", + "\n", + "a->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_1|participant_id->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" ] }, - "execution_count": 14, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_fix_intercept" + "unique_owner_graph = pm.model_to_graphviz(unique_owner_model.pymc_model)\n", + "unique_owner_graph" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "ecfG", + "metadata": { + "marimo": { + "md_prefix": "" + } + }, "source": [ - "## 6. Rule of thumb\n", - "\n", - "* Default to `noncentered=True` for hierarchical models — it usually samples better, and is also bambi's default.\n", - "* Under `noncentered=True`, **the location of a group-specific term lives on the common intercept (or other common terms)**, not on the `mu` argument of the group prior. Set `mu=0` on the group term itself.\n", - "* **Even under `noncentered=False`, set `mu=0` on the group term whenever the formula has a common `Intercept`.** Letting both carry a free location creates a non-identifiable ridge — the data only sees their sum. The two warnings HSSM emits (disconnected `_mu` under non-centered, and location ridge under either parameterization) are saying the same thing from two different angles: mean-zero random effects, full stop.\n", - "* Switch to `noncentered=False` only when you have a specific reason (e.g. you want a hyperprior on the group mean *and* you drop the common intercept so the model stays identifiable).\n", - "* When in doubt, plot the PyMC graph with `pm.model_to_graphviz(model.pymc_model)` and look for floating nodes. HSSM also flags them automatically and prints a warning.\n", - "\n", - "> **Coming from HDDM (\"0 +\" models)?** HDDM's convention is the *centered* specification: drop the common intercept (`v ~ 0 + ...`) and let the subject-specific coefficients come from one hierarchical prior with free `mu` and `sigma`. To reproduce it in HSSM, set `noncentered=False`; the group distributions then appear as `v_1|participant_id_mu`-style nodes. With no common intercept in the formula, a free `mu` on the group term is exactly what keeps the model identifiable — this is the rule-of-thumb exception above in action.\n" + "The group `mu` now has a path through the direct group coefficients to the\n", + "likelihood, and there is no competing `v_Intercept`. HSSM's generated safe\n", + "priors use this centered fallback automatically when one unmatched group\n", + "term is the unique population-location owner." ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "Pvdt", + "metadata": { + "marimo": { + "md_prefix": "r" + } + }, "source": [ - "## 7. Reference\n", + "## 6. Choose one owner for every population location\n", "\n", - "The warnings come from the `\"hssm\"` logger. If you want to silence them (not recommended), set its level to `ERROR`:\n", + "The rule is structural and applies to intercepts and slopes alike:\n", "\n", - "```python\n", - "import logging\n", - "logging.getLogger(\"hssm\").setLevel(logging.ERROR)\n", - "```\n", + "- `x + (0 + x | participant_id)`: common `x` owns the population slope;\n", + " the participant coefficients are zero-mean deviations.\n", + "- `(0 + x | participant_id)` with no common `x`: that single group term may\n", + " own the population slope, so a free location must be effectively centered.\n", + "- `(0 + x | participant_id) + (0 + x | item_id)` with no common `x`: two\n", + " group means compete for the same location. Add common `x` and make both\n", + " group terms zero-mean, or deliberately choose exactly one owner.\n", "\n", - "If you want to inspect the graph programmatically, use `hssm.param.parameterization_check.find_disconnected_free_rvs(model.pymc_model)`.\n", + "An inverse link does not change this logic. Common and group coefficients\n", + "are combined on the linear-predictor scale before the inverse link is\n", + "applied." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ZBYS", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
formula structurelocation ownergroup muHSSM result
0exact common/group matchcommon termfixed at zerocentered or non-centered is valid
1one unmatched group termgroup distributionmay be freeuse effective centering
2repeated unmatched group expressionambiguous until specifieddo not free every meansafe generation rejects; explicit ridges warn
3explicit incompatible NC priorwould be discarded by Bambifree or nonzeropre-build ValueError
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "parameterization_rules = pd.DataFrame(\n", + " [\n", + " {\n", + " \"formula structure\": \"exact common/group match\",\n", + " \"location owner\": \"common term\",\n", + " \"group mu\": \"fixed at zero\",\n", + " \"HSSM result\": \"centered or non-centered is valid\",\n", + " },\n", + " {\n", + " \"formula structure\": \"one unmatched group term\",\n", + " \"location owner\": \"group distribution\",\n", + " \"group mu\": \"may be free\",\n", + " \"HSSM result\": \"use effective centering\",\n", + " },\n", + " {\n", + " \"formula structure\": \"repeated unmatched group expression\",\n", + " \"location owner\": \"ambiguous until specified\",\n", + " \"group mu\": \"do not free every mean\",\n", + " \"HSSM result\": \"safe generation rejects; explicit ridges warn\",\n", + " },\n", + " {\n", + " \"formula structure\": \"explicit incompatible NC prior\",\n", + " \"location owner\": \"would be discarded by Bambi\",\n", + " \"group mu\": \"free or nonzero\",\n", + " \"HSSM result\": \"pre-build ValueError\",\n", + " },\n", + " ]\n", + ")\n", + "parameterization_rules" + ] + }, + { + "cell_type": "markdown", + "id": "aLJB", + "metadata": { + "marimo": { + "md_prefix": "" + } + }, + "source": [ + "## 7. Where to continue\n", + "\n", + "- [Specify hierarchical group priors](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/)\n", + " gives the complete explicit-prior compatibility contract.\n", + "- [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/)\n", + " explains why population locations live on the linear-predictor scale.\n", + "- [Choosing a parameterization per parameter](https://lnccbrown.github.io/HSSM/tutorials/parameterization_per_parameter/)\n", + " shows model, component, and per-prior override precedence.\n", + "- Betancourt's [hierarchical modeling case study](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html)\n", + " develops the posterior geometry behind centered and non-centered choices.\n", "\n", - "**Further reading:** Betancourt's [hierarchical modeling case study](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html) covers the posterior geometry behind these choices in depth. For choosing the parameterization per parameter (the two are mathematically equivalent — in practice you can try both), see [Per-parameter centered vs. non-centered parameterization](https://lnccbrown.github.io/HSSM/tutorials/parameterization_per_parameter/).\n" + "The practical takeaway is simple: parameterization may change coordinates,\n", + "but it must not change who owns a population location. HSSM's preflight\n", + "checks enforce that boundary before Bambi constructs a different graph." ] } ], "metadata": { - "kernelspec": { - "display_name": ".venv", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -1339,10 +1614,16 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.13.2" + "pygments_lexer": "ipython3" + }, + "marimo": { + "app_config": { + "width": "medium" + }, + "header": "# /// script\n# requires-python = \">=3.12,<3.15\"\n# dependencies = [\n# \"bambi==0.20.0\",\n# \"graphviz==0.21\",\n# \"hssm @ git+https://github.com/lnccbrown/HSSM.git@b6a6bdcf68ecd7cf71ffdef5cdda4fb05e8bfaad\",\n# \"marimo==0.24.0\",\n# \"matplotlib==3.11.1\",\n# \"numpy==2.4.6\",\n# \"pandas==3.0.5\",\n# \"pymc==6.3.1\",\n# ]\n# ///\n\n\"\"\"Explain centered and non-centered group effects in current HSSM.\n\nThis construction-only marimo tutorial separates the mathematical\nreparameterization from Bambi's current implementation, demonstrates HSSM's\npre-build compatibility guard, and compares the model graphs for valid and\nproblematic group-location layouts. No sampling is required.\n\nRun the pinned standalone environment locally or in Molab::\n\n uvx marimo edit --sandbox \\\n docs/tutorials/centered_vs_noncentered_basic_logic.py\n\nTo exercise an active HSSM checkout instead, ignore the inline environment::\n\n uv run --group notebook --group docs marimo edit --no-sandbox \\\n docs/tutorials/centered_vs_noncentered_basic_logic.py\n uv run --group notebook --group docs marimo check --strict \\\n docs/tutorials/centered_vs_noncentered_basic_logic.py\n uv run --group notebook --group docs marimo export html --no-sandbox \\\n docs/tutorials/centered_vs_noncentered_basic_logic.py \\\n --output /tmp/centered-vs-noncentered.html --force\n uv run --group notebook --group docs marimo export ipynb --no-sandbox \\\n docs/tutorials/centered_vs_noncentered_basic_logic.py \\\n --output docs/tutorials/centered_vs_noncentered_basic_logic.ipynb \\\n --include-outputs --force\n uv run ruff format \\\n docs/tutorials/centered_vs_noncentered_basic_logic.ipynb\n\"\"\"\n\n# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns)\n", + "marimo_version": "0.24.0" } }, "nbformat": 4, - "nbformat_minor": 4 + "nbformat_minor": 5 } diff --git a/docs/tutorials/centered_vs_noncentered_basic_logic.py b/docs/tutorials/centered_vs_noncentered_basic_logic.py new file mode 100644 index 000000000..c05ac7949 --- /dev/null +++ b/docs/tutorials/centered_vs_noncentered_basic_logic.py @@ -0,0 +1,759 @@ +# /// script +# requires-python = ">=3.12,<3.15" +# dependencies = [ +# "bambi==0.20.0", +# "graphviz==0.21", +# "hssm @ git+https://github.com/lnccbrown/HSSM.git@b6a6bdcf68ecd7cf71ffdef5cdda4fb05e8bfaad", +# "marimo==0.24.0", +# "matplotlib==3.11.1", +# "numpy==2.4.6", +# "pandas==3.0.5", +# "pymc==6.3.1", +# ] +# /// + +"""Explain centered and non-centered group effects in current HSSM. + +This construction-only marimo tutorial separates the mathematical +reparameterization from Bambi's current implementation, demonstrates HSSM's +pre-build compatibility guard, and compares the model graphs for valid and +problematic group-location layouts. No sampling is required. + +Run the pinned standalone environment locally or in Molab:: + + uvx marimo edit --sandbox \ + docs/tutorials/centered_vs_noncentered_basic_logic.py + +To exercise an active HSSM checkout instead, ignore the inline environment:: + + uv run --group notebook --group docs marimo edit --no-sandbox \ + docs/tutorials/centered_vs_noncentered_basic_logic.py + uv run --group notebook --group docs marimo check --strict \ + docs/tutorials/centered_vs_noncentered_basic_logic.py + uv run --group notebook --group docs marimo export html --no-sandbox \ + docs/tutorials/centered_vs_noncentered_basic_logic.py \ + --output /tmp/centered-vs-noncentered.html --force + uv run --group notebook --group docs marimo export ipynb --no-sandbox \ + docs/tutorials/centered_vs_noncentered_basic_logic.py \ + --output docs/tutorials/centered_vs_noncentered_basic_logic.ipynb \ + --include-outputs --force + uv run ruff format \ + docs/tutorials/centered_vs_noncentered_basic_logic.ipynb +""" + +# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns) +import marimo + +__generated_with = "0.24.0" +app = marimo.App(width="medium") + + +@app.cell +def _(): + import logging + import os + import warnings + from tempfile import gettempdir + + # This notebook only constructs models. Molab can expose a CUDA plugin even + # when no usable GPU is present, so select CPU before importing HSSM/JAX. + os.environ["JAX_PLATFORMS"] = "cpu" + os.environ["JAX_SKIP_CUDA_CONSTRAINTS_CHECK"] = "1" + os.environ["XLA_PYTHON_CLIENT_PREALLOCATE"] = "false" + os.environ.setdefault( + "MPLCONFIGDIR", f"{gettempdir()}/hssm-parameterization-matplotlib" + ) + + warnings.filterwarnings("ignore") + logging.getLogger("jax._src.xla_bridge").setLevel(logging.CRITICAL) + logging.getLogger("matplotlib").setLevel(logging.ERROR) + + import bambi as bmb + import marimo as mo + import numpy as np + import pandas as pd + import pymc as pm + from pytensor.graph.traversal import ancestors + + import hssm + + logging.getLogger("hssm").setLevel(logging.WARNING) + hssm.set_floatX("float64") + pd.set_option("display.max_colwidth", 100) + return ancestors, bmb, hssm, logging, mo, np, pd, pm + + +@app.cell +def _(bmb, hssm, mo): + mo.md(f""" + # Centered vs. non-centered parameterizations + + Centering is a choice about **how a hierarchical effect is represented for + computation**. It should not change the statistical model. That distinction + becomes important when a group-specific prior has its own population mean. + + This tutorial uses **HSSM {hssm.__version__}** and **Bambi {bmb.__version__}**. + It constructs models and inspects their PyMC graphs; no MCMC is run. + + By the end, you should be able to: + + 1. distinguish the mathematical non-centered transformation from Bambi's + current shortcut; + 2. recognize a valid zero-mean group deviation in either parameterization; + 3. understand why HSSM rejects an explicit non-centered free group mean; and + 4. choose exactly one owner for each population location. + """) + return + + +@app.cell +def _(mo): + mo.md(r""" + ## 1. Same distribution, different coordinates + + Let $u_g$ be a coefficient for group $g$, with population location + $\mu$ and group scale $\sigma$. + + In the **centered** parameterization we sample the coefficient directly: + + $$ + u_g \sim \mathcal{N}(\mu, \sigma). + $$ + + A mathematically equivalent **non-centered** parameterization samples a + standard-normal coordinate and transforms it: + + $$ + z_g \sim \mathcal{N}(0, 1), + \qquad + u_g = \mu + \sigma z_g. + $$ + + These equations define the same prior distribution for $u_g$. Which + coordinates sample better depends on the amount of information in the data + and the posterior geometry; neither form is universally superior. + """) + return + + +@app.cell +def _(mo): + mo.md(r""" + > **The current Bambi boundary** + + For a non-centered group term, Bambi currently constructs + + $$ + u_g = \sigma z_g, + $$ + + rather than $\mu + \sigma z_g$. This shortcut is faithful when the group + term is a zero-mean deviation: `mu` is absent or fixed entirely to zero. A + free or nonzero `mu` would be dropped. + + HSSM validates explicit group priors before asking Bambi to build the PyMC + model. If that shortcut would discard part of the requested prior, HSSM + raises a `ValueError` instead of silently changing the model. + """) + return + + +@app.cell +def _(ancestors, bmb, hssm, logging, np, pd): + tutorial_data = pd.DataFrame( + { + "rt": 0.38 + 0.015 * np.arange(16), + "response": np.where(np.arange(16) % 2, 1, -1), + "theta": np.linspace(-1.0, 1.0, 16), + "participant_id": np.repeat(np.arange(4), 4), + } + ) + raw_bambi_data = tutorial_data.loc[:, ["participant_id"]].assign( + y=np.array( + [ + -0.20, + 0.05, + 0.10, + -0.05, + 0.30, + 0.35, + 0.20, + 0.40, + -0.35, + -0.20, + -0.10, + -0.25, + 0.15, + 0.25, + 0.05, + 0.30, + ] + ) + ) + + def zero_mean_group_prior(*, noncentered=None): + """Return a fresh zero-mean hierarchical Normal group prior.""" + return hssm.Prior( + "Normal", + mu=0.0, + sigma=hssm.Prior("HalfNormal", sigma=0.5), + noncentered=noncentered, + ) + + def free_location_group_prior(*, noncentered=None): + """Return a fresh hierarchical Normal with a free population mean.""" + return hssm.Prior( + "Normal", + mu=hssm.Prior("Normal", mu=0.0, sigma=0.5), + sigma=hssm.Prior("HalfNormal", sigma=0.5), + noncentered=noncentered, + ) + + def matched_include(group_prior): + """Place a common intercept and matching group intercept in one spec.""" + return [ + { + "name": "v", + "formula": "v ~ 1 + (1|participant_id)", + "prior": { + "Intercept": hssm.Prior("Normal", mu=0.0, sigma=0.5), + "1|participant_id": group_prior, + }, + } + ] + + def hssm_model(*, formula, priors, noncentered=True): + """Build a tiny analytical DDM without sampling or init-value work.""" + return hssm.HSSM( + data=tutorial_data, + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": formula, "prior": priors}], + p_outlier=0.0, + prior_settings=None, + noncentered=noncentered, + process_initvals=False, + initval_jitter=0.0, + ) + + def disconnected_free_rvs(pymc_model): + """Return free-RV names that are not ancestors of observed variables.""" + connected = { + id(variable) + for observed_rv in pymc_model.observed_RVs + for variable in ancestors([observed_rv]) + } + return sorted(rv.name for rv in pymc_model.free_RVs if id(rv) not in connected) + + def capture_hssm_build(builder): + """Build a model while collecting HSSM warning messages.""" + logger = logging.getLogger("hssm") + messages = [] + + class _MessageHandler(logging.Handler): + def emit(self, record): + messages.append(record.getMessage()) + + handler = _MessageHandler(level=logging.WARNING) + previous_handlers = list(logger.handlers) + previous_level = logger.level + previous_propagate = logger.propagate + logger.handlers = [handler] + logger.setLevel(logging.WARNING) + logger.propagate = False + try: + result = builder() + finally: + logger.handlers = previous_handlers + logger.setLevel(previous_level) + logger.propagate = previous_propagate + return result, tuple(messages) + + # Assert the factories return fresh prior trees. HSSM/Bambi attach names + # while preparing priors, so tutorial cases must not share mutable objects. + assert zero_mean_group_prior() is not zero_mean_group_prior() + assert free_location_group_prior() is not free_location_group_prior() + assert isinstance(free_location_group_prior().args["mu"], bmb.Prior) + return ( + capture_hssm_build, + disconnected_free_rvs, + free_location_group_prior, + hssm_model, + matched_include, + raw_bambi_data, + tutorial_data, + zero_mean_group_prior, + ) + + +@app.cell +def _(mo): + mo.md(r""" + ## 2. A valid zero-mean group deviation + + Consider + + ```python + v ~ 1 + (1 | participant_id) + ``` + + The common `Intercept` owns the population location. The participant term + is a deviation around that location, so its prior has `mu=0` and a + hierarchical `sigma`. This statistical model can be represented faithfully + in either parameterization. + """) + return + + +@app.cell +def _( + capture_hssm_build, + disconnected_free_rvs, + hssm, + matched_include, + pd, + tutorial_data, + zero_mean_group_prior, +): + matched_centered_model, matched_centered_messages = capture_hssm_build( + lambda: hssm.HSSM( + data=tutorial_data, + model="ddm", + loglik_kind="analytical", + include=matched_include(zero_mean_group_prior()), + p_outlier=0.0, + prior_settings=None, + noncentered=False, + process_initvals=False, + initval_jitter=0.0, + ) + ) + matched_noncentered_model, matched_noncentered_messages = capture_hssm_build( + lambda: hssm.HSSM( + data=tutorial_data, + model="ddm", + loglik_kind="analytical", + include=matched_include(zero_mean_group_prior()), + p_outlier=0.0, + prior_settings=None, + noncentered=True, + process_initvals=False, + initval_jitter=0.0, + ) + ) + + _centered_names = {rv.name for rv in matched_centered_model.pymc_model.free_RVs} + _noncentered_names = { + rv.name for rv in matched_noncentered_model.pymc_model.free_RVs + } + assert matched_centered_messages == () + assert matched_noncentered_messages == () + assert "v_1|participant_id" in _centered_names + assert "v_1|participant_id_offset" not in _centered_names + assert "v_1|participant_id_offset" in _noncentered_names + assert "v_1|participant_id" not in _noncentered_names + assert "v_1|participant_id_mu" not in _centered_names | _noncentered_names + assert disconnected_free_rvs(matched_centered_model.pymc_model) == [] + assert disconnected_free_rvs(matched_noncentered_model.pymc_model) == [] + + matched_parameterization_table = pd.DataFrame( + [ + { + "effective parameterization": "centered", + "sampled group coordinate": "v_1|participant_id", + "group coefficient in graph": "free RV", + "disconnected free RVs": "none", + }, + { + "effective parameterization": "non-centered", + "sampled group coordinate": "v_1|participant_id_offset", + "group coefficient in graph": "deterministic", + "disconnected free RVs": "none", + }, + ] + ) + matched_parameterization_table + return ( + matched_centered_model, + matched_noncentered_model, + matched_parameterization_table, + ) + + +@app.cell +def _(matched_centered_model, pm): + matched_centered_graph = pm.model_to_graphviz(matched_centered_model.pymc_model) + matched_centered_graph + return (matched_centered_graph,) + + +@app.cell +def _(mo): + mo.md(""" + In the centered graph, `v_1|participant_id` is sampled directly. Its scale + hyperprior points to the group coefficients, and those coefficients point + to the trial-wise drift rate. + """) + return + + +@app.cell +def _(matched_noncentered_model, pm): + matched_noncentered_graph = pm.model_to_graphviz( + matched_noncentered_model.pymc_model + ) + matched_noncentered_graph + return (matched_noncentered_graph,) + + +@app.cell +def _(mo): + mo.md(""" + In the non-centered graph, Bambi samples a standard-normal `offset` and + combines it with `sigma`. The resulting `v_1|participant_id` is a + deterministic node. Because the intended group mean is exactly zero, this + graph represents the same prior as the centered graph above. + """) + return + + +@app.cell +def _(mo): + mo.md(r""" + ## 3. Why HSSM needs a pre-build guard + + To isolate the underlying mechanism, the next example uses Bambi directly + with a tiny Gaussian response. The requested group prior has a free + population location: + + ```python + Normal( + mu=Normal(0, 0.5), + sigma=HalfNormal(0.5), + ) + ``` + + This is a legitimate centered hierarchy. Under Bambi's current + non-centered shortcut, however, the `mu` prior is created and then omitted + from `offset * sigma`. + """) + return + + +@app.cell +def _(bmb, disconnected_free_rvs, pd, raw_bambi_data): + raw_bambi_group_prior = bmb.Prior( + "Normal", + mu=bmb.Prior("Normal", mu=0.0, sigma=0.5), + sigma=bmb.Prior("HalfNormal", sigma=0.5), + ) + raw_bambi_model = bmb.Model( + "y ~ 1 + (1|participant_id)", + raw_bambi_data, + family="gaussian", + priors={"1|participant_id": raw_bambi_group_prior}, + noncentered=True, + ) + raw_bambi_model.build() + raw_bambi_pymc_model = raw_bambi_model.backend.model + raw_bambi_orphans = disconnected_free_rvs(raw_bambi_pymc_model) + + _free_names = [rv.name for rv in raw_bambi_pymc_model.free_RVs] + assert "1|participant_id_mu" in _free_names + assert "1|participant_id_offset" in _free_names + assert raw_bambi_orphans == ["1|participant_id_mu"] + + raw_bambi_summary = pd.DataFrame( + { + "quantity": ["sampled group coordinates", "disconnected free RVs"], + "value": [", ".join(_free_names), ", ".join(raw_bambi_orphans)], + } + ) + raw_bambi_summary + return raw_bambi_model, raw_bambi_orphans, raw_bambi_pymc_model, raw_bambi_summary + + +@app.cell +def _(pm, raw_bambi_pymc_model): + raw_bambi_orphan_graph = pm.model_to_graphviz(raw_bambi_pymc_model) + raw_bambi_orphan_graph + return (raw_bambi_orphan_graph,) + + +@app.cell +def _(mo): + mo.md(""" + > The floating `1|participant_id_mu` node has no path to the observed + > response. Sampling it would consume computation without changing the + > likelihood. This graph documents the **raw Bambi behavior that HSSM + > prevents**; it is not a graph that current HSSM will construct from the + > same explicit prior. + """) + return + + +@app.cell +def _(free_location_group_prior, hssm, matched_include, mo, tutorial_data): + try: + hssm.HSSM( + data=tutorial_data, + model="ddm", + loglik_kind="analytical", + include=matched_include(free_location_group_prior()), + p_outlier=0.0, + prior_settings=None, + noncentered=True, + process_initvals=False, + initval_jitter=0.0, + ) + except ValueError as exc: + hssm_preflight_error = str(exc) + else: + raise AssertionError("HSSM did not reject the incompatible group prior") + + assert "cannot be represented faithfully by bambi" in hssm_preflight_error + assert "1|participant_id" in hssm_preflight_error + assert "omits from `offset * sigma`" in hssm_preflight_error + assert "noncentered=False" in hssm_preflight_error + mo.md(f""" + **Current HSSM result** + + ```text + {hssm_preflight_error} + ``` + """) + return (hssm_preflight_error,) + + +@app.cell +def _(mo): + mo.md(r""" + ## 4. Centering retains the mean—but location ownership still matters + + If we center the explicit hierarchy, Bambi uses its `mu`. With both a + common intercept and a free group mean, the predictor for participant $g$ + contains + + $$ + \eta_g = \beta_0 + u_g, + \qquad + u_g \sim \mathcal{N}(\mu_u, \sigma_u). + $$ + + The likelihood sees $\beta_0 + \mu_u$, not the two locations separately. + Shifting one up and the other down leaves the predictor unchanged. The PyMC + graph is fully connected, but the likelihood has a ridge along that shift + direction. + """) + return + + +@app.cell +def _( + capture_hssm_build, + disconnected_free_rvs, + free_location_group_prior, + hssm, + matched_include, + mo, + tutorial_data, +): + centered_ridge_model, centered_ridge_messages = capture_hssm_build( + lambda: hssm.HSSM( + data=tutorial_data, + model="ddm", + loglik_kind="analytical", + include=matched_include(free_location_group_prior()), + p_outlier=0.0, + prior_settings=None, + noncentered=False, + process_initvals=False, + initval_jitter=0.0, + ) + ) + centered_ridge_warning = next( + message + for message in centered_ridge_messages + if "non-identifiable individually" in message + ) + + assert "posterior will have a ridge" in centered_ridge_warning + assert "common 'Intercept' effect" in centered_ridge_warning + assert disconnected_free_rvs(centered_ridge_model.pymc_model) == [] + assert "v_1|participant_id_mu" in { + rv.name for rv in centered_ridge_model.pymc_model.free_RVs + } + mo.md(f""" + **HSSM's centered location warning** + + ```text + {centered_ridge_warning} + ``` + + The disconnected-node check returns `[]`: this is an identifiability + problem in the likelihood, not an orphan-node problem. + """) + return centered_ridge_model, centered_ridge_warning + + +@app.cell +def _(centered_ridge_model, pm): + centered_ridge_graph = pm.model_to_graphviz(centered_ridge_model.pymc_model) + centered_ridge_graph + return (centered_ridge_graph,) + + +@app.cell +def _(mo): + mo.md(r""" + ## 5. A valid group-owned location + + A free group mean is appropriate when the group term is the unique owner of + that population location. Remove the matching common intercept: + + ```python + v ~ 0 + (1 | participant_id) + ``` + + The outer group prior below carries `noncentered=False`. This per-prior + setting overrides the model-level `noncentered=True`, retaining the + requested $\mu_u$ without changing the parameterization of unrelated + components. + """) + return + + +@app.cell +def _( + capture_hssm_build, + disconnected_free_rvs, + free_location_group_prior, + hssm_model, + pd, +): + unique_owner_model, unique_owner_messages = capture_hssm_build( + lambda: hssm_model( + formula="v ~ 0 + (1|participant_id)", + priors={"1|participant_id": free_location_group_prior(noncentered=False)}, + noncentered=True, + ) + ) + _free_names = {rv.name for rv in unique_owner_model.pymc_model.free_RVs} + assert unique_owner_messages == () + assert "v_Intercept" not in unique_owner_model.pymc_model.named_vars + assert "v_1|participant_id_mu" in _free_names + assert "v_1|participant_id" in _free_names + assert "v_1|participant_id_offset" not in _free_names + assert disconnected_free_rvs(unique_owner_model.pymc_model) == [] + + unique_owner_summary = pd.DataFrame( + [ + { + "model default": "non-centered", + "group-prior override": "centered", + "population-location owner": "1|participant_id mu", + "disconnected free RVs": "none", + } + ] + ) + unique_owner_summary + return unique_owner_model, unique_owner_summary + + +@app.cell +def _(pm, unique_owner_model): + unique_owner_graph = pm.model_to_graphviz(unique_owner_model.pymc_model) + unique_owner_graph + return (unique_owner_graph,) + + +@app.cell +def _(mo): + mo.md(""" + The group `mu` now has a path through the direct group coefficients to the + likelihood, and there is no competing `v_Intercept`. HSSM's generated safe + priors use this centered fallback automatically when one unmatched group + term is the unique population-location owner. + """) + return + + +@app.cell +def _(mo): + mo.md(r""" + ## 6. Choose one owner for every population location + + The rule is structural and applies to intercepts and slopes alike: + + - `x + (0 + x | participant_id)`: common `x` owns the population slope; + the participant coefficients are zero-mean deviations. + - `(0 + x | participant_id)` with no common `x`: that single group term may + own the population slope, so a free location must be effectively centered. + - `(0 + x | participant_id) + (0 + x | item_id)` with no common `x`: two + group means compete for the same location. Add common `x` and make both + group terms zero-mean, or deliberately choose exactly one owner. + + An inverse link does not change this logic. Common and group coefficients + are combined on the linear-predictor scale before the inverse link is + applied. + """) + return + + +@app.cell +def _(pd): + parameterization_rules = pd.DataFrame( + [ + { + "formula structure": "exact common/group match", + "location owner": "common term", + "group mu": "fixed at zero", + "HSSM result": "centered or non-centered is valid", + }, + { + "formula structure": "one unmatched group term", + "location owner": "group distribution", + "group mu": "may be free", + "HSSM result": "use effective centering", + }, + { + "formula structure": "repeated unmatched group expression", + "location owner": "ambiguous until specified", + "group mu": "do not free every mean", + "HSSM result": "safe generation rejects; explicit ridges warn", + }, + { + "formula structure": "explicit incompatible NC prior", + "location owner": "would be discarded by Bambi", + "group mu": "free or nonzero", + "HSSM result": "pre-build ValueError", + }, + ] + ) + parameterization_rules + return (parameterization_rules,) + + +@app.cell +def _(mo): + mo.md(""" + ## 7. Where to continue + + - [Specify hierarchical group priors](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) + gives the complete explicit-prior compatibility contract. + - [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/) + explains why population locations live on the linear-predictor scale. + - [Choosing a parameterization per parameter](https://lnccbrown.github.io/HSSM/tutorials/parameterization_per_parameter/) + shows model, component, and per-prior override precedence. + - Betancourt's [hierarchical modeling case study](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html) + develops the posterior geometry behind centered and non-centered choices. + + The practical takeaway is simple: parameterization may change coordinates, + but it must not change who owns a population location. HSSM's preflight + checks enforce that boundary before Bambi constructs a different graph. + """) + return + + +if __name__ == "__main__": + app.run() From b68309aa5a2986515f05f87908a61fbf8e46ce2e Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 17:11:20 -0400 Subject: [PATCH 2/3] docs: migrate per-parameterization tutorial (#1273) --- .../parameterization_per_parameter.ipynb | 1945 ++++++++++++++--- .../parameterization_per_parameter.py | 883 ++++++++ 2 files changed, 2493 insertions(+), 335 deletions(-) create mode 100644 docs/tutorials/parameterization_per_parameter.py diff --git a/docs/tutorials/parameterization_per_parameter.ipynb b/docs/tutorials/parameterization_per_parameter.ipynb index 60fa77d0c..69bc0c8ce 100644 --- a/docs/tutorials/parameterization_per_parameter.ipynb +++ b/docs/tutorials/parameterization_per_parameter.ipynb @@ -1,490 +1,1759 @@ { "cells": [ { - "cell_type": "markdown", - "id": "e93eb83e", + "cell_type": "code", + "execution_count": null, + "id": "Hbol", "metadata": {}, + "outputs": [], "source": [ - "# Choosing a parameterization per parameter\n", - "\n", - "Hierarchical models can express a group-specific (random) effect in one of two\n", - "equivalent ways:\n", + "import logging\n", + "import os\n", + "import warnings\n", + "from contextlib import redirect_stderr, redirect_stdout\n", + "from io import StringIO\n", + "from tempfile import gettempdir\n", "\n", - "- **Centered:** $\\theta_g \\sim \\mathrm{Normal}(\\mu, \\sigma)$.\n", - "- **Non-centered:** $\\theta_g = \\mu + \\sigma \\, z_g$, with $z_g \\sim \\mathrm{Normal}(0, 1)$.\n", + "# This notebook only constructs graphs. Molab and some development hosts can\n", + "# expose an unusable CUDA plugin, so prevent JAX from probing it.\n", + "os.environ[\"JAX_PLATFORMS\"] = \"cpu\"\n", + "os.environ[\"JAX_SKIP_CUDA_CONSTRAINTS_CHECK\"] = \"1\"\n", + "os.environ.setdefault(\"MPLCONFIGDIR\", f\"{gettempdir()}/hssm-per-parameter-matplotlib\")\n", + "warnings.filterwarnings(\"ignore\")\n", + "logging.getLogger(\"jax._src.xla_bridge\").setLevel(logging.CRITICAL)\n", "\n", - "They imply the same model but sample very differently. Non-centered usually\n", - "samples better when a group is only weakly informed by the data (it avoids\n", - "Neal's funnel); centered can be better when a group is strongly informed. In a\n", - "model with several parameters, the better choice can **differ from parameter to\n", - "parameter** — so a single global switch is often too coarse.\n", + "import bambi as bmb\n", + "import jax\n", + "import marimo as mo\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pymc as pm\n", "\n", - "For the full mechanism — how HSSM names the PyMC nodes under each parameterization, the warnings it emits, and the identifiability footgun — see [Centered vs. non-centered parameterizations](https://lnccbrown.github.io/HSSM/tutorials/centered_vs_noncentered_basic_logic/). The two forms describe the same model **when the group-level `mu` is zero** — HSSM's non-centered form is `offset * sigma`, so a non-zero group `mu` is ignored (put the location on the common `Intercept` instead); given that, what differs is posterior geometry ([Betancourt's case study](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html) treats this in depth), so in practice you can try both — per parameter.\n", + "# HSSM reports backend configuration and registry details during import and\n", + "# setup. Suppress that incidental output; later helpers capture HSSM's\n", + "# parameterization warnings explicitly and display them in the relevant cell.\n", + "with redirect_stdout(StringIO()), redirect_stderr(StringIO()):\n", + " import hssm\n", + " from hssm.param.parameterization_check import find_disconnected_free_rvs\n", "\n", - "HSSM surfaces bambi's per-parameter control, letting you pick the\n", - "parameterization **per parameter**, or even **per term**.\n", + " hssm.set_floatX(\"float64\")\n", "\n", - "> **Requirements.** The per-parameter forms below (a `noncentered` *dict* and the\n", - "> per-prior `noncentered` field) require **bambi ≥ 0.19**\n", - "> ([PR #983](https://github.com/bambinos/bambi/pull/983)), which HSSM depends on.\n", - "> The plain `noncentered=True/False` form works on any supported bambi." + "assert jax.default_backend() == \"cpu\"\n", + "pd.set_option(\"display.max_colwidth\", 100)" ] }, { - "cell_type": "markdown", - "id": "8503894c", - "metadata": {}, + "cell_type": "code", + "execution_count": null, + "id": "MJUe", + "metadata": { + "jupyter": { + "source_hidden": true + }, + "marimo": { + "config": { + "hide_code": true + } + }, + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "text/markdown": [ + "# Choosing a parameterization per parameter\n", + "\n", + "Hierarchical models often sample better when different parameters—or even\n", + "different group terms—use different parameterizations. This tutorial shows\n", + "how HSSM resolves those choices and how to verify the resulting PyMC graph.\n", + "\n", + "By the end, you will be able to:\n", + "\n", + "1. choose a model-wide or component-specific default;\n", + "2. override one explicit group prior safely;\n", + "3. recognize HSSM's pre-build errors and location-ridge warnings; and\n", + "4. decide which formula term owns each population location.\n", + "\n", + "Everything below is structural: the models are built but never sampled.\n", + "\n", + "**Environment:** HSSM `0.4.0`, Bambi `0.20.0`,\n", + "JAX `0.11.1` on `cpu`." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "## Setup" + "mo.md(f\"\"\"\n", + "# Choosing a parameterization per parameter\n", + "\n", + "Hierarchical models often sample better when different parameters—or even\n", + "different group terms—use different parameterizations. This tutorial shows\n", + "how HSSM resolves those choices and how to verify the resulting PyMC graph.\n", + "\n", + "By the end, you will be able to:\n", + "\n", + "1. choose a model-wide or component-specific default;\n", + "2. override one explicit group prior safely;\n", + "3. recognize HSSM's pre-build errors and location-ridge warnings; and\n", + "4. decide which formula term owns each population location.\n", + "\n", + "Everything below is structural: the models are built but never sampled.\n", + "\n", + "**Environment:** HSSM `{hssm.__version__}`, Bambi `{bmb.__version__}`,\n", + "JAX `{jax.__version__}` on `{jax.default_backend()}`.\n", + "\"\"\")" ] }, { - "cell_type": "code", - "execution_count": 1, - "id": "e3da97f5", + "cell_type": "markdown", + "id": "vblA", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:15.011969Z", - "iopub.status.busy": "2026-07-13T00:03:15.011851Z", - "iopub.status.idle": "2026-07-13T00:03:16.789634Z", - "shell.execute_reply": "2026-07-13T00:03:16.789261Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "r" } }, + "source": [ + "## Centering is both geometry and model structure\n", + "\n", + "A textbook hierarchical Normal can be written in centered form,\n", + "\n", + "$$b_g \\sim \\mathcal N(\\mu, \\sigma),$$\n", + "\n", + "or in the mathematically equivalent non-centered form,\n", + "\n", + "$$z_g \\sim \\mathcal N(0,1), \\qquad b_g = \\mu + \\sigma z_g.$$\n", + "\n", + "The likelihood can be identical while the posterior geometry—and therefore\n", + "sampling efficiency—changes. Non-centering is often helpful for weakly\n", + "informed groups; centering can be better for strongly informed groups.\n", + "\n", + "There is one important implementation boundary. Current Bambi constructs a\n", + "non-centered group term as\n", + "\n", + "$$b_g = \\sigma z_g,$$\n", + "\n", + "so this route faithfully represents a built-in Normal group prior only when\n", + "its location is absent or fixed entirely to zero and `sigma` is hierarchical.\n", + "HSSM checks explicit priors before asking Bambi to build the model. A free or\n", + "nonzero `mu` under effective non-centering now raises an actionable error; it\n", + "is not silently accepted and then discarded.\n", + "\n", + "The [hierarchical group-prior guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/)\n", + "gives the complete compatibility table. For the scale on which these effects\n", + "combine, see [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bkHC", + "metadata": {}, "outputs": [ { - "name": "stderr", - "output_type": "stream", - "text": [ - "You supplied a model 'lba4', which is currently not supported in the ssm_simulators package. An error will be thrown when sampling from the random variable or when using any posterior or prior predictive sampling methods.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Setting PyTensor floatX type to float32.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Setting \"jax_enable_x64\" to False. If this is not intended, please set `jax` to False.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "hssm 0.4.0 | bambi 0.19.0\n" - ] + "data": { + "text/html": [ + "
rtresponsethetaparticipant_idconf
00.42-1-1.00low
10.4310.00high
20.44-11.00low
30.451-1.00high
40.46-10.00low
50.4711.00high
60.48-1-1.01low
70.4910.01high
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "import logging\n", - "import warnings\n", + "_trial = np.arange(24)\n", + "tutorial_data = pd.DataFrame(\n", + " {\n", + " \"rt\": 0.42 + 0.01 * _trial,\n", + " \"response\": np.where(_trial % 2, 1, -1),\n", + " \"theta\": np.tile([-1.0, 0.0, 1.0], 8),\n", + " \"participant_id\": np.repeat(np.arange(4), 6),\n", + " \"conf\": np.tile([\"low\", \"high\"], 12),\n", + " }\n", + ")\n", + "tutorial_data.head(8)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "lEQa", + "metadata": {}, + "outputs": [], + "source": [ + "_base_model_kwargs = {\n", + " \"data\": tutorial_data,\n", + " \"model\": \"ddm\",\n", + " \"loglik_kind\": \"analytical\",\n", + " \"p_outlier\": 0.0,\n", + " \"prior_settings\": \"safe\",\n", + " \"process_initvals\": False,\n", + " \"initval_jitter\": 0.0,\n", + " \"z\": 0.5,\n", + " \"t\": 0.2,\n", + "}\n", "\n", - "warnings.simplefilter(\"ignore\")\n", - "logging.getLogger(\"hssm\").setLevel(logging.ERROR) # quiet info/warnings for the demo\n", "\n", - "import bambi as bmb\n", + "def build_model(include, **kwargs):\n", + " \"\"\"Build quietly while returning only HSSM warning messages.\"\"\"\n", + " _log_stream = StringIO()\n", + " _handler = logging.StreamHandler(_log_stream)\n", + " _handler.setFormatter(logging.Formatter(\"%(message)s\"))\n", + " _logger = logging.getLogger(\"hssm\")\n", + " _old_handlers = list(_logger.handlers)\n", + " _old_level = _logger.level\n", + " _old_propagate = _logger.propagate\n", + " _logger.handlers = [_handler]\n", + " _logger.setLevel(logging.WARNING)\n", + " _logger.propagate = False\n", + " try:\n", + " with redirect_stdout(StringIO()), redirect_stderr(StringIO()):\n", + " _model = hssm.HSSM(\n", + " **_base_model_kwargs,\n", + " include=include,\n", + " **kwargs,\n", + " )\n", + " finally:\n", + " _logger.handlers = _old_handlers\n", + " _logger.setLevel(_old_level)\n", + " _logger.propagate = _old_propagate\n", + " _messages = tuple(\n", + " _line.strip() for _line in _log_stream.getvalue().splitlines() if _line\n", + " )\n", + " return _model, _messages\n", + "\n", + "\n", + "def expect_model_error(include, **kwargs):\n", + " \"\"\"Return the expected pre-build ValueError as stable tutorial evidence.\"\"\"\n", + " try:\n", + " build_model(include, **kwargs)\n", + " except ValueError as _error:\n", + " return str(_error)\n", + " raise AssertionError(\"HSSM unexpectedly built an incompatible model\")\n", + "\n", + "\n", + "def free_rv_names(model):\n", + " \"\"\"Return exact free-RV names; do not infer component names by splitting.\"\"\"\n", + " return {variable.name for variable in model.pymc_model.free_RVs}\n", + "\n", "\n", - "import hssm\n", - "from hssm import Prior\n", + "def group_term_structure(model, parameter, term):\n", + " \"\"\"Summarize one exact group key in the built PyMC graph.\"\"\"\n", + " _prefix = f\"{parameter}_{term}\"\n", + " _free = free_rv_names(model)\n", + " _prior = model.params[parameter].prior[term]\n", + " return {\n", + " \"parameter\": parameter,\n", + " \"group term\": term,\n", + " \"prior override\": getattr(_prior, \"noncentered\", None),\n", + " \"effective form\": (\n", + " \"non-centered\" if f\"{_prefix}_offset\" in _free else \"centered\"\n", + " ),\n", + " \"direct group RV\": _prefix in _free,\n", + " \"offset RV\": f\"{_prefix}_offset\" in _free,\n", + " \"free mu RV\": f\"{_prefix}_mu\" in _free,\n", + " \"free sigma RV\": f\"{_prefix}_sigma\" in _free,\n", + " \"disconnected RVs\": \", \".join(find_disconnected_free_rvs(model.pymc_model))\n", + " or \"none\",\n", + " }\n", "\n", - "hssm.set_floatX(\"float32\")\n", - "print(\"hssm\", hssm.__version__, \"| bambi\", bmb.__version__)" + "\n", + "def assert_connected(model):\n", + " \"\"\"Make every successful example double as a graph regression.\"\"\"\n", + " assert find_disconnected_free_rvs(model.pymc_model) == []\n", + "\n", + "\n", + "def model_graph(model):\n", + " \"\"\"Render a compact construction-only PyMC graph.\"\"\"\n", + " return pm.model_to_graphviz(\n", + " model.pymc_model,\n", + " graph_attr={\"bgcolor\": \"white\", \"rankdir\": \"LR\"},\n", + " )" ] }, { - "cell_type": "code", - "execution_count": 2, - "id": "ce06658c", + "cell_type": "markdown", + "id": "PKri", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:16.790769Z", - "iopub.status.busy": "2026-07-13T00:03:16.790697Z", - "iopub.status.idle": "2026-07-13T00:03:16.798778Z", - "shell.execute_reply": "2026-07-13T00:03:16.798415Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" } }, + "source": [ + "## Three levels of control\n", + "\n", + "HSSM passes a scalar `noncentered=True` or `False` to Bambi as the default\n", + "for every group term. A dictionary selects defaults by HSSM parameter name.\n", + "Missing dictionary keys fall back to `True`, not to the value of another\n", + "component. Finally, `noncentered` on an explicit prior wins for that one term.\n", + "\n", + "The formulas below contain matching common and group intercepts. The common\n", + "`Intercept` owns the population location, so each generated group intercept is\n", + "a mean-zero deviation and either parameterization is faithful." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "Xref", + "metadata": {}, + "outputs": [], + "source": [ + "hierarchical_specs = [\n", + " {\"name\": \"v\", \"formula\": \"v ~ 1 + (1 | participant_id)\"},\n", + " {\"name\": \"a\", \"formula\": \"a ~ 1 + (1 | participant_id)\"},\n", + "]\n", + "\n", + "_scalar_noncentered_model, _scalar_nc_messages = build_model(\n", + " hierarchical_specs,\n", + " noncentered=True,\n", + ")\n", + "_scalar_centered_model, _scalar_c_messages = build_model(\n", + " hierarchical_specs,\n", + " noncentered=False,\n", + ")\n", + "_component_dict_model, _component_messages = build_model(\n", + " hierarchical_specs,\n", + " noncentered={\"v\": False},\n", + ")\n", + "\n", + "parameterization_models = {\n", + " \"scalar True\": _scalar_noncentered_model,\n", + " \"scalar False\": _scalar_centered_model,\n", + " \"dict: v=False; a omitted\": _component_dict_model,\n", + "}\n", + "assert not (_scalar_nc_messages or _scalar_c_messages or _component_messages)\n", + "for _model in parameterization_models.values():\n", + " assert_connected(_model)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "SFPL", + "metadata": {}, "outputs": [ { "data": { "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
participant_idstimrtresponsethetadbsconf
00LL1.211.00.6562751HC
10WL1.631.0-0.3278891LC
20WW1.031.0-0.4802851HC
30WL2.771.01.9274271LC
40WW1.14-1.0-0.2132361HC
\n", - "
" - ], - "text/plain": [ - " participant_id stim rt response theta dbs conf\n", - "0 0 LL 1.21 1.0 0.656275 1 HC\n", - "1 0 WL 1.63 1.0 -0.327889 1 LC\n", - "2 0 WW 1.03 1.0 -0.480285 1 HC\n", - "3 0 WL 2.77 1.0 1.927427 1 LC\n", - "4 0 WW 1.14 -1.0 -0.213236 1 HC" + "
model settingparametereffective formdirect group RVoffset RVfree mu RVdisconnected RVs
0scalar Truevnon-centeredFalseTrueFalsenone
1scalar Trueanon-centeredFalseTrueFalsenone
2scalar FalsevcenteredTrueFalseFalsenone
3scalar FalseacenteredTrueFalseFalsenone
4dict: v=False; a omittedvcenteredTrueFalseFalsenone
5dict: v=False; a omittedanon-centeredFalseTrueFalsenone
" ] }, - "execution_count": 2, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "data = hssm.load_data(\"cavanagh_theta\")\n", - "data.head()" + "_rows = []\n", + "for _setting, _model in parameterization_models.items():\n", + " for _parameter in (\"v\", \"a\"):\n", + " _row = group_term_structure(_model, _parameter, \"1|participant_id\")\n", + " _row = {\"model setting\": _setting, **_row}\n", + " _rows.append(_row)\n", + "\n", + "parameterization_table = pd.DataFrame(_rows)[\n", + " [\n", + " \"model setting\",\n", + " \"parameter\",\n", + " \"effective form\",\n", + " \"direct group RV\",\n", + " \"offset RV\",\n", + " \"free mu RV\",\n", + " \"disconnected RVs\",\n", + " ]\n", + "]\n", + "assert parameterization_table.loc[\n", + " parameterization_table[\"model setting\"] == \"scalar True\", \"offset RV\"\n", + "].all()\n", + "assert parameterization_table.loc[\n", + " parameterization_table[\"model setting\"] == \"scalar False\", \"direct group RV\"\n", + "].all()\n", + "_dict_rows = parameterization_table[\n", + " parameterization_table[\"model setting\"] == \"dict: v=False; a omitted\"\n", + "].set_index(\"parameter\")\n", + "assert _dict_rows.loc[\"v\", \"effective form\"] == \"centered\"\n", + "assert _dict_rows.loc[\"a\", \"effective form\"] == \"non-centered\"\n", + "parameterization_table" ] }, { - "cell_type": "markdown", - "id": "fcdd93a8", + "cell_type": "code", + "execution_count": null, + "id": "BYtC", "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", + "\n", + "\n", + "clusterparticipant_id__factor_dim (4)\n", + "\n", + "participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (24)\n", + "\n", + "__obs__ (24)\n", + "\n", + "\n", + "cluster__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "a_1|participant_id_sigma\n", + "\n", + "a_1|participant_id_sigma\n", + "~\n", + "Weibull\n", + "\n", + "\n", + "\n", + "a_1|participant_id\n", + "\n", + "a_1|participant_id\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "a_1|participant_id_sigma->a_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "a_Intercept\n", + "\n", + "a_Intercept\n", + "~\n", + "Truncated\n", + "\n", + "\n", + "\n", + "a\n", + "\n", + "a\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "a_Intercept->a\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Weibull\n", + "\n", + "\n", + "\n", + "v_1|participant_id\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma->v_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_Intercept\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_Intercept->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "a_1|participant_id_offset\n", + "\n", + "a_1|participant_id_offset\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "a_1|participant_id_offset->a_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "a_1|participant_id->a\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "v->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "a->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "We will detect the parameterization **structurally**, without sampling: bambi\n", - "creates a standard-normal `_..._offset` random variable for each\n", - "*non-centered* group term (and a `Deterministic` equal to `offset * sigma`). A\n", - "*centered* term instead has a direct `_` random variable. So the\n", - "set of parameters that own an `*_offset` free RV is exactly the set of\n", - "non-centered parameters." + "mixed_parameterization_graph = model_graph(\n", + " parameterization_models[\"dict: v=False; a omitted\"]\n", + ")\n", + "mixed_parameterization_graph" ] }, { - "cell_type": "code", - "execution_count": 3, - "id": "b1060547", + "cell_type": "markdown", + "id": "RGSE", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:16.799736Z", - "iopub.status.busy": "2026-07-13T00:03:16.799683Z", - "iopub.status.idle": "2026-07-13T00:03:16.801245Z", - "shell.execute_reply": "2026-07-13T00:03:16.800917Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" } }, - "outputs": [], "source": [ - "def noncentered_params(model):\n", - " \"\"\"Return HSSM parameters whose group term uses the non-centered form.\"\"\"\n", - " names = [rv.name for rv in model.pymc_model.free_RVs]\n", - " return sorted({n.split(\"_\")[0] for n in names if \"_offset\" in n})" + "In the dictionary case, `v_1|participant_id` is sampled directly because\n", + "`v` is centered. The missing `a` key takes Bambi's default `True`, so\n", + "`a_1|participant_id_offset` is sampled and multiplied by its scale. Neither\n", + "graph contains a group `mu`: these are zero-mean deviations around their\n", + "matching common intercepts.\n", + "\n", + "The graph below fixes the mixed dictionary case as the same static view in\n", + "marimo, Molab, and the rendered documentation." ] }, { "cell_type": "markdown", - "id": "2ce0d69b", + "id": "Kclp", + "metadata": { + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" + } + }, + "source": [ + "## A per-prior override wins\n", + "\n", + "The finest control lives on the group prior itself. Here the component\n", + "dictionary requests centered `v`, while this one explicit prior requests\n", + "non-centering. Its `mu=0` and hierarchical `sigma` satisfy Bambi's current\n", + "non-centered contract." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "emfo", + "metadata": {}, + "outputs": [], + "source": [ + "zero_mean_noncentered_prior = hssm.Prior(\n", + " \"Normal\",\n", + " mu=0.0,\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + " noncentered=True,\n", + ")\n", + "_override_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + (1 | participant_id)\",\n", + " \"prior\": {\"1|participant_id\": zero_mean_noncentered_prior},\n", + " }\n", + "]\n", + "override_model, override_messages = build_model(\n", + " _override_include,\n", + " a=1.5,\n", + " noncentered={\"v\": False},\n", + ")\n", + "_override_free = free_rv_names(override_model)\n", + "assert not override_messages\n", + "assert \"v_1|participant_id_offset\" in _override_free\n", + "assert \"v_1|participant_id_mu\" not in _override_free\n", + "assert_connected(override_model)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "Hstk", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
parametergroup termprior overrideeffective formdirect group RVoffset RVfree mu RVfree sigma RVdisconnected RVs
0v1|participant_idTruenon-centeredFalseTrueFalseTruenone
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "override_structure_table = pd.DataFrame(\n", + " [group_term_structure(override_model, \"v\", \"1|participant_id\")]\n", + ")\n", + "override_structure_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "nWHF", "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterv_1|participant_id__factor_dim (4)\n", + "\n", + "v_1|participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (24)\n", + "\n", + "__obs__ (24)\n", + "\n", + "\n", + "cluster__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma\n", + "\n", + "v_1|participant_id_sigma\n", + "~\n", + "Halfnormal\n", + "\n", + "\n", + "\n", + "v_1|participant_id\n", + "\n", + "v_1|participant_id\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_1|participant_id_sigma->v_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_Intercept\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_Intercept->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_1|participant_id_offset\n", + "\n", + "v_1|participant_id_offset\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_1|participant_id_offset->v_1|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "v->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "## A hierarchical model — the default\n", + "faithful_noncentered_graph = model_graph(override_model)\n", + "faithful_noncentered_graph" + ] + }, + { + "cell_type": "markdown", + "id": "iLit", + "metadata": { + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" + } + }, + "source": [ + "The graph contains the group scale and standard-normal offset, but no\n", + "`v_1|participant_id_mu`. This is not an orphan-removal trick: zero is the\n", + "intended location because the common `v_Intercept` owns the population mean.\n", "\n", - "We give both `v` (drift rate) and `a` (boundary separation) a by-participant\n", - "random intercept. bambi's default is **non-centered everywhere**, so both\n", - "parameters get an `*_offset`." + "The complete resolution order is:" ] }, { "cell_type": "code", - "execution_count": 4, - "id": "e7f26c9f", + "execution_count": null, + "id": "ZHCJ", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
prioritycontrolscoperule
01per-prior noncenteredone explicit group termwins when True or False
12model-level component dictionaryall group terms for that HSSM parameternamed key wins; missing key defaults to True
23model-level scalarall group termsTrue is the default; False requests centering
3safe-policy overridegenerated unique group-only prioronly the generated location-owning termHSSM centers it to preserve its location
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "precedence_table = pd.DataFrame(\n", + " [\n", + " {\n", + " \"priority\": 1,\n", + " \"control\": \"per-prior noncentered\",\n", + " \"scope\": \"one explicit group term\",\n", + " \"rule\": \"wins when True or False\",\n", + " },\n", + " {\n", + " \"priority\": 2,\n", + " \"control\": \"model-level component dictionary\",\n", + " \"scope\": \"all group terms for that HSSM parameter\",\n", + " \"rule\": \"named key wins; missing key defaults to True\",\n", + " },\n", + " {\n", + " \"priority\": 3,\n", + " \"control\": \"model-level scalar\",\n", + " \"scope\": \"all group terms\",\n", + " \"rule\": \"True is the default; False requests centering\",\n", + " },\n", + " {\n", + " \"priority\": \"safe-policy override\",\n", + " \"control\": \"generated unique group-only prior\",\n", + " \"scope\": \"only the generated location-owning term\",\n", + " \"rule\": \"HSSM centers it to preserve its location\",\n", + " },\n", + " ]\n", + ")\n", + "precedence_table" + ] + }, + { + "cell_type": "markdown", + "id": "ROlb", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:16.802127Z", - "iopub.status.busy": "2026-07-13T00:03:16.802073Z", - "iopub.status.idle": "2026-07-13T00:03:16.936114Z", - "shell.execute_reply": "2026-07-13T00:03:16.935770Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" } }, + "source": [ + "The last row is deliberately separate from user precedence. HSSM may adapt\n", + "a prior that it generated itself, but an explicit user prior remains\n", + "authoritative and is never rewritten.\n", + "\n", + "## Incompatible explicit non-centered locations fail before build\n", + "\n", + "A free group mean and a fixed nonzero group mean would both be omitted from\n", + "Bambi's `offset * sigma` construction. HSSM therefore rejects both before a\n", + "PyMC model exists." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "qnkX", + "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model initialized successfully.\n" - ] + "data": { + "text/html": [ + "
explicit group mueffective settingHSSM resultreason
0Normal hyperprior (free)non-centeredValueError before Bambi/PyMC buildBambi would create mu but omit it from offset * sigma
11.0 (fixed nonzero)non-centeredValueError before Bambi/PyMC buildBambi would ignore the requested location
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "_free_mu_prior = hssm.Prior(\n", + " \"Normal\",\n", + " mu=hssm.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + ")\n", + "_nonzero_mu_prior = hssm.Prior(\n", + " \"Normal\",\n", + " mu=1.0,\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + ")\n", + "\n", + "_free_mu_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + (1 | participant_id)\",\n", + " \"prior\": {\"1|participant_id\": _free_mu_prior},\n", + " }\n", + "]\n", + "_nonzero_mu_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + (1 | participant_id)\",\n", + " \"prior\": {\"1|participant_id\": _nonzero_mu_prior},\n", + " }\n", + "]\n", + "free_mu_error = expect_model_error(\n", + " _free_mu_include,\n", + " a=1.5,\n", + " noncentered=True,\n", + ")\n", + "nonzero_mu_error = expect_model_error(\n", + " _nonzero_mu_include,\n", + " a=1.5,\n", + " noncentered=True,\n", + ")\n", + "\n", + "assert free_mu_error.startswith(\n", + " \"Explicit group-specific prior specification(s) cannot be represented\"\n", + ")\n", + "assert \"mu` hyperprior\" in free_mu_error\n", + "assert \"disconnected node\" in free_mu_error\n", + "assert \"not fixed entirely to zero\" in nonzero_mu_error\n", + "assert \"silently ignored\" in nonzero_mu_error\n", + "\n", + "incompatible_location_table = pd.DataFrame(\n", + " [\n", + " {\n", + " \"explicit group mu\": \"Normal hyperprior (free)\",\n", + " \"effective setting\": \"non-centered\",\n", + " \"HSSM result\": \"ValueError before Bambi/PyMC build\",\n", + " \"reason\": \"Bambi would create mu but omit it from offset * sigma\",\n", + " },\n", + " {\n", + " \"explicit group mu\": \"1.0 (fixed nonzero)\",\n", + " \"effective setting\": \"non-centered\",\n", + " \"HSSM result\": \"ValueError before Bambi/PyMC build\",\n", + " \"reason\": \"Bambi would ignore the requested location\",\n", + " },\n", + " ]\n", + ")\n", + "incompatible_location_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "TqIu", + "metadata": { + "jupyter": { + "source_hidden": true }, + "marimo": { + "config": { + "hide_code": true + } + }, + "tags": [ + "remove-input" + ] + }, + "outputs": [ { "data": { - "text/plain": [ - "['a', 'v']" + "text/markdown": [ + "HSSM reports both the underlying limitation and term-specific repairs. The\n", + "free-mean case says:\n", + "\n", + "```text\n", + "Explicit group-specific prior specification(s) cannot be represented faithfully by bambi:\n", + "- User prior for group term '1|participant_id' on parameter 'v' is incompatible with the effective parameterization: outer prior supplies a `mu` hyperprior that bambi creates and then omits from `offset * sigma`, leaving a disconnected node. Continuing would either fail in bambi or change the requested prior. Keep the common formula term '1' and use a plain built-in Normal with hierarchical `sigma`, absent or fixed-all-zero `mu`, and no additional arguments for its zero-mean group deviation. To retain the explicit prior instead, set `noncentered=False` on this prior and on any nested hierarchical hyperpriors (or remove their overrides and make the effective component setting for 'v' centered); remove the common effect as well if a free group mean should own the population location.\n", + "```\n", + "\n", + "For a matched common/group expression, keep the common effect and use a\n", + "zero-mean group deviation. If the group distribution should own a free\n", + "population location, remove the matching common term and center that group\n", + "prior intentionally." ] }, - "execution_count": 4, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "hierarchical = [\n", - " {\"name\": \"v\", \"formula\": \"v ~ 1 + (1|participant_id)\"},\n", - " {\"name\": \"a\", \"formula\": \"a ~ 1 + (1|participant_id)\"},\n", - "]\n", + "mo.md(f\"\"\"\n", + "HSSM reports both the underlying limitation and term-specific repairs. The\n", + "free-mean case says:\n", "\n", - "model_default = hssm.HSSM(data=data, model=\"ddm\", include=hierarchical, p_outlier=0.0)\n", - "noncentered_params(model_default) # -> ['a', 'v']" + "```text\n", + "{free_mu_error}\n", + "```\n", + "\n", + "For a matched common/group expression, keep the common effect and use a\n", + "zero-mean group deviation. If the group distribution should own a free\n", + "population location, remove the matching common term and center that group\n", + "prior intentionally.\n", + "\"\"\")" ] }, { "cell_type": "markdown", - "id": "46a0b5cf", - "metadata": {}, + "id": "Vxnm", + "metadata": { + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" + } + }, "source": [ - "## Choose per parameter with a dict\n", + "## A unique group-only term owns its location\n", "\n", - "Pass `noncentered` as a `dict` keyed by **HSSM parameter name**. Here we make\n", - "`v` centered while keeping `a` non-centered. Only `a` keeps an `*_offset`." + "Now `theta` appears only inside `(0 + theta | participant_id)`. The group\n", + "distribution must estimate the population slope: fixing its mean to zero\n", + "would change the scientific model. With `prior_settings=\"safe\"`, HSSM keeps\n", + "that generated location and centers just this term—even when the model-level\n", + "request is non-centered." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "DnEU", + "metadata": {}, + "outputs": [], + "source": [ + "_group_only_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + (0 + theta | participant_id)\",\n", + " }\n", + "]\n", + "generated_owner_model, generated_owner_messages = build_model(\n", + " _group_only_include,\n", + " a=1.5,\n", + " noncentered=True,\n", + ")\n", + "_owner_prior = generated_owner_model.params[\"v\"].prior[\"theta|participant_id\"]\n", + "_owner_free = free_rv_names(generated_owner_model)\n", + "assert _owner_prior.noncentered is False\n", + "assert \"v_theta|participant_id\" in _owner_free\n", + "assert \"v_theta|participant_id_mu\" in _owner_free\n", + "assert \"v_theta|participant_id_offset\" not in _owner_free\n", + "assert len(generated_owner_messages) == 1\n", + "assert \"generated location-bearing group-only term\" in generated_owner_messages[0]\n", + "assert \"Explicit priors were not changed\" in generated_owner_messages[0]\n", + "assert_connected(generated_owner_model)" ] }, { "cell_type": "code", - "execution_count": 5, - "id": "96b882f9", + "execution_count": null, + "id": "ulZA", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:16.937624Z", - "iopub.status.busy": "2026-07-13T00:03:16.937300Z", - "iopub.status.idle": "2026-07-13T00:03:17.003536Z", - "shell.execute_reply": "2026-07-13T00:03:17.003088Z" - } + "jupyter": { + "source_hidden": true + }, + "marimo": { + "config": { + "hide_code": true + } + }, + "tags": [ + "remove-input" + ] }, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model initialized successfully.\n" - ] - }, + "data": { + "text/markdown": [ + "HSSM makes the generated fallback visible:\n", + "\n", + "```text\n", + "Safe priors for parameter 'v' generated location-bearing group-only term(s) [\"'theta|participant_id' (expression 'theta')\"]. The effective model/component setting requested noncentered=True, but Bambi's current non-centered construction cannot retain these group locations. HSSM set noncentered=False on these generated priors, preserving their locations on the response/parameter scale. Explicit priors were not changed. To use non-centering, add the exact common formula term(s) ['theta'] so the group effects become zero-mean deviations.\n", + "```\n", + "\n", + "In the graph, `v_theta|participant_id_mu` and the scale both feed the direct\n", + "group coefficient. There is no offset and no disconnected node." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mo.md(f\"\"\"\n", + "HSSM makes the generated fallback visible:\n", + "\n", + "```text\n", + "{generated_owner_messages[0]}\n", + "```\n", + "\n", + "In the graph, `v_theta|participant_id_mu` and the scale both feed the direct\n", + "group coefficient. There is no offset and no disconnected node.\n", + "\"\"\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecfG", + "metadata": {}, + "outputs": [ { "data": { - "text/plain": [ - "['a']" + "text/html": [ + "
parametergroup termprior overrideeffective formdirect group RVoffset RVfree mu RVfree sigma RVdisconnected RVs
0vtheta|participant_idFalsecenteredTrueFalseTrueTruenone
" ] }, - "execution_count": 5, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "model_mixed = hssm.HSSM(\n", - " data=data,\n", - " model=\"ddm\",\n", - " include=hierarchical,\n", - " p_outlier=0.0,\n", - " noncentered={\"v\": False, \"a\": True},\n", + "generated_owner_table = pd.DataFrame(\n", + " [group_term_structure(generated_owner_model, \"v\", \"theta|participant_id\")]\n", ")\n", - "noncentered_params(model_mixed) # -> ['a'] (v is now centered)" + "generated_owner_table" ] }, { - "cell_type": "markdown", - "id": "b407990f", + "cell_type": "code", + "execution_count": null, + "id": "Pvdt", "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterv_theta|participant_id__factor_dim (4)\n", + "\n", + "v_theta|participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (24)\n", + "\n", + "__obs__ (24)\n", + "\n", + "\n", + "cluster__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_mu\n", + "\n", + "v_theta|participant_id_mu\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_theta|participant_id\n", + "\n", + "v_theta|participant_id\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_mu->v_theta|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_sigma\n", + "\n", + "v_theta|participant_id_sigma\n", + "~\n", + "Weibull\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_sigma->v_theta|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_Intercept\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_Intercept->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_theta|participant_id->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "v->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "## Override a single term with a per-prior `noncentered`\n", - "\n", - "For the finest control, attach `noncentered` directly to a term's prior. It\n", - "**overrides** the model-level setting for just that term. The precedence is:\n", + "generated_owner_graph = model_graph(generated_owner_model)\n", + "generated_owner_graph" + ] + }, + { + "cell_type": "markdown", + "id": "ZBYS", + "metadata": { + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" + } + }, + "source": [ + "## Centered does not automatically mean identifiable\n", "\n", - "```\n", - "per-prior noncentered > model-level noncentered dict > default (True)\n", + "Centering makes a free group `mu` part of the model, but the formula still\n", + "needs exactly one population-location owner. If common `theta` and the mean\n", + "of `theta|participant_id` are both free, the likelihood sees only their sum.\n", + "HSSM can build this model faithfully, so it warns about the location ridge\n", + "instead of rejecting the prior." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aLJB", + "metadata": {}, + "outputs": [], + "source": [ + "_matched_free_location = hssm.Prior(\n", + " \"Normal\",\n", + " mu=hssm.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", + ")\n", + "_matched_ridge_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": \"v ~ 1 + theta + (0 + theta | participant_id)\",\n", + " \"prior\": {\"theta|participant_id\": _matched_free_location},\n", + " }\n", + "]\n", + "matched_ridge_model, matched_ridge_messages = build_model(\n", + " _matched_ridge_include,\n", + " a=1.5,\n", + " noncentered=False,\n", + ")\n", + "assert len(matched_ridge_messages) == 1\n", + "assert \"non-identifiable\" in matched_ridge_messages[0]\n", + "assert \"common 'theta'\" in matched_ridge_messages[0]\n", + "assert \"disconnected\" not in matched_ridge_messages[0].lower()\n", + "assert_connected(matched_ridge_model)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "nHfw", + "metadata": { + "jupyter": { + "source_hidden": true + }, + "marimo": { + "config": { + "hide_code": true + } + }, + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "text/markdown": [ + "```text\n", + "User prior for 'theta|participant_id' on parameter 'v' has a free `mu`, and its Formulae expression 'theta' also occurs as a common effect under the effective centered parameterization. The data only constrains their sum; the common and group locations are non-identifiable individually and the posterior will have a ridge along the anti-diagonal. Keep the common 'theta' effect and set `mu=0` on the matching group term, or remove that common effect if the group-level mean should own the location.\n", + "```\n", + "\n", + "The graph is connected, but connectivity is not identifiability. Both\n", + "`v_theta` and `v_theta|participant_id_mu` shift the same predictor. Keep\n", + "common `theta` and set the group `mu=0`, or remove common `theta` and let the\n", + "centered group distribution own the location." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mo.md(f\"\"\"\n", + "```text\n", + "{matched_ridge_messages[0]}\n", "```\n", "\n", - "This works whether you specify the prior as a plain `dict` or as an\n", - "`hssm.Prior` object. Below, the model-level dict asks for centered `v`, but the\n", - "per-prior field wins and makes it non-centered." + "The graph is connected, but connectivity is not identifiability. Both\n", + "`v_theta` and `v_theta|participant_id_mu` shift the same predictor. Keep\n", + "common `theta` and set the group `mu=0`, or remove common `theta` and let the\n", + "centered group distribution own the location.\n", + "\"\"\")" ] }, { "cell_type": "code", - "execution_count": 6, - "id": "559925f4", + "execution_count": null, + "id": "xXTn", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
parametergroup termprior overrideeffective formdirect group RVoffset RVfree mu RVfree sigma RVdisconnected RVs
0vtheta|participant_idNonecenteredTrueFalseTrueTruenone
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "matched_ridge_table = pd.DataFrame(\n", + " [group_term_structure(matched_ridge_model, \"v\", \"theta|participant_id\")]\n", + ")\n", + "matched_ridge_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "AjVT", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterv_theta|participant_id__factor_dim (4)\n", + "\n", + "v_theta|participant_id__factor_dim (4)\n", + "\n", + "\n", + "cluster__obs__ (24)\n", + "\n", + "__obs__ (24)\n", + "\n", + "\n", + "cluster__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "__obs__ (24) x rt,response_extra_dim_0 (2)\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_mu\n", + "\n", + "v_theta|participant_id_mu\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_theta|participant_id\n", + "\n", + "v_theta|participant_id\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_mu->v_theta|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_theta\n", + "\n", + "v_theta\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v\n", + "\n", + "v\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "v_theta->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_sigma\n", + "\n", + "v_theta|participant_id_sigma\n", + "~\n", + "Halfnormal\n", + "\n", + "\n", + "\n", + "v_theta|participant_id_sigma->v_theta|participant_id\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_Intercept\n", + "\n", + "v_Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "v_Intercept->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "v_theta|participant_id->v\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "rt,response\n", + "\n", + "rt,response\n", + "~\n", + "Ddm_RV\n", + "\n", + "\n", + "\n", + "v->rt,response\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "matched_ridge_graph = model_graph(matched_ridge_model)\n", + "matched_ridge_graph" + ] + }, + { + "cell_type": "markdown", + "id": "pHFh", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:17.004607Z", - "iopub.status.busy": "2026-07-13T00:03:17.004469Z", - "iopub.status.idle": "2026-07-13T00:03:17.038070Z", - "shell.execute_reply": "2026-07-13T00:03:17.037782Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" } }, + "source": [ + "### Repeated group-only owners have the same ridge\n", + "\n", + "A different ambiguity appears when the same unmatched expression has free\n", + "means under multiple grouping factors. Shifting every participant effect up\n", + "and every `conf` effect down leaves the predictor unchanged. Explicit priors\n", + "remain authoritative, so HSSM builds the centered model and emits one\n", + "aggregated warning." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "NCOB", + "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model initialized successfully.\n" - ] - }, { "data": { - "text/plain": [ - "['v']" + "text/html": [ + "
parametergroup termprior overrideeffective formdirect group RVoffset RVfree mu RVfree sigma RVdisconnected RVs
0vtheta|participant_idNonecenteredTrueFalseTrueTruenone
1vtheta|confNonecenteredTrueFalseTrueTruenone
" ] }, - "execution_count": 6, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "prior = {\n", - " \"1|participant_id\": Prior(\n", + "def _fresh_free_location():\n", + " return hssm.Prior(\n", " \"Normal\",\n", - " mu=0.0,\n", - " sigma=Prior(\"HalfNormal\", sigma=1.0),\n", - " noncentered=True, # overrides the model-level setting for this term\n", + " mu=hssm.Prior(\"Normal\", mu=0.0, sigma=0.5),\n", + " sigma=hssm.Prior(\"HalfNormal\", sigma=0.5),\n", " )\n", - "}\n", "\n", - "model_override = hssm.HSSM(\n", - " data=data,\n", - " model=\"ddm\",\n", - " include=[{\"name\": \"v\", \"formula\": \"v ~ 1 + (1|participant_id)\", \"prior\": prior}],\n", - " p_outlier=0.0,\n", - " noncentered={\"v\": False}, # model says \"centered\"; the per-prior field wins\n", + "\n", + "_repeated_owner_include = [\n", + " {\n", + " \"name\": \"v\",\n", + " \"formula\": (\"v ~ 1 + (0 + theta | participant_id) + (0 + theta | conf)\"),\n", + " \"prior\": {\n", + " \"theta|participant_id\": _fresh_free_location(),\n", + " \"theta|conf\": _fresh_free_location(),\n", + " },\n", + " }\n", + "]\n", + "repeated_owner_model, repeated_owner_messages = build_model(\n", + " _repeated_owner_include,\n", + " a=1.5,\n", + " noncentered=False,\n", + ")\n", + "assert len(repeated_owner_messages) == 1\n", + "assert \"identified only by the priors\" in repeated_owner_messages[0]\n", + "assert \"theta|participant_id\" in repeated_owner_messages[0]\n", + "assert \"theta|conf\" in repeated_owner_messages[0]\n", + "assert_connected(repeated_owner_model)\n", + "\n", + "repeated_owner_table = pd.DataFrame(\n", + " [\n", + " group_term_structure(\n", + " repeated_owner_model,\n", + " \"v\",\n", + " _term,\n", + " )\n", + " for _term in (\"theta|participant_id\", \"theta|conf\")\n", + " ]\n", ")\n", - "noncentered_params(model_override) # -> ['v']" + "repeated_owner_table" ] }, { - "cell_type": "markdown", - "id": "62316c4c", - "metadata": {}, + "cell_type": "code", + "execution_count": null, + "id": "aqbW", + "metadata": { + "jupyter": { + "source_hidden": true + }, + "marimo": { + "config": { + "hide_code": true + } + }, + "tags": [ + "remove-input" + ] + }, + "outputs": [ + { + "data": { + "text/markdown": [ + "```text\n", + "User priors for group terms ['theta|conf', 'theta|participant_id'] on parameter 'v' each have a free `mu` under the effective centered parameterization, and their exact Formulae expression 'theta' has no common effect. The likelihood is invariant when one group location is shifted up and another is shifted down; their decomposition is identified only by the priors, and the likelihood has a location ridge. Proper priors may still yield a proper posterior. Add the exact common formula term 'theta' and set `mu=0` on every matching group deviation, or choose exactly one group term to own the free population location and fix the other group locations intentionally.\n", + "```\n", + "\n", + "A cleaner formula adds common `theta` and uses mean-zero deviations for both\n", + "grouping factors. Alternatively, choose exactly one centered group term as\n", + "the location owner and fix the other group location intentionally." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "## Caveats\n", - "\n", - "- **Hierarchical terms only.** `noncentered` affects group-specific terms\n", - " (e.g. `... + (1|participant_id)`). Setting it for a parameter with no group\n", - " term is a silent no-op.\n", - "- **Unknown keys fail fast.** A typo'd parameter name in the dict raises at\n", - " construction, listing the valid names.\n", - "- **`mu` is dropped under non-centered.** The non-centered reparameterization is\n", - " `offset * sigma`; a non-zero `mu` on a group `Normal` prior is ignored when\n", - " that term is non-centered. Put the location on the common `Intercept` instead — see [the explanation page](https://lnccbrown.github.io/HSSM/tutorials/centered_vs_noncentered_basic_logic/) for the full footgun anatomy.\n", - "- **Bounded group priors.** A bounded/truncated prior on a group term is\n", - " incompatible with bambi's group-hyperprior requirement regardless of\n", - " `noncentered` — unrelated to this feature, but worth knowing.\n", + "mo.md(f\"\"\"\n", + "```text\n", + "{repeated_owner_messages[0]}\n", + "```\n", "\n", - "The unknown-key guard in action:" + "A cleaner formula adds common `theta` and uses mean-zero deviations for both\n", + "grouping factors. Alternatively, choose exactly one centered group term as\n", + "the location owner and fix the other group location intentionally.\n", + "\"\"\")" ] }, { - "cell_type": "code", - "execution_count": 7, - "id": "66dbce4a", + "cell_type": "markdown", + "id": "TRpd", "metadata": { - "execution": { - "iopub.execute_input": "2026-07-13T00:03:17.039344Z", - "iopub.status.busy": "2026-07-13T00:03:17.039220Z", - "iopub.status.idle": "2026-07-13T00:03:17.044738Z", - "shell.execute_reply": "2026-07-13T00:03:17.044333Z" + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" } }, + "source": [ + "## Configuration errors fail fast too\n", + "\n", + "Component dictionaries are keyed by HSSM parameter names. A typo does not\n", + "silently fall back to a default: HSSM/Bambi list the valid names during\n", + "construction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "TXez", + "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Unknown component name(s) in `noncentered`: ['vv']. Valid component names for this model: ['a', 't', 'v', 'z'].\n" - ] + "data": { + "text/html": [ + "
inputresultmessage
0noncentered={'vv': True}ValueError at constructionUnknown component name(s) in `noncentered`: ['vv']. Valid component names for this model: ['a', 't', 'v', 'z'].
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "try:\n", - " hssm.HSSM(\n", - " data=data,\n", - " model=\"ddm\",\n", - " include=[{\"name\": \"v\", \"formula\": \"v ~ 1 + (1|participant_id)\"}],\n", - " p_outlier=0.0,\n", - " noncentered={\"vv\": True}, # typo: there is no parameter \"vv\"\n", - " )\n", - "except ValueError as err:\n", - " print(err)" + "_unknown_key_include = [{\"name\": \"v\", \"formula\": \"v ~ 1 + (1 | participant_id)\"}]\n", + "unknown_key_error = expect_model_error(\n", + " _unknown_key_include,\n", + " a=1.5,\n", + " noncentered={\"vv\": True},\n", + ")\n", + "assert \"Unknown component name(s)\" in unknown_key_error\n", + "assert \"['a', 't', 'v', 'z']\" in unknown_key_error\n", + "unknown_key_table = pd.DataFrame(\n", + " [\n", + " {\n", + " \"input\": \"noncentered={'vv': True}\",\n", + " \"result\": \"ValueError at construction\",\n", + " \"message\": unknown_key_error,\n", + " }\n", + " ]\n", + ")\n", + "unknown_key_table" + ] + }, + { + "cell_type": "markdown", + "id": "dNNg", + "metadata": { + "marimo": { + "config": { + "hide_code": true + }, + "md_prefix": "" + } + }, + "source": [ + "## Practical rules\n", + "\n", + "| Formula/prior situation | Recommended action |\n", + "| --- | --- |\n", + "| Common and group expressions match | Let the common term own the location; use group `mu=0`; choose centered or non-centered for sampling geometry. |\n", + "| One generated group-only expression | Let that group distribution own the location; HSSM's safe prior centers it automatically and explains the fallback. |\n", + "| Explicit group-only prior with a free location | Set that prior effectively centered so its `mu` is retained. |\n", + "| Explicit free or nonzero `mu` under non-centering | Change the formula/prior ownership or center the term; HSSM rejects the incompatible specification. |\n", + "| Two free centered owners for one expression | Add the exact common term and use zero-mean deviations, or choose exactly one owner. |\n", + "\n", + "Two final boundaries are worth remembering:\n", + "\n", + "- `noncentered` only affects group-specific terms. A component with no group\n", + " term is unchanged.\n", + "- HSSM's truncated group-prior wrapper is incompatible with Bambi's\n", + " hierarchical group-prior contract under either parameterization. Use an\n", + " untruncated hierarchical coefficient prior and a support-respecting link.\n", + "\n", + "Continue with:\n", + "\n", + "- [Specify hierarchical group priors](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) for the exact compatibility and location-ownership rules;\n", + "- [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/) for predictor versus response scale; and\n", + "- [Centered vs. non-centered parameterizations](https://lnccbrown.github.io/HSSM/tutorials/centered_vs_noncentered_basic_logic/) for the underlying graph transformation.\n", + "\n", + "For the sampling geometry behind this choice, Michael Betancourt's\n", + "[Hierarchical Modeling](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html)\n", + "case study provides a detailed treatment." ] } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -494,8 +1763,14 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.14.6" + "pygments_lexer": "ipython3" + }, + "marimo": { + "app_config": { + "width": "medium" + }, + "header": "# /// script\n# requires-python = \">=3.12,<3.15\"\n# dependencies = [\n# \"bambi==0.20.0\",\n# \"graphviz==0.21\",\n# \"hssm @ git+https://github.com/lnccbrown/HSSM.git@b6a6bdcf68ecd7cf71ffdef5cdda4fb05e8bfaad\",\n# \"marimo==0.24.0\",\n# \"numpy==2.4.6\",\n# \"pandas==3.0.5\",\n# \"pymc==6.3.1\",\n# ]\n# ///\n\n\"\"\"Choose centered or non-centered group effects at each HSSM control level.\n\nThis construction-only marimo tutorial is the source of truth for the rendered\nJupyter artifact published by MkDocs. It uses tiny synthetic data, inspects the\nPyMC graphs, and exercises HSSM's current group-prior preflight without sampling.\n\nRun the pinned standalone environment locally or in Molab::\n\n uvx marimo edit --sandbox docs/tutorials/parameterization_per_parameter.py\n\nTo exercise an active HSSM checkout instead, ignore the inline environment::\n\n uv run --group notebook --group docs marimo edit --no-sandbox \\\n docs/tutorials/parameterization_per_parameter.py\n uv run --group notebook --group docs marimo check --strict \\\n docs/tutorials/parameterization_per_parameter.py\n uv run --group notebook --group docs marimo export html --no-sandbox \\\n docs/tutorials/parameterization_per_parameter.py \\\n --output /tmp/parameterization-per-parameter.html --force\n uv run --group notebook --group docs marimo export ipynb --no-sandbox \\\n docs/tutorials/parameterization_per_parameter.py \\\n --output docs/tutorials/parameterization_per_parameter.ipynb \\\n --include-outputs --force\n uv run ruff format docs/tutorials/parameterization_per_parameter.ipynb\n\nGraph cells require the Graphviz ``dot`` executable (for example,\n``brew install graphviz`` on macOS). No sampling is performed.\n\"\"\"\n\n# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns)\n", + "marimo_version": "0.24.0" } }, "nbformat": 4, diff --git a/docs/tutorials/parameterization_per_parameter.py b/docs/tutorials/parameterization_per_parameter.py new file mode 100644 index 000000000..d61350540 --- /dev/null +++ b/docs/tutorials/parameterization_per_parameter.py @@ -0,0 +1,883 @@ +# /// script +# requires-python = ">=3.12,<3.15" +# dependencies = [ +# "bambi==0.20.0", +# "graphviz==0.21", +# "hssm @ git+https://github.com/lnccbrown/HSSM.git@b6a6bdcf68ecd7cf71ffdef5cdda4fb05e8bfaad", +# "marimo==0.24.0", +# "numpy==2.4.6", +# "pandas==3.0.5", +# "pymc==6.3.1", +# ] +# /// + +"""Choose centered or non-centered group effects at each HSSM control level. + +This construction-only marimo tutorial is the source of truth for the rendered +Jupyter artifact published by MkDocs. It uses tiny synthetic data, inspects the +PyMC graphs, and exercises HSSM's current group-prior preflight without sampling. + +Run the pinned standalone environment locally or in Molab:: + + uvx marimo edit --sandbox docs/tutorials/parameterization_per_parameter.py + +To exercise an active HSSM checkout instead, ignore the inline environment:: + + uv run --group notebook --group docs marimo edit --no-sandbox \ + docs/tutorials/parameterization_per_parameter.py + uv run --group notebook --group docs marimo check --strict \ + docs/tutorials/parameterization_per_parameter.py + uv run --group notebook --group docs marimo export html --no-sandbox \ + docs/tutorials/parameterization_per_parameter.py \ + --output /tmp/parameterization-per-parameter.html --force + uv run --group notebook --group docs marimo export ipynb --no-sandbox \ + docs/tutorials/parameterization_per_parameter.py \ + --output docs/tutorials/parameterization_per_parameter.ipynb \ + --include-outputs --force + uv run ruff format docs/tutorials/parameterization_per_parameter.ipynb + +Graph cells require the Graphviz ``dot`` executable (for example, +``brew install graphviz`` on macOS). No sampling is performed. +""" + +# ruff: noqa: B018, D401, E501, PLR1711 (generated marimo notebook: prose, cell display expressions, and bare returns) +import marimo + +__generated_with = "0.24.0" +app = marimo.App(width="medium") + + +@app.cell +def _(): + import logging + import os + import warnings + from contextlib import redirect_stderr, redirect_stdout + from io import StringIO + from tempfile import gettempdir + + # This notebook only constructs graphs. Molab and some development hosts can + # expose an unusable CUDA plugin, so prevent JAX from probing it. + os.environ["JAX_PLATFORMS"] = "cpu" + os.environ["JAX_SKIP_CUDA_CONSTRAINTS_CHECK"] = "1" + os.environ.setdefault( + "MPLCONFIGDIR", f"{gettempdir()}/hssm-per-parameter-matplotlib" + ) + warnings.filterwarnings("ignore") + logging.getLogger("jax._src.xla_bridge").setLevel(logging.CRITICAL) + + import bambi as bmb + import jax + import marimo as mo + import numpy as np + import pandas as pd + import pymc as pm + + # HSSM reports backend configuration and registry details during import and + # setup. Suppress that incidental output; later helpers capture HSSM's + # parameterization warnings explicitly and display them in the relevant cell. + with redirect_stdout(StringIO()), redirect_stderr(StringIO()): + import hssm + from hssm.param.parameterization_check import find_disconnected_free_rvs + + hssm.set_floatX("float64") + + assert jax.default_backend() == "cpu" + pd.set_option("display.max_colwidth", 100) + return ( + StringIO, + bmb, + find_disconnected_free_rvs, + hssm, + jax, + logging, + mo, + np, + pd, + pm, + redirect_stderr, + redirect_stdout, + ) + + +@app.cell(hide_code=True) +def _(bmb, hssm, jax, mo): + mo.md(f""" + # Choosing a parameterization per parameter + + Hierarchical models often sample better when different parameters—or even + different group terms—use different parameterizations. This tutorial shows + how HSSM resolves those choices and how to verify the resulting PyMC graph. + + By the end, you will be able to: + + 1. choose a model-wide or component-specific default; + 2. override one explicit group prior safely; + 3. recognize HSSM's pre-build errors and location-ridge warnings; and + 4. decide which formula term owns each population location. + + Everything below is structural: the models are built but never sampled. + + **Environment:** HSSM `{hssm.__version__}`, Bambi `{bmb.__version__}`, + JAX `{jax.__version__}` on `{jax.default_backend()}`. + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Centering is both geometry and model structure + + A textbook hierarchical Normal can be written in centered form, + + \[ + b_g \sim \mathcal N(\mu, \sigma), + \] + + or in the mathematically equivalent non-centered form, + + \[ + z_g \sim \mathcal N(0,1), \qquad b_g = \mu + \sigma z_g. + \] + + The likelihood can be identical while the posterior geometry—and therefore + sampling efficiency—changes. Non-centering is often helpful for weakly + informed groups; centering can be better for strongly informed groups. + + There is one important implementation boundary. Current Bambi constructs a + non-centered group term as + + \[ + b_g = \sigma z_g, + \] + + so this route faithfully represents a built-in Normal group prior only when + its location is absent or fixed entirely to zero and `sigma` is hierarchical. + HSSM checks explicit priors before asking Bambi to build the model. A free or + nonzero `mu` under effective non-centering now raises an actionable error; it + is not silently accepted and then discarded. + + The [hierarchical group-prior guide](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) + gives the complete compatibility table. For the scale on which these effects + combine, see [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/). + """) + return + + +@app.cell +def _(np, pd): + _trial = np.arange(24) + tutorial_data = pd.DataFrame( + { + "rt": 0.42 + 0.01 * _trial, + "response": np.where(_trial % 2, 1, -1), + "theta": np.tile([-1.0, 0.0, 1.0], 8), + "participant_id": np.repeat(np.arange(4), 6), + "conf": np.tile(["low", "high"], 12), + } + ) + tutorial_data.head(8) + return (tutorial_data,) + + +@app.cell +def _( + StringIO, + find_disconnected_free_rvs, + hssm, + logging, + pm, + redirect_stderr, + redirect_stdout, + tutorial_data, +): + _base_model_kwargs = { + "data": tutorial_data, + "model": "ddm", + "loglik_kind": "analytical", + "p_outlier": 0.0, + "prior_settings": "safe", + "process_initvals": False, + "initval_jitter": 0.0, + "z": 0.5, + "t": 0.2, + } + + def build_model(include, **kwargs): + """Build quietly while returning only HSSM warning messages.""" + _log_stream = StringIO() + _handler = logging.StreamHandler(_log_stream) + _handler.setFormatter(logging.Formatter("%(message)s")) + _logger = logging.getLogger("hssm") + _old_handlers = list(_logger.handlers) + _old_level = _logger.level + _old_propagate = _logger.propagate + _logger.handlers = [_handler] + _logger.setLevel(logging.WARNING) + _logger.propagate = False + try: + with redirect_stdout(StringIO()), redirect_stderr(StringIO()): + _model = hssm.HSSM( + **_base_model_kwargs, + include=include, + **kwargs, + ) + finally: + _logger.handlers = _old_handlers + _logger.setLevel(_old_level) + _logger.propagate = _old_propagate + _messages = tuple( + _line.strip() for _line in _log_stream.getvalue().splitlines() if _line + ) + return _model, _messages + + def expect_model_error(include, **kwargs): + """Return the expected pre-build ValueError as stable tutorial evidence.""" + try: + build_model(include, **kwargs) + except ValueError as _error: + return str(_error) + raise AssertionError("HSSM unexpectedly built an incompatible model") + + def free_rv_names(model): + """Return exact free-RV names; do not infer component names by splitting.""" + return {variable.name for variable in model.pymc_model.free_RVs} + + def group_term_structure(model, parameter, term): + """Summarize one exact group key in the built PyMC graph.""" + _prefix = f"{parameter}_{term}" + _free = free_rv_names(model) + _prior = model.params[parameter].prior[term] + return { + "parameter": parameter, + "group term": term, + "prior override": getattr(_prior, "noncentered", None), + "effective form": ( + "non-centered" if f"{_prefix}_offset" in _free else "centered" + ), + "direct group RV": _prefix in _free, + "offset RV": f"{_prefix}_offset" in _free, + "free mu RV": f"{_prefix}_mu" in _free, + "free sigma RV": f"{_prefix}_sigma" in _free, + "disconnected RVs": ", ".join(find_disconnected_free_rvs(model.pymc_model)) + or "none", + } + + def assert_connected(model): + """Make every successful example double as a graph regression.""" + assert find_disconnected_free_rvs(model.pymc_model) == [] + + def model_graph(model): + """Render a compact construction-only PyMC graph.""" + return pm.model_to_graphviz( + model.pymc_model, + graph_attr={"bgcolor": "white", "rankdir": "LR"}, + ) + + return ( + assert_connected, + build_model, + expect_model_error, + free_rv_names, + group_term_structure, + model_graph, + ) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## Three levels of control + + HSSM passes a scalar `noncentered=True` or `False` to Bambi as the default + for every group term. A dictionary selects defaults by HSSM parameter name. + Missing dictionary keys fall back to `True`, not to the value of another + component. Finally, `noncentered` on an explicit prior wins for that one term. + + The formulas below contain matching common and group intercepts. The common + `Intercept` owns the population location, so each generated group intercept is + a mean-zero deviation and either parameterization is faithful. + """) + return + + +@app.cell +def _(assert_connected, build_model): + hierarchical_specs = [ + {"name": "v", "formula": "v ~ 1 + (1 | participant_id)"}, + {"name": "a", "formula": "a ~ 1 + (1 | participant_id)"}, + ] + + _scalar_noncentered_model, _scalar_nc_messages = build_model( + hierarchical_specs, + noncentered=True, + ) + _scalar_centered_model, _scalar_c_messages = build_model( + hierarchical_specs, + noncentered=False, + ) + _component_dict_model, _component_messages = build_model( + hierarchical_specs, + noncentered={"v": False}, + ) + + parameterization_models = { + "scalar True": _scalar_noncentered_model, + "scalar False": _scalar_centered_model, + "dict: v=False; a omitted": _component_dict_model, + } + assert not (_scalar_nc_messages or _scalar_c_messages or _component_messages) + for _model in parameterization_models.values(): + assert_connected(_model) + return hierarchical_specs, parameterization_models + + +@app.cell +def _(group_term_structure, parameterization_models, pd): + _rows = [] + for _setting, _model in parameterization_models.items(): + for _parameter in ("v", "a"): + _row = group_term_structure(_model, _parameter, "1|participant_id") + _row = {"model setting": _setting, **_row} + _rows.append(_row) + + parameterization_table = pd.DataFrame(_rows)[ + [ + "model setting", + "parameter", + "effective form", + "direct group RV", + "offset RV", + "free mu RV", + "disconnected RVs", + ] + ] + assert parameterization_table.loc[ + parameterization_table["model setting"] == "scalar True", "offset RV" + ].all() + assert parameterization_table.loc[ + parameterization_table["model setting"] == "scalar False", "direct group RV" + ].all() + _dict_rows = parameterization_table[ + parameterization_table["model setting"] == "dict: v=False; a omitted" + ].set_index("parameter") + assert _dict_rows.loc["v", "effective form"] == "centered" + assert _dict_rows.loc["a", "effective form"] == "non-centered" + parameterization_table + return (parameterization_table,) + + +@app.cell +def _(model_graph, parameterization_models): + mixed_parameterization_graph = model_graph( + parameterization_models["dict: v=False; a omitted"] + ) + mixed_parameterization_graph + return (mixed_parameterization_graph,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + In the dictionary case, `v_1|participant_id` is sampled directly because + `v` is centered. The missing `a` key takes Bambi's default `True`, so + `a_1|participant_id_offset` is sampled and multiplied by its scale. Neither + graph contains a group `mu`: these are zero-mean deviations around their + matching common intercepts. + + The graph below fixes the mixed dictionary case as the same static view in + marimo, Molab, and the rendered documentation. + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## A per-prior override wins + + The finest control lives on the group prior itself. Here the component + dictionary requests centered `v`, while this one explicit prior requests + non-centering. Its `mu=0` and hierarchical `sigma` satisfy Bambi's current + non-centered contract. + """) + return + + +@app.cell +def _(assert_connected, build_model, free_rv_names, hssm): + zero_mean_noncentered_prior = hssm.Prior( + "Normal", + mu=0.0, + sigma=hssm.Prior("HalfNormal", sigma=0.5), + noncentered=True, + ) + _override_include = [ + { + "name": "v", + "formula": "v ~ 1 + (1 | participant_id)", + "prior": {"1|participant_id": zero_mean_noncentered_prior}, + } + ] + override_model, override_messages = build_model( + _override_include, + a=1.5, + noncentered={"v": False}, + ) + _override_free = free_rv_names(override_model) + assert not override_messages + assert "v_1|participant_id_offset" in _override_free + assert "v_1|participant_id_mu" not in _override_free + assert_connected(override_model) + return override_model, zero_mean_noncentered_prior + + +@app.cell +def _(group_term_structure, override_model, pd): + override_structure_table = pd.DataFrame( + [group_term_structure(override_model, "v", "1|participant_id")] + ) + override_structure_table + return (override_structure_table,) + + +@app.cell +def _(model_graph, override_model): + faithful_noncentered_graph = model_graph(override_model) + faithful_noncentered_graph + return (faithful_noncentered_graph,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + The graph contains the group scale and standard-normal offset, but no + `v_1|participant_id_mu`. This is not an orphan-removal trick: zero is the + intended location because the common `v_Intercept` owns the population mean. + + The complete resolution order is: + """) + return + + +@app.cell +def _(pd): + precedence_table = pd.DataFrame( + [ + { + "priority": 1, + "control": "per-prior noncentered", + "scope": "one explicit group term", + "rule": "wins when True or False", + }, + { + "priority": 2, + "control": "model-level component dictionary", + "scope": "all group terms for that HSSM parameter", + "rule": "named key wins; missing key defaults to True", + }, + { + "priority": 3, + "control": "model-level scalar", + "scope": "all group terms", + "rule": "True is the default; False requests centering", + }, + { + "priority": "safe-policy override", + "control": "generated unique group-only prior", + "scope": "only the generated location-owning term", + "rule": "HSSM centers it to preserve its location", + }, + ] + ) + precedence_table + return (precedence_table,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + The last row is deliberately separate from user precedence. HSSM may adapt + a prior that it generated itself, but an explicit user prior remains + authoritative and is never rewritten. + + ## Incompatible explicit non-centered locations fail before build + + A free group mean and a fixed nonzero group mean would both be omitted from + Bambi's `offset * sigma` construction. HSSM therefore rejects both before a + PyMC model exists. + """) + return + + +@app.cell +def _(expect_model_error, hssm, pd): + _free_mu_prior = hssm.Prior( + "Normal", + mu=hssm.Prior("Normal", mu=0.0, sigma=0.5), + sigma=hssm.Prior("HalfNormal", sigma=0.5), + ) + _nonzero_mu_prior = hssm.Prior( + "Normal", + mu=1.0, + sigma=hssm.Prior("HalfNormal", sigma=0.5), + ) + + _free_mu_include = [ + { + "name": "v", + "formula": "v ~ 1 + (1 | participant_id)", + "prior": {"1|participant_id": _free_mu_prior}, + } + ] + _nonzero_mu_include = [ + { + "name": "v", + "formula": "v ~ 1 + (1 | participant_id)", + "prior": {"1|participant_id": _nonzero_mu_prior}, + } + ] + free_mu_error = expect_model_error( + _free_mu_include, + a=1.5, + noncentered=True, + ) + nonzero_mu_error = expect_model_error( + _nonzero_mu_include, + a=1.5, + noncentered=True, + ) + + assert free_mu_error.startswith( + "Explicit group-specific prior specification(s) cannot be represented" + ) + assert "mu` hyperprior" in free_mu_error + assert "disconnected node" in free_mu_error + assert "not fixed entirely to zero" in nonzero_mu_error + assert "silently ignored" in nonzero_mu_error + + incompatible_location_table = pd.DataFrame( + [ + { + "explicit group mu": "Normal hyperprior (free)", + "effective setting": "non-centered", + "HSSM result": "ValueError before Bambi/PyMC build", + "reason": "Bambi would create mu but omit it from offset * sigma", + }, + { + "explicit group mu": "1.0 (fixed nonzero)", + "effective setting": "non-centered", + "HSSM result": "ValueError before Bambi/PyMC build", + "reason": "Bambi would ignore the requested location", + }, + ] + ) + incompatible_location_table + return free_mu_error, incompatible_location_table, nonzero_mu_error + + +@app.cell(hide_code=True) +def _(free_mu_error, mo): + mo.md(f""" + HSSM reports both the underlying limitation and term-specific repairs. The + free-mean case says: + + ```text + {free_mu_error} + ``` + + For a matched common/group expression, keep the common effect and use a + zero-mean group deviation. If the group distribution should own a free + population location, remove the matching common term and center that group + prior intentionally. + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## A unique group-only term owns its location + + Now `theta` appears only inside `(0 + theta | participant_id)`. The group + distribution must estimate the population slope: fixing its mean to zero + would change the scientific model. With `prior_settings="safe"`, HSSM keeps + that generated location and centers just this term—even when the model-level + request is non-centered. + """) + return + + +@app.cell +def _(assert_connected, build_model, free_rv_names): + _group_only_include = [ + { + "name": "v", + "formula": "v ~ 1 + (0 + theta | participant_id)", + } + ] + generated_owner_model, generated_owner_messages = build_model( + _group_only_include, + a=1.5, + noncentered=True, + ) + _owner_prior = generated_owner_model.params["v"].prior["theta|participant_id"] + _owner_free = free_rv_names(generated_owner_model) + assert _owner_prior.noncentered is False + assert "v_theta|participant_id" in _owner_free + assert "v_theta|participant_id_mu" in _owner_free + assert "v_theta|participant_id_offset" not in _owner_free + assert len(generated_owner_messages) == 1 + assert "generated location-bearing group-only term" in generated_owner_messages[0] + assert "Explicit priors were not changed" in generated_owner_messages[0] + assert_connected(generated_owner_model) + return generated_owner_messages, generated_owner_model + + +@app.cell(hide_code=True) +def _(generated_owner_messages, mo): + mo.md(f""" + HSSM makes the generated fallback visible: + + ```text + {generated_owner_messages[0]} + ``` + + In the graph, `v_theta|participant_id_mu` and the scale both feed the direct + group coefficient. There is no offset and no disconnected node. + """) + return + + +@app.cell +def _(generated_owner_model, group_term_structure, pd): + generated_owner_table = pd.DataFrame( + [group_term_structure(generated_owner_model, "v", "theta|participant_id")] + ) + generated_owner_table + return (generated_owner_table,) + + +@app.cell +def _(generated_owner_model, model_graph): + generated_owner_graph = model_graph(generated_owner_model) + generated_owner_graph + return (generated_owner_graph,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## Centered does not automatically mean identifiable + + Centering makes a free group `mu` part of the model, but the formula still + needs exactly one population-location owner. If common `theta` and the mean + of `theta|participant_id` are both free, the likelihood sees only their sum. + HSSM can build this model faithfully, so it warns about the location ridge + instead of rejecting the prior. + """) + return + + +@app.cell +def _(assert_connected, build_model, hssm): + _matched_free_location = hssm.Prior( + "Normal", + mu=hssm.Prior("Normal", mu=0.0, sigma=0.5), + sigma=hssm.Prior("HalfNormal", sigma=0.5), + ) + _matched_ridge_include = [ + { + "name": "v", + "formula": "v ~ 1 + theta + (0 + theta | participant_id)", + "prior": {"theta|participant_id": _matched_free_location}, + } + ] + matched_ridge_model, matched_ridge_messages = build_model( + _matched_ridge_include, + a=1.5, + noncentered=False, + ) + assert len(matched_ridge_messages) == 1 + assert "non-identifiable" in matched_ridge_messages[0] + assert "common 'theta'" in matched_ridge_messages[0] + assert "disconnected" not in matched_ridge_messages[0].lower() + assert_connected(matched_ridge_model) + return matched_ridge_messages, matched_ridge_model + + +@app.cell(hide_code=True) +def _(matched_ridge_messages, mo): + mo.md(f""" + ```text + {matched_ridge_messages[0]} + ``` + + The graph is connected, but connectivity is not identifiability. Both + `v_theta` and `v_theta|participant_id_mu` shift the same predictor. Keep + common `theta` and set the group `mu=0`, or remove common `theta` and let the + centered group distribution own the location. + """) + return + + +@app.cell +def _(group_term_structure, matched_ridge_model, pd): + matched_ridge_table = pd.DataFrame( + [group_term_structure(matched_ridge_model, "v", "theta|participant_id")] + ) + matched_ridge_table + return (matched_ridge_table,) + + +@app.cell +def _(matched_ridge_model, model_graph): + matched_ridge_graph = model_graph(matched_ridge_model) + matched_ridge_graph + return (matched_ridge_graph,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ### Repeated group-only owners have the same ridge + + A different ambiguity appears when the same unmatched expression has free + means under multiple grouping factors. Shifting every participant effect up + and every `conf` effect down leaves the predictor unchanged. Explicit priors + remain authoritative, so HSSM builds the centered model and emits one + aggregated warning. + """) + return + + +@app.cell +def _(assert_connected, build_model, group_term_structure, hssm, pd): + def _fresh_free_location(): + return hssm.Prior( + "Normal", + mu=hssm.Prior("Normal", mu=0.0, sigma=0.5), + sigma=hssm.Prior("HalfNormal", sigma=0.5), + ) + + _repeated_owner_include = [ + { + "name": "v", + "formula": ("v ~ 1 + (0 + theta | participant_id) + (0 + theta | conf)"), + "prior": { + "theta|participant_id": _fresh_free_location(), + "theta|conf": _fresh_free_location(), + }, + } + ] + repeated_owner_model, repeated_owner_messages = build_model( + _repeated_owner_include, + a=1.5, + noncentered=False, + ) + assert len(repeated_owner_messages) == 1 + assert "identified only by the priors" in repeated_owner_messages[0] + assert "theta|participant_id" in repeated_owner_messages[0] + assert "theta|conf" in repeated_owner_messages[0] + assert_connected(repeated_owner_model) + + repeated_owner_table = pd.DataFrame( + [ + group_term_structure( + repeated_owner_model, + "v", + _term, + ) + for _term in ("theta|participant_id", "theta|conf") + ] + ) + repeated_owner_table + return repeated_owner_messages, repeated_owner_model, repeated_owner_table + + +@app.cell(hide_code=True) +def _(mo, repeated_owner_messages): + mo.md(f""" + ```text + {repeated_owner_messages[0]} + ``` + + A cleaner formula adds common `theta` and uses mean-zero deviations for both + grouping factors. Alternatively, choose exactly one centered group term as + the location owner and fix the other group location intentionally. + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## Configuration errors fail fast too + + Component dictionaries are keyed by HSSM parameter names. A typo does not + silently fall back to a default: HSSM/Bambi list the valid names during + construction. + """) + return + + +@app.cell +def _(expect_model_error, pd): + _unknown_key_include = [{"name": "v", "formula": "v ~ 1 + (1 | participant_id)"}] + unknown_key_error = expect_model_error( + _unknown_key_include, + a=1.5, + noncentered={"vv": True}, + ) + assert "Unknown component name(s)" in unknown_key_error + assert "['a', 't', 'v', 'z']" in unknown_key_error + unknown_key_table = pd.DataFrame( + [ + { + "input": "noncentered={'vv': True}", + "result": "ValueError at construction", + "message": unknown_key_error, + } + ] + ) + unknown_key_table + return unknown_key_error, unknown_key_table + + +@app.cell(hide_code=True) +def _(mo): + mo.md(""" + ## Practical rules + + | Formula/prior situation | Recommended action | + | --- | --- | + | Common and group expressions match | Let the common term own the location; use group `mu=0`; choose centered or non-centered for sampling geometry. | + | One generated group-only expression | Let that group distribution own the location; HSSM's safe prior centers it automatically and explains the fallback. | + | Explicit group-only prior with a free location | Set that prior effectively centered so its `mu` is retained. | + | Explicit free or nonzero `mu` under non-centering | Change the formula/prior ownership or center the term; HSSM rejects the incompatible specification. | + | Two free centered owners for one expression | Add the exact common term and use zero-mean deviations, or choose exactly one owner. | + + Two final boundaries are worth remembering: + + - `noncentered` only affects group-specific terms. A component with no group + term is unchanged. + - HSSM's truncated group-prior wrapper is incompatible with Bambi's + hierarchical group-prior contract under either parameterization. Use an + untruncated hierarchical coefficient prior and a support-respecting link. + + Continue with: + + - [Specify hierarchical group priors](https://lnccbrown.github.io/HSSM/how_to/specify_group_priors/) for the exact compatibility and location-ownership rules; + - [Link functions and safe priors](https://lnccbrown.github.io/HSSM/tutorials/link_functions/) for predictor versus response scale; and + - [Centered vs. non-centered parameterizations](https://lnccbrown.github.io/HSSM/tutorials/centered_vs_noncentered_basic_logic/) for the underlying graph transformation. + + For the sampling geometry behind this choice, Michael Betancourt's + [Hierarchical Modeling](https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html) + case study provides a detailed treatment. + """) + return + + +if __name__ == "__main__": + app.run() From 4b26477a93af07c6faf4bd51f7eb10d6473d851f Mon Sep 17 00:00:00 2001 From: Alexander Fengler Date: Sun, 30 Aug 2026 17:11:26 -0400 Subject: [PATCH 3/3] docs: publish marimo parameterization sources (#1273) --- mkdocs.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mkdocs.yml b/mkdocs.yml index a7ad85480..c9d0ebd47 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -13,9 +13,11 @@ validation: # phase-2 consolidation targets, deliberately unlisted until then not_in_nav: | tutorials/attentional_ddm.py + tutorials/centered_vs_noncentered_basic_logic.py tutorials/main_tutorial.py tutorials/main_tutorial_scenic_route.py tutorials/link_functions.py + tutorials/parameterization_per_parameter.py tutorials/poisson_race.py tutorials/random_slope_safe_priors.py @@ -139,9 +141,11 @@ plugins: - mkdocs-jupyter: ignore: - tutorials/attentional_ddm.py + - tutorials/centered_vs_noncentered_basic_logic.py - tutorials/main_tutorial.py - tutorials/main_tutorial_scenic_route.py - tutorials/link_functions.py + - tutorials/parameterization_per_parameter.py - tutorials/poisson_race.py - tutorials/random_slope_safe_priors.py execute: true