Hierarchical Implementation: The h_m01 Stan Model

Jeff Helzner

Hierarchical Implementation: The h_m01 Stan Model

Foundational Report 9

foundations

implementation

h_m01

Concrete implementation of the hierarchical SEU sensitivity model in Stan. Details the data structure, non-centered parameterization, cell-specific belief formation, and generated quantities for model checking.

Author

Jeff Helzner

Published

May 12, 2026

0.1 Introduction

Report 8 established the abstract formulation of the hierarchical SEU sensitivity model: a regression on $\log(\alpha_j)$ links cell-level sensitivity to experimental factors, with cell-specific beliefs and shared utilities. Report 2 walked through the base m_01 implementation in Stan.

This report bridges the two by presenting the complete h_m01 Stan implementation block by block. We follow the same expository pattern as Report 2: each code block is presented as a listing, annotated with a table mapping mathematical symbols to Stan variables, and followed by interpretive commentary.

Reading This Report

This report assumes familiarity with the base model implementation (Report 2) and the hierarchical formulation (Report 8). We focus on what is new in h_m01 relative to m_01, noting unchanged components briefly.

0.2 Data Structure

The h_m01 data block extends m_01 with cell-level indexing and a design matrix:

h_m01.stan (data block)

data {
  // --- Dimensions ---
  int<lower=1> J;                         // number of experimental cells
  int<lower=2> K;                         // number of consequences
  int<lower=1> D;                         // embedding dimensions per alternative
  int<lower=2> R;                         // number of distinct alternatives (shared pool)
  int<lower=1> P;                         // number of predictors in design matrix (excluding intercept)

  // --- Shared alternatives ---
  array[R] vector[D] w;                   // feature vectors for alternatives (shared across cells)

  // --- Stacked observations ---
  int<lower=1> M_total;                   // total observations across all cells
  array[M_total] int<lower=1,upper=J> cell; // cell membership for each observation
  array[M_total, R] int<lower=0,upper=1> I; // indicator: I[m,r]=1 if alt r available in obs m
  array[M_total] int<lower=1> y;          // observed choices (1-indexed within active set)

  // --- Cell-level design matrix ---
  matrix[J, P] X;                         // predictor matrix (centered/coded), no intercept column

  // --- Per-cell observation counts (for bookkeeping/validation) ---
  array[J] int<lower=1> M_per_cell;       // M_per_cell[j] = number of obs in cell j
}

Notation: Hierarchical Data Structure

Symbol	Stan Variable	New?	Description
$J$	`J`	✓	Number of experimental cells
$P$	`P`	✓	Number of predictors in design matrix
$M_{\text{total}}$	`M_total`	✓	Total observations across all cells
$\text{cell}[m]$	`cell[m]`	✓	Cell membership for observation $m$
$M_j$	`M_per_cell[j]`	✓	Observation count in cell $j$
$\mathbf{X}$	`X[J, P]`	✓	Design matrix (no intercept)
$K, D, R$	`K`, `D`, `R`	✗	Unchanged from m_01
$\mathbf{w}_r$	`w[r]`	✗	Feature vectors (shared across cells)
$I_{m,r}$	`I[m, r]`	✗	Availability indicators (now stacked)
$y_m$	`y[m]`	✗	Observed choices (now stacked)

0.3 Transformed Data: Validation and Flattening

The transformed data block performs two functions: data validation and feature flattening.

h_m01.stan (transformed data block)

transformed data {
  // Validate stacked structure
  {
    array[J] int cell_count = rep_array(0, J);
    for (m in 1:M_total) {
      cell_count[cell[m]] += 1;
    }
    for (j in 1:J) {
      if (cell_count[j] != M_per_cell[j])
        reject("cell_count[", j, "] = ", cell_count[j],
               " but M_per_cell[", j, "] = ", M_per_cell[j]);
    }
  }

  // Calculate number of alternatives per observation
  array[M_total] int<lower=2> N_obs;
  int total_alts = 0;
  for (m in 1:M_total) {
    N_obs[m] = sum(I[m]);
    total_alts += N_obs[m];
    if (y[m] > N_obs[m])
      reject("y[", m, "] = ", y[m], " must be <= N_obs[", m, "] = ", N_obs[m]);
  }

  // Flatten feature vectors based on I (same pattern as m_01)
  array[total_alts] vector[D] x_flat;
  {
    int pos = 1;
    for (m in 1:M_total) {
      for (r in 1:R) {
        if (I[m, r] == 1) {
          x_flat[pos] = w[r];
          pos += 1;
        }
      }
    }
  }
}

The cell-count validation is new: it verifies that the cell vector is consistent with M_per_cell, catching data preparation errors at model instantiation rather than during sampling. The alternative-count computation (N_obs) and feature flattening (x_flat) follow the same pattern as m_01, applied to the stacked data.

0.4 Parameters

The parameters block organizes parameters into three groups: regression (new), cell-specific (extended), and shared (unchanged).

h_m01.stan (parameters block)

parameters {
  // --- Regression on log(α) ---
  real gamma0;                            // intercept (grand mean of log-alpha)
  vector[P] gamma;                        // predictor coefficients
  real<lower=0> sigma_cell;               // residual cell-level SD on log scale
  vector[J] z_alpha;                      // non-centered cell deviations (standard normal)

  // --- Per-cell belief parameters ---
  array[J] matrix[K, D] beta;            // cell-specific feature-to-probability mappings

  // --- Shared utility ---
  simplex[K-1] delta;                     // utility increments (shared)
}

The regression parameters (gamma0, gamma, sigma_cell, z_alpha) are entirely new. The belief parameters beta are extended from a single matrix to an array of $J$ matrices, one per cell. The utility increments delta are unchanged.

Non-Centered Parameterization

The model parameterizes cell deviations as $z_j \sim \mathcal{N}(0, 1)$ and constructs $\log(\alpha_j) = \gamma_0 + \mathbf{X}_j \boldsymbol{\gamma} + \sigma_{\text{cell}} \cdot z_j$, rather than placing a normal prior directly on $\log(\alpha_j)$. These are mathematically equivalent but the non-centered form performs better when $\sigma_{\text{cell}}$ is small (weak hierarchy), because it avoids the “funnel” geometry in the $(\sigma_{\text{cell}}, z_j)$ joint posterior.

When $\sigma_{\text{cell}}$ is well-identified and not too small, the centered parameterization $\log(\alpha_j) \sim \mathcal{N}(\gamma_0 + \mathbf{X}_j \boldsymbol{\gamma}, \sigma_{\text{cell}})$ can be more efficient. The choice is purely computational — the induced prior on $\alpha_j$ is identical.

0.5 Transformed Parameters

The transformed parameters block constructs the quantities needed for the likelihood from the raw parameters.

h_m01.stan (transformed parameters block)

transformed parameters {
  // Cell-specific log-alpha via non-centered parameterization
  vector[J] log_alpha;
  for (j in 1:J) {
    log_alpha[j] = gamma0 + X[j] * gamma + sigma_cell * z_alpha[j];
  }
  vector<lower=0>[J] alpha = exp(log_alpha);

  // Shared ordered utilities
  ordered[K] upsilon = cumulative_sum(append_row(0, delta));

  // Compute choice probabilities for all observations
  vector[total_alts] eta;                 // expected utilities (flat)
  {
    int pos = 1;
    for (m in 1:M_total) {
      int j = cell[m];
      for (idx in 1:N_obs[m]) {
        // Subjective probabilities from cell j's beta
        vector[K] psi_i = softmax(beta[j] * x_flat[pos]);
        eta[pos] = dot_product(psi_i, upsilon);
        pos += 1;
      }
    }
  }
}

The parameter flow is:

Regression → cell-specific $\alpha_j$: The regression coefficients $(\gamma_0, \boldsymbol{\gamma})$ combine with the design matrix $\mathbf{X}$ and cell deviations $z_j$ to produce $\log(\alpha_j)$, then exponentiated.
Shared $\boldsymbol{\delta}$ → ordered utilities: cumulative_sum(append_row(0, delta)) constructs $\boldsymbol{\upsilon}$ — unchanged from m_01.
Cell-specific $\boldsymbol{\beta}_j$ → expected utilities: For each observation $m$ in cell $j$, the belief parameters $\boldsymbol{\beta}_j$ map features to subjective probabilities $\boldsymbol{\psi}$, which are then dotted with utilities to produce expected utilities $\eta$.

The key difference from m_01 is the indexing: beta[j] (cell $j$’s belief matrix) replaces the single beta, and alpha[cell[m]] replaces the scalar alpha.

Show code

fig, ax = plt.subplots(figsize=(11, 6))
ax.set_xlim(0, 11)
ax.set_ylim(0, 6.5)
ax.axis('off')

# Box drawing helper
def draw_box(ax, x, y, w, h, text, color='white', edgecolor='black', fontsize=9, bold=False):
    rect = plt.Rectangle((x, y), w, h, facecolor=color, edgecolor=edgecolor, linewidth=1.5, zorder=2)
    ax.add_patch(rect)
    weight = 'bold' if bold else 'normal'
    ax.text(x + w/2, y + h/2, text, ha='center', va='center', fontsize=fontsize, fontweight=weight, zorder=3)

def draw_arrow(ax, x1, y1, x2, y2, color='black'):
    ax.annotate('', xy=(x2, y2), xytext=(x1, y1),
                arrowprops=dict(arrowstyle='->', color=color, lw=1.5), zorder=1)

# --- Regression parameters (top left) ---
draw_box(ax, 0.5, 5.2, 1.5, 0.8, '$\\gamma_0$\nintercept', '#d4e6f1')
draw_box(ax, 2.3, 5.2, 1.5, 0.8, '$\\boldsymbol{\\gamma}$\ncoefficients', '#d4e6f1')
draw_box(ax, 4.1, 5.2, 1.5, 0.8, '$\\sigma_{\\mathrm{cell}}$\nresidual SD', '#d4e6f1')

# Design matrix and z (top right)
draw_box(ax, 0.5, 3.8, 1.5, 0.8, '$\\mathbf{X}$\ndesign matrix', '#fdebd0')
draw_box(ax, 4.1, 3.8, 1.5, 0.8, '$z_j$\ndeviations', '#d4e6f1')

# log(alpha) -> alpha (middle)
draw_box(ax, 1.8, 2.4, 2.2, 0.8, '$\\log(\\alpha_j)$ → $\\alpha_j$\ncell sensitivity', '#abebc6', bold=True)

# Arrows into log(alpha)
draw_arrow(ax, 1.25, 5.2, 2.2, 3.2)
draw_arrow(ax, 3.05, 5.2, 2.9, 3.2)
draw_arrow(ax, 4.85, 5.2, 4.85, 4.6)
draw_arrow(ax, 4.85, 3.8, 3.6, 3.1)
draw_arrow(ax, 1.25, 3.8, 2.2, 3.1)

# --- Belief parameters (right) ---
draw_box(ax, 6.5, 4.0, 2.0, 0.8, '$\\boldsymbol{\\beta}_j$\ncell beliefs', '#d4e6f1')
draw_box(ax, 6.5, 5.2, 2.0, 0.8, '$\\mathbf{w}_r$\nfeatures', '#fdebd0')

# psi (middle right)
draw_box(ax, 6.5, 2.4, 2.0, 0.8, '$\\boldsymbol{\\psi}_{r,j}$\nsubjective probs', '#abebc6')

# Arrows into psi
draw_arrow(ax, 7.5, 5.2, 7.5, 4.8)
draw_arrow(ax, 7.5, 4.0, 7.5, 3.2)

# --- Shared utility (far right) ---
draw_box(ax, 9.3, 4.0, 1.5, 0.8, '$\\boldsymbol{\\delta}$\nincrements', '#d4e6f1')
draw_box(ax, 9.3, 2.4, 1.5, 0.8, '$\\boldsymbol{\\upsilon}$\nutilities', '#abebc6')

# Arrow into upsilon
draw_arrow(ax, 10.05, 4.0, 10.05, 3.2)

# --- Expected utility (bottom center) ---
draw_box(ax, 4.0, 0.8, 2.5, 0.8, '$\\eta_{r,m}$\nexpected utilities', '#abebc6', bold=True)

# Arrows into eta
draw_arrow(ax, 7.5, 2.4, 5.8, 1.6)
draw_arrow(ax, 10.05, 2.4, 6.2, 1.5)

# --- Choice probability (bottom) ---
draw_box(ax, 7.5, 0.8, 3.0, 0.8, '$P(y_m | \\alpha_{\\mathrm{cell}[m]}, \\eta)$\nchoice probability', '#f9e79f', bold=True)

# Arrows into choice prob
draw_arrow(ax, 6.5, 0.8 + 0.4, 7.5, 0.8 + 0.4)
draw_arrow(ax, 2.9, 2.4, 4.3, 1.6)

# Labels
ax.text(2.2, 6.2, 'Regression (new)', fontsize=10, fontstyle='italic', color='#2874a6')
ax.text(6.5, 6.2, 'Beliefs (per-cell)', fontsize=10, fontstyle='italic', color='#2874a6')
ax.text(9.3, 6.2, 'Utility (shared)', fontsize=10, fontstyle='italic', color='#2874a6')

ax.set_title('Parameter Flow in h_m01', fontsize=13, pad=10)
plt.tight_layout()
plt.show()

Figure 1: Parameter flow in h_m01. Regression parameters (γ₀, γ, σ_cell) combine with design matrix X and cell deviations z to produce cell-specific α. Cell-specific β maps features to beliefs. Shared δ determines utilities. All three combine to produce choice probabilities.

0.6 Prior Distributions

The model block specifies priors and the likelihood.

h_m01.stan (model block — priors)

model {
  // --- Priors ---
  // Regression (tightened 2026-04-23 after parameter recovery showed soft
  // identifiability of alpha with heavy right tail under the original priors)
  gamma0 ~ normal(2.5, 0.5);             // prior on intercept: median alpha ~ 12, 97.5th pct ~ 33
  gamma ~ normal(0, 0.5);                // stronger shrinkage toward no effect
  sigma_cell ~ normal(0, 0.3);           // half-normal: modest cell variation (J=6 weakly identifies)
  z_alpha ~ std_normal();                // non-centered parameterization

  // Beliefs: per-cell, same prior as m_01
  for (j in 1:J) {
    to_vector(beta[j]) ~ std_normal();
  }

  // Utilities: shared
  delta ~ dirichlet(rep_vector(1, K-1));
  // ...
}

Each prior is motivated by its role in the model:

Table 1: Prior distributions in h_m01 and their m_01 analogues.

Parameter	Prior	m_01 Analogue	Rationale
$\gamma_0$	$\mathcal{N}(2.5, 0.5)$	$\alpha \sim \text{LogNormal}(3.0, 0.75)$	Tightened intercept centred so that the prior on $\alpha$ has median $\approx 12$ and 97.5th percentile $\approx 33$, controlling the heavy right tail seen under the original $\mathcal{N}(3.0, 1.0)$ prior
$\boldsymbol{\gamma}$	$\mathcal{N}(0, 0.5)$	— (new)	Stronger shrinkage toward no effect; a priori, experimental factors don’t shift $\alpha$
$\sigma_{\text{cell}}$	$\text{HalfNormal}(0, 0.3)$	— (new)	Modest cell-level variation; 0.3 on log scale $\approx$ factor of 1.35. The $J=6$ design weakly identifies $\sigma_{\text{cell}}$, motivating the tighter scale
$z_j$	$\mathcal{N}(0, 1)$	— (new)	Non-centered parameterization; computational device
$\boldsymbol{\beta}_j$	$\mathcal{N}(0, 1)$ per element	$\boldsymbol{\beta} \sim \mathcal{N}(0, 1)$ per element	Identical to m_01, applied independently per cell
$\boldsymbol{\delta}$	$\text{Dirichlet}(\mathbf{1})$	$\boldsymbol{\delta} \sim \text{Dirichlet}(\mathbf{1})$	Unchanged

Connection to m_01 Prior Calibration

The m_01 prior $\alpha \sim \text{LogNormal}(3.0, 0.75)$ was calibrated via prior predictive analysis (Report 3) to yield a prior-implied SEU-max rate of approximately 0.78. The h_m01 priors above are the tightened values that match the current models/h_m01.stan (commit on 2026-04-24). They were tightened from an earlier $\gamma_0 \sim \mathcal{N}(3.0, 1.0)$, $\boldsymbol{\gamma} \sim \mathcal{N}(0, 1.0)$, $\sigma_{\text{cell}} \sim \text{HalfNormal}(0, 0.5)$ specification after parameter recovery (Report 11) showed soft identifiability of $\alpha$ with a heavy right tail. The hierarchical prior predictive (Report 10), parameter recovery (Report 11), and SBC validation (Report 12) all use the tightened priors shown here.

0.7 Likelihood

The likelihood is structurally identical to m_01, with one change: alpha is indexed by cell[m].

h_m01.stan (model block — likelihood)

  // --- Likelihood ---
  {
    int pos = 1;
    for (m in 1:M_total) {
      vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]);
      y[m] ~ categorical(softmax(alpha[cell[m]] * problem_eta));
      pos += N_obs[m];
    }
  }

Compare to m_01:

m_01.stan (likelihood)

      y[m] ~ categorical(softmax(alpha * problem_eta));

The only change is alpha[cell[m]] replacing the scalar alpha. This minimal modification connects the hierarchical regression structure to the observation model: each observation is generated according to its cell’s sensitivity parameter.

0.8 Generated Quantities

The generated quantities block computes diagnostics for model checking and comparison.

h_m01.stan (generated quantities — log-likelihood and predictions)

generated quantities {
  // Log-likelihood per observation
  vector[M_total] log_lik;
  {
    int pos = 1;
    for (m in 1:M_total) {
      vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]);
      log_lik[m] = categorical_lpmf(y[m] | softmax(alpha[cell[m]] * problem_eta));
      pos += N_obs[m];
    }
  }

  // Posterior predictive samples
  array[M_total] int y_pred;
  {
    int pos = 1;
    for (m in 1:M_total) {
      vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]);
      y_pred[m] = categorical_rng(softmax(alpha[cell[m]] * problem_eta));
      pos += N_obs[m];
    }
  }

Four categories of generated quantities are produced:

Log-likelihood (log_lik): Pointwise, for LOO-CV via the loo package. Same structure as m_01 but with cell-indexed alpha.
Posterior predictive samples (y_pred): For posterior predictive checks. Same structure as m_01.
PPC statistics: Log-likelihood discrepancy, modal choice accuracy, sum of chosen probabilities — all computed globally across all cells, following the same pattern as m_01.
Per-cell log-likelihoods (log_lik_cell): New — aggregates pointwise log-likelihoods by cell.

h_m01.stan (generated quantities — per-cell log-likelihoods)

  // Per-cell log-likelihoods (for cell-level model comparison)
  vector[J] log_lik_cell = rep_vector(0, J);
  for (m in 1:M_total) {
    log_lik_cell[cell[m]] += log_lik[m];
  }
}

The log_lik_cell vector enables diagnosing cells where the model fits poorly — a form of cell-level residual analysis that is unique to the hierarchical setting.

0.9 The Simulation Model (`h_m01_sim.stan`)

The simulation model generates synthetic choice data from the h_m01 prior, for use in prior predictive analysis and parameter recovery studies. All parameter generation occurs in generated quantities with fixed_param=True (no sampling needed).

h_m01_sim.stan (generated quantities — parameter generation)

generated quantities {
  // Draw regression parameters
  real gamma0 = normal_rng(gamma0_mean, gamma0_sd);
  vector[P] gamma;
  for (p in 1:P) {
    gamma[p] = normal_rng(0, gamma_sd);
  }
  real<lower=0> sigma_cell = abs(normal_rng(0, sigma_cell_sd));

  // Draw cell-level alphas
  vector[J] log_alpha;
  vector[J] alpha;
  for (j in 1:J) {
    real z_j = normal_rng(0, 1);
    log_alpha[j] = gamma0 + X[j] * gamma + sigma_cell * z_j;
    alpha[j] = exp(log_alpha[j]);
  }

  // Draw per-cell betas
  array[J] matrix[K, D] beta;
  for (j in 1:J) {
    for (k in 1:K) {
      for (d in 1:D) {
        beta[j][k, d] = normal_rng(0, beta_sd);
      }
    }
  }

  // Draw shared utilities
  simplex[K-1] delta = dirichlet_rng(rep_vector(1.0, K-1));
  vector[K] upsilon = cumulative_sum(append_row(0, delta));

  // Compute expected utilities and generate choices
  array[M_total] int y;
  {
    int pos = 1;
    for (m in 1:M_total) {
      int j = cell[m];
      vector[N_obs[m]] problem_eta;
      for (idx in 1:N_obs[m]) {
        vector[K] psi_i = softmax(beta[j] * x_flat[pos]);
        problem_eta[idx] = dot_product(psi_i, upsilon);
        pos += 1;
      }
      y[m] = categorical_rng(softmax(alpha[j] * problem_eta));
    }
  }
}

Key design features:

Hyperparameter controls (gamma0_mean, gamma0_sd, gamma_sd, sigma_cell_sd, beta_sd) are passed as data, allowing the prior to be varied without recompiling the model.
The structure mirrors Report 2’s m_01_sim.stan, extended for the hierarchical regression and per-cell beliefs.

0.10 The SBC Model (`h_m01_sbc.stan`)

The SBC (simulation-based calibration) model (Talts et al. 2018) combines parameter generation, data generation, and model fitting in a single Stan program. True parameters are drawn in transformed data (so they are fixed per chain), synthetic data are generated from those parameters, and the model is fit in the model block. The generated quantities block produces rank statistics for calibration checking.

h_m01_sbc.stan (generated quantities — rank statistics)

generated quantities {
  // Copy generated data
  array[M_total] int y_ = y;

  vector[2 + P + J + (K - 1)] pars_;
  vector[2 + P + J + (K - 1)] ranks_;
  {
    int idx = 1;

    // gamma0
    pars_[idx] = gamma0_;
    ranks_[idx] = (gamma0 > gamma0_) ? 1 : 0;
    idx += 1;

    // gamma
    for (p in 1:P) {
      pars_[idx] = gamma_[p];
      ranks_[idx] = (gamma[p] > gamma_[p]) ? 1 : 0;
      idx += 1;
    }

    // sigma_cell
    pars_[idx] = sigma_cell_;
    ranks_[idx] = (sigma_cell > sigma_cell_) ? 1 : 0;
    idx += 1;

    // alpha (per cell)
    for (j in 1:J) {
      pars_[idx] = alpha_[j];
      ranks_[idx] = (alpha[j] > alpha_[j]) ? 1 : 0;
      idx += 1;
    }

    // delta
    for (k in 1:(K-1)) {
      pars_[idx] = delta_[k];
      ranks_[idx] = (delta[k] > delta_[k]) ? 1 : 0;
      idx += 1;
    }
  }

The parameter ordering in pars_/ranks_ is: $\gamma_0$, $\gamma_1, \ldots, \gamma_P$, $\sigma_{\text{cell}}$, $\alpha_1, \ldots, \alpha_J$, $\delta_1, \ldots, \delta_{K-1}$.

Design Decision: Which Parameters to Track

Per-cell $\boldsymbol{\beta}_j$ matrices are not tracked in SBC — with $J \times K \times D$ elements, rank histograms would be unreadable. $\boldsymbol{\beta}$ calibration can be spot-checked separately. The tracked parameters are those with direct scientific interpretation: the regression structure ($\gamma_0, \boldsymbol{\gamma}, \sigma_{\text{cell}}$), cell-level sensitivities ($\alpha_j$), and utilities ($\boldsymbol{\delta}$).

Note that the pars_/ranks_ vector size is computed inline as 2 + P + J + (K - 1) rather than stored in a local variable. This is a necessary workaround: Stan does not permit non-data variables in top-level size declarations.

0.11 The Study Design Utility

The HierarchicalStudyDesign class in utils/study_design_hierarchical.py generates stacked data structures for h_m01. The typical workflow is:

from utils.study_design_hierarchical import HierarchicalStudyDesign

design = HierarchicalStudyDesign(
    J=6, K=3, D=2, R=10, P=2,
    M_per_cell=20,
    X=np.array([[0,0], [0,1], [1,0], [1,1], [1,0], [0,1]])
)
design.generate()           # Generate feature vectors and indicator arrays
data = design.get_data_dict()  # Stan-compatible data dictionary
design.save("design.json")  # Persist for reproducibility

Key features:

Constructor accepts all design parameters: J, K, D, R, P, M_per_cell (scalar or per-cell list), design matrix X, and feature generation settings.
generate() creates shared feature vectors $\mathbf{w}_r$ and the stacked indicator array with cell membership vector.
get_data_dict() returns a dictionary matching the h_m01.stan data block, including default simulation hyperparameters (gamma0_mean=2.5, gamma0_sd=0.5, gamma_sd=0.5, sigma_cell_sd=0.3, beta_sd=1.0). These tightened values were calibrated via the prior-predictive, parameter-recovery, and SBC analyses in Reports 10–12.
save() / load() serialize the design to JSON with metadata (design name, timestamp, feature distribution parameters) for reproducibility.

0.12 Summary

The following table summarizes how each m_01 component maps to h_m01:

Table 2: Mapping from m_01 to h_m01 components.

m_01 Component	h_m01 Component	Change
`alpha` (scalar)	`alpha[J]` (vector)	Scalar → $J$-dimensional, via regression structure
`beta` (matrix)	`beta[J]` (array of matrices)	Single → per-cell
`delta` (simplex)	`delta` (simplex)	Unchanged
$M$ observations	$M_{\text{total}}$ stacked observations + `cell[m]`	Single group → stacked multi-cell
$\text{LogNormal}(3.0, 0.75)$ prior on $\alpha$	$\mathcal{N}(3.0, 1.0)$ on $\gamma_0$ + regression + half-normal on $\sigma_{\text{cell}}$	Prior on $\alpha$ mediated by regression

The implementation changes are minimal — the likelihood differs only in indexing alpha[cell[m]] rather than a scalar alpha — but the inferential consequences are substantial: formal regression coefficients, partial pooling, and principled shrinkage across experimental cells.

Connection to Earlier Reports

Report 1: Softmax choice properties hold within each cell
Report 2: Base m_01 implementation that h_m01 extends
Reports 5–7: Risky-choice extensions could also be hierarchicalized
Report 8: Abstract formulation implemented here

References

Talts, Sean, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. 2018. “Validating Bayesian Inference Algorithms with Simulation-Based Calibration.” arXiv Preprint arXiv:1804.06788.

Reuse

CC BY-SA 4.0

Citation

BibTeX citation:

@online{helzner2026,
  author = {Helzner, Jeff},
  title = {Hierarchical {Implementation:} {The} h\_m01 {Stan} {Model}},
  date = {2026-05-12},
  url = {https://jeffhelzner.github.io/seu-sensitivity/foundations/09_hierarchical_implementation.html},
  langid = {en}
}

For attribution, please cite this work as:

Helzner, Jeff. 2026. “Hierarchical Implementation: The h_m01 Stan Model.” SEU Sensitivity Project, May 12. https://jeffhelzner.github.io/seu-sensitivity/foundations/09_hierarchical_implementation.html.

--- title: "Hierarchical Implementation: The h_m01 Stan Model" subtitle: "Foundational Report 9" description: | Concrete implementation of the hierarchical SEU sensitivity model in Stan. Details the data structure, non-centered parameterization, cell-specific belief formation, and generated quantities for model checking. categories: [foundations, implementation, h_m01] execute: cache: true --- ```{python} #| label: setup #| include: false import sys import os sys.path.insert(0, os.path.join(os.getcwd(), '..')) import numpy as np import matplotlib.pyplot as plt from scipy.special import softmax ``` ## Introduction [Report 8](08_hierarchical_formulation.qmd) established the abstract formulation of the hierarchical SEU sensitivity model: a regression on $\log(\alpha_j)$ links cell-level sensitivity to experimental factors, with cell-specific beliefs and shared utilities. [Report 2](02_concrete_implementation.qmd) walked through the base `m_01` implementation in Stan. This report bridges the two by presenting the complete `h_m01` Stan implementation block by block. We follow the same expository pattern as Report 2: each code block is presented as a listing, annotated with a table mapping mathematical symbols to Stan variables, and followed by interpretive commentary. ::: {.callout-note} ## Reading This Report This report assumes familiarity with the base model implementation ([Report 2](02_concrete_implementation.qmd)) and the hierarchical formulation ([Report 8](08_hierarchical_formulation.qmd)). We focus on what is **new** in `h_m01` relative to `m_01`, noting unchanged components briefly. ::: ## Data Structure The `h_m01` data block extends `m_01` with cell-level indexing and a design matrix: ```{.stan filename="h_m01.stan (data block)"} data { // --- Dimensions --- int<lower=1> J; // number of experimental cells int<lower=2> K; // number of consequences int<lower=1> D; // embedding dimensions per alternative int<lower=2> R; // number of distinct alternatives (shared pool) int<lower=1> P; // number of predictors in design matrix (excluding intercept) // --- Shared alternatives --- array[R] vector[D] w; // feature vectors for alternatives (shared across cells) // --- Stacked observations --- int<lower=1> M_total; // total observations across all cells array[M_total] int<lower=1,upper=J> cell; // cell membership for each observation array[M_total, R] int<lower=0,upper=1> I; // indicator: I[m,r]=1 if alt r available in obs m array[M_total] int<lower=1> y; // observed choices (1-indexed within active set) // --- Cell-level design matrix --- matrix[J, P] X; // predictor matrix (centered/coded), no intercept column // --- Per-cell observation counts (for bookkeeping/validation) --- array[J] int<lower=1> M_per_cell; // M_per_cell[j] = number of obs in cell j } ``` ::: {.callout-note} ## Notation: Hierarchical Data Structure | Symbol | Stan Variable | New? | Description | |--------|---------------|:----:|-------------| | $J$ | `J` | ✓ | Number of experimental cells | | $P$ | `P` | ✓ | Number of predictors in design matrix | | $M_{\text{total}}$ | `M_total` | ✓ | Total observations across all cells | | $\text{cell}[m]$ | `cell[m]` | ✓ | Cell membership for observation $m$ | | $M_j$ | `M_per_cell[j]` | ✓ | Observation count in cell $j$ | | $\mathbf{X}$ | `X[J, P]` | ✓ | Design matrix (no intercept) | | $K, D, R$ | `K`, `D`, `R` | ✗ | Unchanged from m_01 | | $\mathbf{w}_r$ | `w[r]` | ✗ | Feature vectors (shared across cells) | | $I_{m,r}$ | `I[m, r]` | ✗ | Availability indicators (now stacked) | | $y_m$ | `y[m]` | ✗ | Observed choices (now stacked) | ::: ## Transformed Data: Validation and Flattening The `transformed data` block performs two functions: data validation and feature flattening. ```{.stan filename="h_m01.stan (transformed data block)"} transformed data { // Validate stacked structure { array[J] int cell_count = rep_array(0, J); for (m in 1:M_total) { cell_count[cell[m]] += 1; } for (j in 1:J) { if (cell_count[j] != M_per_cell[j]) reject("cell_count[", j, "] = ", cell_count[j], " but M_per_cell[", j, "] = ", M_per_cell[j]); } } // Calculate number of alternatives per observation array[M_total] int<lower=2> N_obs; int total_alts = 0; for (m in 1:M_total) { N_obs[m] = sum(I[m]); total_alts += N_obs[m]; if (y[m] > N_obs[m]) reject("y[", m, "] = ", y[m], " must be <= N_obs[", m, "] = ", N_obs[m]); } // Flatten feature vectors based on I (same pattern as m_01) array[total_alts] vector[D] x_flat; { int pos = 1; for (m in 1:M_total) { for (r in 1:R) { if (I[m, r] == 1) { x_flat[pos] = w[r]; pos += 1; } } } } } ``` The **cell-count validation** is new: it verifies that the `cell` vector is consistent with `M_per_cell`, catching data preparation errors at model instantiation rather than during sampling. The alternative-count computation (`N_obs`) and feature flattening (`x_flat`) follow the same pattern as `m_01`, applied to the stacked data. ## Parameters The `parameters` block organizes parameters into three groups: regression (new), cell-specific (extended), and shared (unchanged). ```{.stan filename="h_m01.stan (parameters block)"} parameters { // --- Regression on log(α) --- real gamma0; // intercept (grand mean of log-alpha) vector[P] gamma; // predictor coefficients real<lower=0> sigma_cell; // residual cell-level SD on log scale vector[J] z_alpha; // non-centered cell deviations (standard normal) // --- Per-cell belief parameters --- array[J] matrix[K, D] beta; // cell-specific feature-to-probability mappings // --- Shared utility --- simplex[K-1] delta; // utility increments (shared) } ``` The regression parameters (`gamma0`, `gamma`, `sigma_cell`, `z_alpha`) are entirely new. The belief parameters `beta` are extended from a single matrix to an **array of $J$ matrices**, one per cell. The utility increments `delta` are unchanged. ::: {.callout-important} ## Non-Centered Parameterization The model parameterizes cell deviations as $z_j \sim \mathcal{N}(0, 1)$ and constructs $\log(\alpha_j) = \gamma_0 + \mathbf{X}_j \boldsymbol{\gamma} + \sigma_{\text{cell}} \cdot z_j$, rather than placing a normal prior directly on $\log(\alpha_j)$. These are mathematically equivalent but the non-centered form performs better when $\sigma_{\text{cell}}$ is small (weak hierarchy), because it avoids the "funnel" geometry in the $(\sigma_{\text{cell}}, z_j)$ joint posterior. When $\sigma_{\text{cell}}$ is well-identified and not too small, the centered parameterization $\log(\alpha_j) \sim \mathcal{N}(\gamma_0 + \mathbf{X}_j \boldsymbol{\gamma}, \sigma_{\text{cell}})$ can be more efficient. The choice is purely computational — the induced prior on $\alpha_j$ is identical. ::: ## Transformed Parameters The `transformed parameters` block constructs the quantities needed for the likelihood from the raw parameters. ```{.stan filename="h_m01.stan (transformed parameters block)"} transformed parameters { // Cell-specific log-alpha via non-centered parameterization vector[J] log_alpha; for (j in 1:J) { log_alpha[j] = gamma0 + X[j] * gamma + sigma_cell * z_alpha[j]; } vector<lower=0>[J] alpha = exp(log_alpha); // Shared ordered utilities ordered[K] upsilon = cumulative_sum(append_row(0, delta)); // Compute choice probabilities for all observations vector[total_alts] eta; // expected utilities (flat) { int pos = 1; for (m in 1:M_total) { int j = cell[m]; for (idx in 1:N_obs[m]) { // Subjective probabilities from cell j's beta vector[K] psi_i = softmax(beta[j] * x_flat[pos]); eta[pos] = dot_product(psi_i, upsilon); pos += 1; } } } } ``` The parameter flow is: 1. **Regression → cell-specific $\alpha_j$**: The regression coefficients $(\gamma_0, \boldsymbol{\gamma})$ combine with the design matrix $\mathbf{X}$ and cell deviations $z_j$ to produce $\log(\alpha_j)$, then exponentiated. 2. **Shared $\boldsymbol{\delta}$ → ordered utilities**: `cumulative_sum(append_row(0, delta))` constructs $\boldsymbol{\upsilon}$ — unchanged from `m_01`. 3. **Cell-specific $\boldsymbol{\beta}_j$ → expected utilities**: For each observation $m$ in cell $j$, the belief parameters $\boldsymbol{\beta}_j$ map features to subjective probabilities $\boldsymbol{\psi}$, which are then dotted with utilities to produce expected utilities $\eta$. The key difference from `m_01` is the indexing: `beta[j]` (cell $j$'s belief matrix) replaces the single `beta`, and `alpha[cell[m]]` replaces the scalar `alpha`. ```{python} #| label: fig-parameter-flow #| fig-cap: "Parameter flow in h_m01. Regression parameters (γ₀, γ, σ_cell) combine with design matrix X and cell deviations z to produce cell-specific α. Cell-specific β maps features to beliefs. Shared δ determines utilities. All three combine to produce choice probabilities." fig, ax = plt.subplots(figsize=(11, 6)) ax.set_xlim(0, 11) ax.set_ylim(0, 6.5) ax.axis('off') # Box drawing helper def draw_box(ax, x, y, w, h, text, color='white', edgecolor='black', fontsize=9, bold=False): rect = plt.Rectangle((x, y), w, h, facecolor=color, edgecolor=edgecolor, linewidth=1.5, zorder=2) ax.add_patch(rect) weight = 'bold' if bold else 'normal' ax.text(x + w/2, y + h/2, text, ha='center', va='center', fontsize=fontsize, fontweight=weight, zorder=3) def draw_arrow(ax, x1, y1, x2, y2, color='black'): ax.annotate('', xy=(x2, y2), xytext=(x1, y1), arrowprops=dict(arrowstyle='->', color=color, lw=1.5), zorder=1) # --- Regression parameters (top left) --- draw_box(ax, 0.5, 5.2, 1.5, 0.8, '$\\gamma_0$\nintercept', '#d4e6f1') draw_box(ax, 2.3, 5.2, 1.5, 0.8, '$\\boldsymbol{\\gamma}$\ncoefficients', '#d4e6f1') draw_box(ax, 4.1, 5.2, 1.5, 0.8, '$\\sigma_{\\mathrm{cell}}$\nresidual SD', '#d4e6f1') # Design matrix and z (top right) draw_box(ax, 0.5, 3.8, 1.5, 0.8, '$\\mathbf{X}$\ndesign matrix', '#fdebd0') draw_box(ax, 4.1, 3.8, 1.5, 0.8, '$z_j$\ndeviations', '#d4e6f1') # log(alpha) -> alpha (middle) draw_box(ax, 1.8, 2.4, 2.2, 0.8, '$\\log(\\alpha_j)$ → $\\alpha_j$\ncell sensitivity', '#abebc6', bold=True) # Arrows into log(alpha) draw_arrow(ax, 1.25, 5.2, 2.2, 3.2) draw_arrow(ax, 3.05, 5.2, 2.9, 3.2) draw_arrow(ax, 4.85, 5.2, 4.85, 4.6) draw_arrow(ax, 4.85, 3.8, 3.6, 3.1) draw_arrow(ax, 1.25, 3.8, 2.2, 3.1) # --- Belief parameters (right) --- draw_box(ax, 6.5, 4.0, 2.0, 0.8, '$\\boldsymbol{\\beta}_j$\ncell beliefs', '#d4e6f1') draw_box(ax, 6.5, 5.2, 2.0, 0.8, '$\\mathbf{w}_r$\nfeatures', '#fdebd0') # psi (middle right) draw_box(ax, 6.5, 2.4, 2.0, 0.8, '$\\boldsymbol{\\psi}_{r,j}$\nsubjective probs', '#abebc6') # Arrows into psi draw_arrow(ax, 7.5, 5.2, 7.5, 4.8) draw_arrow(ax, 7.5, 4.0, 7.5, 3.2) # --- Shared utility (far right) --- draw_box(ax, 9.3, 4.0, 1.5, 0.8, '$\\boldsymbol{\\delta}$\nincrements', '#d4e6f1') draw_box(ax, 9.3, 2.4, 1.5, 0.8, '$\\boldsymbol{\\upsilon}$\nutilities', '#abebc6') # Arrow into upsilon draw_arrow(ax, 10.05, 4.0, 10.05, 3.2) # --- Expected utility (bottom center) --- draw_box(ax, 4.0, 0.8, 2.5, 0.8, '$\\eta_{r,m}$\nexpected utilities', '#abebc6', bold=True) # Arrows into eta draw_arrow(ax, 7.5, 2.4, 5.8, 1.6) draw_arrow(ax, 10.05, 2.4, 6.2, 1.5) # --- Choice probability (bottom) --- draw_box(ax, 7.5, 0.8, 3.0, 0.8, '$P(y_m | \\alpha_{\\mathrm{cell}[m]}, \\eta)$\nchoice probability', '#f9e79f', bold=True) # Arrows into choice prob draw_arrow(ax, 6.5, 0.8 + 0.4, 7.5, 0.8 + 0.4) draw_arrow(ax, 2.9, 2.4, 4.3, 1.6) # Labels ax.text(2.2, 6.2, 'Regression (new)', fontsize=10, fontstyle='italic', color='#2874a6') ax.text(6.5, 6.2, 'Beliefs (per-cell)', fontsize=10, fontstyle='italic', color='#2874a6') ax.text(9.3, 6.2, 'Utility (shared)', fontsize=10, fontstyle='italic', color='#2874a6') ax.set_title('Parameter Flow in h_m01', fontsize=13, pad=10) plt.tight_layout() plt.show() ``` ## Prior Distributions The `model` block specifies priors and the likelihood. ```{.stan filename="h_m01.stan (model block — priors)"} model { // --- Priors --- // Regression (tightened 2026-04-23 after parameter recovery showed soft // identifiability of alpha with heavy right tail under the original priors) gamma0 ~ normal(2.5, 0.5); // prior on intercept: median alpha ~ 12, 97.5th pct ~ 33 gamma ~ normal(0, 0.5); // stronger shrinkage toward no effect sigma_cell ~ normal(0, 0.3); // half-normal: modest cell variation (J=6 weakly identifies) z_alpha ~ std_normal(); // non-centered parameterization // Beliefs: per-cell, same prior as m_01 for (j in 1:J) { to_vector(beta[j]) ~ std_normal(); } // Utilities: shared delta ~ dirichlet(rep_vector(1, K-1)); // ... } ``` Each prior is motivated by its role in the model: | Parameter | Prior | m_01 Analogue | Rationale | |-----------|-------|---------------|-----------| | $\gamma_0$ | $\mathcal{N}(2.5, 0.5)$ | $\alpha \sim \text{LogNormal}(3.0, 0.75)$ | Tightened intercept centred so that the prior on $\alpha$ has median $\approx 12$ and 97.5th percentile $\approx 33$, controlling the heavy right tail seen under the original $\mathcal{N}(3.0, 1.0)$ prior | | $\boldsymbol{\gamma}$ | $\mathcal{N}(0, 0.5)$ | — (new) | Stronger shrinkage toward no effect; a priori, experimental factors don't shift $\alpha$ | | $\sigma_{\text{cell}}$ | $\text{HalfNormal}(0, 0.3)$ | — (new) | Modest cell-level variation; 0.3 on log scale $\approx$ factor of 1.35. The $J=6$ design weakly identifies $\sigma_{\text{cell}}$, motivating the tighter scale | | $z_j$ | $\mathcal{N}(0, 1)$ | — (new) | Non-centered parameterization; computational device | | $\boldsymbol{\beta}_j$ | $\mathcal{N}(0, 1)$ per element | $\boldsymbol{\beta} \sim \mathcal{N}(0, 1)$ per element | Identical to m_01, applied independently per cell | | $\boldsymbol{\delta}$ | $\text{Dirichlet}(\mathbf{1})$ | $\boldsymbol{\delta} \sim \text{Dirichlet}(\mathbf{1})$ | Unchanged | : Prior distributions in h_m01 and their m_01 analogues. {#tbl-priors} ::: {.callout-note} ## Connection to m_01 Prior Calibration The m_01 prior $\alpha \sim \text{LogNormal}(3.0, 0.75)$ was calibrated via prior predictive analysis ([Report 3](03_prior_analysis.qmd)) to yield a prior-implied SEU-max rate of approximately 0.78. The h_m01 priors above are the **tightened** values that match the current `models/h_m01.stan` (commit on 2026-04-24). They were tightened from an earlier $\gamma_0 \sim \mathcal{N}(3.0, 1.0)$, $\boldsymbol{\gamma} \sim \mathcal{N}(0, 1.0)$, $\sigma_{\text{cell}} \sim \text{HalfNormal}(0, 0.5)$ specification after parameter recovery ([Report 11](11_hierarchical_parameter_recovery.qmd)) showed soft identifiability of $\alpha$ with a heavy right tail. The hierarchical prior predictive ([Report 10](10_hierarchical_prior_analysis.qmd)), parameter recovery ([Report 11](11_hierarchical_parameter_recovery.qmd)), and SBC validation ([Report 12](12_hierarchical_sbc_validation.qmd)) all use the tightened priors shown here. ::: ## Likelihood The likelihood is structurally identical to `m_01`, with one change: `alpha` is indexed by `cell[m]`. ```{.stan filename="h_m01.stan (model block — likelihood)"} // --- Likelihood --- { int pos = 1; for (m in 1:M_total) { vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]); y[m] ~ categorical(softmax(alpha[cell[m]] * problem_eta)); pos += N_obs[m]; } } ``` Compare to `m_01`: ```{.stan filename="m_01.stan (likelihood)"} y[m] ~ categorical(softmax(alpha * problem_eta)); ``` The only change is `alpha[cell[m]]` replacing the scalar `alpha`. This minimal modification connects the hierarchical regression structure to the observation model: each observation is generated according to its cell's sensitivity parameter. ## Generated Quantities The `generated quantities` block computes diagnostics for model checking and comparison. ```{.stan filename="h_m01.stan (generated quantities — log-likelihood and predictions)"} generated quantities { // Log-likelihood per observation vector[M_total] log_lik; { int pos = 1; for (m in 1:M_total) { vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]); log_lik[m] = categorical_lpmf(y[m] | softmax(alpha[cell[m]] * problem_eta)); pos += N_obs[m]; } } // Posterior predictive samples array[M_total] int y_pred; { int pos = 1; for (m in 1:M_total) { vector[N_obs[m]] problem_eta = segment(eta, pos, N_obs[m]); y_pred[m] = categorical_rng(softmax(alpha[cell[m]] * problem_eta)); pos += N_obs[m]; } } ``` Four categories of generated quantities are produced: 1. **Log-likelihood** (`log_lik`): Pointwise, for LOO-CV via the `loo` package. Same structure as `m_01` but with cell-indexed alpha. 2. **Posterior predictive samples** (`y_pred`): For posterior predictive checks. Same structure as `m_01`. 3. **PPC statistics**: Log-likelihood discrepancy, modal choice accuracy, sum of chosen probabilities — all computed globally across all cells, following the same pattern as `m_01`. 4. **Per-cell log-likelihoods** (`log_lik_cell`): **New** — aggregates pointwise log-likelihoods by cell. ```{.stan filename="h_m01.stan (generated quantities — per-cell log-likelihoods)"} // Per-cell log-likelihoods (for cell-level model comparison) vector[J] log_lik_cell = rep_vector(0, J); for (m in 1:M_total) { log_lik_cell[cell[m]] += log_lik[m]; } } ``` The `log_lik_cell` vector enables diagnosing cells where the model fits poorly — a form of cell-level residual analysis that is unique to the hierarchical setting. ## The Simulation Model (`h_m01_sim.stan`) {#sec-sim-model} The simulation model generates synthetic choice data from the `h_m01` prior, for use in prior predictive analysis and parameter recovery studies. All parameter generation occurs in `generated quantities` with `fixed_param=True` (no sampling needed). ```{.stan filename="h_m01_sim.stan (generated quantities — parameter generation)"} generated quantities { // Draw regression parameters real gamma0 = normal_rng(gamma0_mean, gamma0_sd); vector[P] gamma; for (p in 1:P) { gamma[p] = normal_rng(0, gamma_sd); } real<lower=0> sigma_cell = abs(normal_rng(0, sigma_cell_sd)); // Draw cell-level alphas vector[J] log_alpha; vector[J] alpha; for (j in 1:J) { real z_j = normal_rng(0, 1); log_alpha[j] = gamma0 + X[j] * gamma + sigma_cell * z_j; alpha[j] = exp(log_alpha[j]); } // Draw per-cell betas array[J] matrix[K, D] beta; for (j in 1:J) { for (k in 1:K) { for (d in 1:D) { beta[j][k, d] = normal_rng(0, beta_sd); } } } // Draw shared utilities simplex[K-1] delta = dirichlet_rng(rep_vector(1.0, K-1)); vector[K] upsilon = cumulative_sum(append_row(0, delta)); // Compute expected utilities and generate choices array[M_total] int y; { int pos = 1; for (m in 1:M_total) { int j = cell[m]; vector[N_obs[m]] problem_eta; for (idx in 1:N_obs[m]) { vector[K] psi_i = softmax(beta[j] * x_flat[pos]); problem_eta[idx] = dot_product(psi_i, upsilon); pos += 1; } y[m] = categorical_rng(softmax(alpha[j] * problem_eta)); } } } ``` Key design features: - **Hyperparameter controls** (`gamma0_mean`, `gamma0_sd`, `gamma_sd`, `sigma_cell_sd`, `beta_sd`) are passed as data, allowing the prior to be varied without recompiling the model. - The structure mirrors [Report 2](02_concrete_implementation.qmd)'s `m_01_sim.stan`, extended for the hierarchical regression and per-cell beliefs. ## The SBC Model (`h_m01_sbc.stan`) {#sec-sbc-model} The SBC (simulation-based calibration) model [@talts2018] combines parameter generation, data generation, and model fitting in a single Stan program. True parameters are drawn in `transformed data` (so they are fixed per chain), synthetic data are generated from those parameters, and the model is fit in the `model` block. The `generated quantities` block produces rank statistics for calibration checking. ```{.stan filename="h_m01_sbc.stan (generated quantities — rank statistics)"} generated quantities { // Copy generated data array[M_total] int y_ = y; vector[2 + P + J + (K - 1)] pars_; vector[2 + P + J + (K - 1)] ranks_; { int idx = 1; // gamma0 pars_[idx] = gamma0_; ranks_[idx] = (gamma0 > gamma0_) ? 1 : 0; idx += 1; // gamma for (p in 1:P) { pars_[idx] = gamma_[p]; ranks_[idx] = (gamma[p] > gamma_[p]) ? 1 : 0; idx += 1; } // sigma_cell pars_[idx] = sigma_cell_; ranks_[idx] = (sigma_cell > sigma_cell_) ? 1 : 0; idx += 1; // alpha (per cell) for (j in 1:J) { pars_[idx] = alpha_[j]; ranks_[idx] = (alpha[j] > alpha_[j]) ? 1 : 0; idx += 1; } // delta for (k in 1:(K-1)) { pars_[idx] = delta_[k]; ranks_[idx] = (delta[k] > delta_[k]) ? 1 : 0; idx += 1; } } ``` The parameter ordering in `pars_`/`ranks_` is: $\gamma_0$, $\gamma_1, \ldots, \gamma_P$, $\sigma_{\text{cell}}$, $\alpha_1, \ldots, \alpha_J$, $\delta_1, \ldots, \delta_{K-1}$. ::: {.callout-note} ## Design Decision: Which Parameters to Track Per-cell $\boldsymbol{\beta}_j$ matrices are **not** tracked in SBC — with $J \times K \times D$ elements, rank histograms would be unreadable. $\boldsymbol{\beta}$ calibration can be spot-checked separately. The tracked parameters are those with direct scientific interpretation: the regression structure ($\gamma_0, \boldsymbol{\gamma}, \sigma_{\text{cell}}$), cell-level sensitivities ($\alpha_j$), and utilities ($\boldsymbol{\delta}$). ::: Note that the `pars_`/`ranks_` vector size is computed inline as `2 + P + J + (K - 1)` rather than stored in a local variable. This is a necessary workaround: Stan does not permit non-data variables in top-level size declarations. ## The Study Design Utility {#sec-study-design} The `HierarchicalStudyDesign` class in `utils/study_design_hierarchical.py` generates stacked data structures for `h_m01`. The typical workflow is: ```python from utils.study_design_hierarchical import HierarchicalStudyDesign design = HierarchicalStudyDesign( J=6, K=3, D=2, R=10, P=2, M_per_cell=20, X=np.array([[0,0], [0,1], [1,0], [1,1], [1,0], [0,1]]) ) design.generate() # Generate feature vectors and indicator arrays data = design.get_data_dict() # Stan-compatible data dictionary design.save("design.json") # Persist for reproducibility ``` Key features: - **Constructor** accepts all design parameters: `J`, `K`, `D`, `R`, `P`, `M_per_cell` (scalar or per-cell list), design matrix `X`, and feature generation settings. - **`generate()`** creates shared feature vectors $\mathbf{w}_r$ and the stacked indicator array with cell membership vector. - **`get_data_dict()`** returns a dictionary matching the `h_m01.stan` data block, including default simulation hyperparameters (`gamma0_mean=2.5`, `gamma0_sd=0.5`, `gamma_sd=0.5`, `sigma_cell_sd=0.3`, `beta_sd=1.0`). These tightened values were calibrated via the prior-predictive, parameter-recovery, and SBC analyses in [Reports 10–12](10_hierarchical_prior_analysis.qmd). - **`save()` / `load()`** serialize the design to JSON with metadata (design name, timestamp, feature distribution parameters) for reproducibility. ## Summary The following table summarizes how each `m_01` component maps to `h_m01`: | m_01 Component | h_m01 Component | Change | |----------------|-----------------|--------| | `alpha` (scalar) | `alpha[J]` (vector) | Scalar → $J$-dimensional, via regression structure | | `beta` (matrix) | `beta[J]` (array of matrices) | Single → per-cell | | `delta` (simplex) | `delta` (simplex) | Unchanged | | $M$ observations | $M_{\text{total}}$ stacked observations + `cell[m]` | Single group → stacked multi-cell | | $\text{LogNormal}(3.0, 0.75)$ prior on $\alpha$ | $\mathcal{N}(3.0, 1.0)$ on $\gamma_0$ + regression + half-normal on $\sigma_{\text{cell}}$ | Prior on $\alpha$ mediated by regression | : Mapping from m_01 to h_m01 components. {#tbl-component-mapping} The implementation changes are minimal — the likelihood differs only in indexing `alpha[cell[m]]` rather than a scalar `alpha` — but the inferential consequences are substantial: formal regression coefficients, partial pooling, and principled shrinkage across experimental cells. ::: {.callout-note} ## Connection to Earlier Reports - [Report 1](01_abstract_formulation.qmd): Softmax choice properties hold within each cell - [Report 2](02_concrete_implementation.qmd): Base `m_01` implementation that `h_m01` extends - [Reports 5–7](05_adding_risky_choices.qmd): Risky-choice extensions could also be hierarchicalized - [Report 8](08_hierarchical_formulation.qmd): Abstract formulation implemented here :::

Parameter	Prior	m_01 Analogue	Rationale
\(\gamma_0\)	\(\mathcal{N}(2.5, 0.5)\)	\(\alpha \sim \text{LogNormal}(3.0, 0.75)\)	Tightened intercept centred so that the prior on \(\alpha\) has median \(\approx 12\) and 97.5th percentile \(\approx 33\), controlling the heavy right tail seen under the original \(\mathcal{N}(3.0, 1.0)\) prior
\(\boldsymbol{\gamma}\)	\(\mathcal{N}(0, 0.5)\)	— (new)	Stronger shrinkage toward no effect; a priori, experimental factors don’t shift \(\alpha\)
\(\sigma_{\text{cell}}\)	\(\text{HalfNormal}(0, 0.3)\)	— (new)	Modest cell-level variation; 0.3 on log scale \(\approx\) factor of 1.35. The \(J=6\) design weakly identifies \(\sigma_{\text{cell}}\), motivating the tighter scale
\(z_j\)	\(\mathcal{N}(0, 1)\)	— (new)	Non-centered parameterization; computational device
\(\boldsymbol{\beta}_j\)	\(\mathcal{N}(0, 1)\) per element	\(\boldsymbol{\beta} \sim \mathcal{N}(0, 1)\) per element	Identical to m_01, applied independently per cell
\(\boldsymbol{\delta}\)	\(\text{Dirichlet}(\mathbf{1})\)	\(\boldsymbol{\delta} \sim \text{Dirichlet}(\mathbf{1})\)	Unchanged