SPICE Analyses

Two layers of downstream analysis sit on top of a fitted SpiceEstimator:

1. Generative benchmarking (weinhardt2026/utils/task.py) — simulate new behavior by running the fitted model through the task environment, for comparison against real data (analysis_generative_comparison.py, one command per study).
2. Cross-study analysis pipelines (weinhardt2026/analysis/) — model evaluation, morphing, coefficient-level statistics, clustering. These operate on a fitted SpiceEstimator (and optionally its generated behavior) and are shared across all studies in weinhardt2026/studies/.

See training.md for the model/training internals these analyses consume, and studies.md for how individual studies wire them together.

---

Generative Benchmarking (weinhardt2026/utils/task.py)

Generative benchmarking simulates new behavioral data by running a fitted model through the original task environment. This produces synthetic datasets that can be compared against the original human data — the basis for the generative-comparison and behavioral-clustering analyses below.

Architecture

task.py                          # Shared infrastructure
├── Env (base class)             # Abstract task environment
└── generate_behavior()          # Batched trial-by-trial generation loop

studies/<study>/benchmarking_<study>.py  # Per-study file
├── get_dataset()                # Load & split data
├── BenchmarkModel (nn.Module)   # Hand-coded cognitive model (e.g. GQLModel)
├── Environment<Study>(Env)      # Study-specific reward mechanics
└── generate_behavior()          # Thin wrapper → calls shared _generate_behavior

Env Base Class

All task environments subclass Env and implement batched reset() + step():

class Env:
    def __init__(self, n_actions: int, n_participants: int, n_blocks: int):
        ...

    @property
    def n_sessions(self) -> int:
        return self.n_participants * self.n_blocks

    def reset(self, block_ids: torch.Tensor, participant_ids: torch.Tensor = None) -> None:
        """Set up per-session environment state from dataset metadata."""
        ...

    def step(self, action: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
        """One trial for all sessions in parallel.
        Args:    action: (n_sessions,) integer action indices.
        Returns: (reward, terminated) — both (n_sessions,) tensors.
        """
        ...

Shared generate_behavior() Flow

1. Extract block_ids and participant_ids from dataset metadata (xs[:, 0, 0, -3] and xs[:, 0, 0, -1])
2. Call environment.reset(block_ids, participant_ids)
3. Resolve the model’s device (model.device for a SpiceEstimator, else its parameters); the model is called as model(obs, state)
4. Draw the first action uniformly at random (the model has seen no trial yet)
5. For each trial t in range(n_trials):
   • reward, _ = environment.step(action_idx) — environment gives reward for current action
   • Build observation: one-hot action + partial-feedback reward (NaN for unchosen) + metadata from original dataset
   • Forward pass: rnn(obs, state) for BaseModel, rnn(obs, state) otherwise
   • Normalize logits: 5D (E,B,T,W,A) → mean over ensemble → 4D → extract (B, A)
   • Sample next action: multinomial(softmax(logits))
6. Restore NaN padding for variable-length sessions (matching original dataset structure)

use_sindy (SPICE-EQ vs SPICE-RNN) is set by the caller. Wrappers that unwrap a SpiceEstimator (bruckner2025, weber2024) must not call BaseModel.eval() on it: its default use_sindy=True would silently turn every SPICE-RNN run into SPICE-EQ.

Replay what the data record. Where the CSV stores the outcome of every option on every trial, the environment replays it instead of simulating a reward process, so a model is tested on the trials the participant actually faced: ganesh2024a (reward_left/reward_right; the contrast–reward contingency changes per block and is not stored), eckstein2026 (payout_1..4, scaled like normalize_rewards), bruckner2025 (outcomes, change points) and weber2024 (laser trajectories). dezfouli2019, braun2018 and bustamante2023 simulate the task’s reward process. kolff2025 replays the partner’s acts and the timing of the focal ape’s acts and is closed-loop only on the modelled decision (groom or not; the non-grooming act is the real one or, where the real act was grooming, drawn from the observed non-grooming acts).

Result: a synthetic SpiceDataset with the same shape/metadata as the input, but actions (and derived rewards) drawn from the fitted model instead of the human participant. Feed this into analysis_generative_comparison.py or analysis_behavioral_clustering.py below.

Study-Specific Environment Pattern

class EnvironmentMyStudy(Env):
    REWARD_PROBS = torch.tensor([...])  # Task-specific reward structure

    def __init__(self, n_actions, n_participants, n_blocks):
        super().__init__(n_actions, n_participants, n_blocks)

    def reset(self, block_ids, participant_ids=None):
        # Map block IDs to per-session reward parameters
        self.session_reward_probs = self.REWARD_PROBS[block_ids]

    def step(self, action):
        # Sample rewards based on task mechanics
        probs = self.session_reward_probs[torch.arange(len(action)), action]
        reward = torch.bernoulli(probs)
        return reward, torch.zeros(len(action), dtype=torch.bool)

Study-Specific generate_behavior Wrapper

def generate_behavior(model, path_data=None, dataset=None, save_dataset=None):
    if dataset is None:
        dataset, _, _ = get_dataset(path_data=path_data)
    environment = EnvironmentMyStudy(
        n_actions=dataset.n_actions,
        n_participants=dataset.n_participants,
        n_blocks=N_BLOCKS,
    )
    return _generate_behavior(dataset=dataset, model=model, environment=environment, save_dataset=save_dataset)

---

Cross-Study Analysis Pipelines (weinhardt2026/analysis/)

All analyses below are study-agnostic functions that take a fitted SpiceEstimator (loaded from a .pkl checkpoint or passed in-memory) plus a SpiceDataset, and write CSVs/plots to an output_dir. Each study wires these into its own data paths from weinhardt2026/studies/<study>/<study>.py — there are no per-study analysis runner scripts; anything reusable belongs here instead.

Analyses that only read coefficients do not need the training data. weinhardt2026/utils/checkpoints.py loads a checkpoint on its own:

load_estimator(model_path, spice_class, spice_config, n_actions=2,
               n_participants=None, polynomial_degree=2, model_kwargs=None)

inferring ensemble_size and n_participants from the saved tensors. It also provides select_lowest_bic_checkpoint(hpscan_csv, params_dir), participant_label_map(data_path, label_col) (SPICE’s 0-based participant index → a label column) and log_trials_per_participant(data_path) (the data-volume covariate that presence regressions adjust for).

Model Evaluation — analysis_model_evaluation.py

Quantifies how well a fitted SPICE model (RNN and SINDy-equation form) fits the training data and predicts held-out behavior, and compares against a hand-coded benchmark model and/or a plain GRU baseline. Both splits are evaluated in one call and reported in one table.

analysis_model_evaluation(
    dataset_train, dataset_test, spice_model=None, benchmark_model=None, gru_model=None,
    spice_selection='best-on-train', verbose=False, output_dir=None, trial_filter=None,
    benchmark_group_model=False,
)

• Training split: per-trial likelihood, NLL, AIC and BIC. Test split: per-trial likelihood and NLL only — held-out likelihood already measures generalization, so an information criterion on it would charge for parsimony twice.
• Number of parameters per model: n_parameters (mean ± std over the training split’s (participant, experiment) groups; SPICE-EQ counts the active coefficients) and n_parameters (total), the parameter count of the joint model.
• AIC/BIC treat the whole training split as one model (joint_information_criteria): the NLL is summed over all N trials of all participants, and K = n_parameters (total) counts every participant’s own parameters — per-participant parameters summed over participants (SPICE-EQ: over its distinct (participant, experiment) coefficient sets; benchmark: count_parameters() × participants), shared parameters once (GRU, SPICE-RNN, a SPICE group model with n_participants=1, or the benchmark with benchmark_group_model=True). AIC = 2 NLL + 2K, BIC = 2 NLL + K log N.
• AIC/BIC are reported only for the interpretable models (Benchmark, SPICE-EQ); GRU and SPICE-RNN get NaN (printed as -).
• spice_selection sets how the SPICE ensemble members enter the scores (SPICE-EQ and SPICE-RNN separately):
  • 'best-on-train' (default, used by every study): each member is scored on its own; the member with the highest training likelihood is reported, on train and test.
  • 'averaged-model': each member is scored on its own (its own logits and coefficient count); the scores are averaged over members.
  • 'averaged-logit': the members’ logits are averaged at every time step and that averaged model is scored (SPICE-EQ counts the terms active in any member).
• When SPICE and the benchmark are both evaluated, the benchmark’s row reports a per-participant sign test against SPICE-EQ for each split (Train … / Test … columns, participant_win_counts).
• trial_filter lets you restrict evaluation to a subset of trials (applied to both splits, e.g. excluding “wait” responses in a DDM-style task).
• The printed table is transposed (models as columns), numbers rounded to 3 decimals (integers without decimals), NaN as -; the returned DataFrame and the CSV keep full precision.
• analysis_model_evaluation_mse — parallel pipeline for continuous-output (regression) models with the same dataset_train/dataset_test/spice_selection interface: MSE, RMSE, MAE, R² on both splits, AIC/BIC on the training split from a joint Gaussian likelihood (one MLE variance over all trials; same parameter counting; Benchmark and SPICE-EQ only); best-on-train selects the member with the lowest training MSE. Writes model_evaluation_mse.csv and a plot of the test metrics.
• Result: a DataFrame (saved as model_evaluation.csv in output_dir) with one row per model (Benchmark, GRU, SPICE-RNN, SPICE-EQ) and columns n_parameters, n_parameters (total), Train Trial Lik., Train NLL, AIC, BIC, Test Trial Lik., Test NLL (plus stds) — the standard model-comparison table used to argue SPICE fits at least as well as hand-crafted models.

One-Step-Ahead Diagnosis — analysis_sindy_onestepahead.py

Diagnoses why the fitted SINDy equations underperform the RNN at test time: per-step approximation error vs. autoregressive error accumulation.

analysis_sindy_onestepahead(dataset, spice_model) -> (summary_df, per_trial_df)

Compares three evaluation modes: SPICE-RNN (autoregressive, ensemble mean), SPICE-SINDy (autoregressive) (equations feed their own past predictions forward), SPICE-SINDy (one-step-ahead) (equations receive the RNN’s true state at each trial). If one-step-ahead ≈ RNN but autoregressive ≪ RNN → the equations are locally accurate but their errors compound over trials. If one-step-ahead ≪ RNN → the equations poorly approximate the RNN even locally.

Model Morphing — analysis_morphing.py

Traces how SINDy equation coefficients change along a continuous behavioral/embedding axis (e.g. from healthy to depressed, or low to high reward rate), rather than treating participant groups as discrete.

run_morphing(
    estimator, dataset, metric_values, n_steps=20, morphing_range_sd=1.0,
    n_pruning_rounds=20, pruning_threshold=0.05, save_dir=None,
) -> dict  # {'member_results': [...]}

1. _find_morphing_direction — finds the direction in each ensemble member’s participant-embedding space that best predicts metric_values (e.g. via LinearRegression).
2. _create_morphed_dataset / _create_morphed_estimator — steps n_steps points along that direction, refitting SINDy coefficients at each point with a fast ridge → prune → ridge → prune → … → ridge cycle (closed-form solves, no SGD — orders of magnitude faster than full retraining).
3. Aggregates across ensemble members (each has its own RNN/embedding space, so morphing directions are found independently per member) into mean ± SE coefficient curves.
4. get_morphed_coefficients(result) extracts the coefficient trajectories for plotting.

analysis_morphing(...) wraps all of the above and is the entry point a study calls:

analysis_morphing(
    estimator, dataset, metric_values, output_dir,
    environment=None,           # a weinhardt2026.utils.task.Env; enables validation
    n_steps=20, morphing_range_sd=1.0, save_dir=None,
    n_validation_runs=5, prefix='morphing',
)

It writes <prefix>_coefficients.npz and, when an environment is given, calls validate_morphing() — which generates behavior from the morphed model at each step and checks that the metric actually follows the axis (monotonicity count + step/metric correlation), writing <prefix>_validation.npz. The environment is the only task-specific piece, so it is passed in rather than imported.

• Result: for each SINDy term, a curve of its coefficient value as a function of position along the morphing axis — showing e.g. a nonlinear exploration term’s coefficient collapsing toward zero as depression severity increases. This is how SPICE reveals structural (not just parametric) individual differences.

Generative Comparison — analysis_generative_comparison.py

Simulates every model in the study’s task and compares the simulated behavior with the real participants’. Runnable per study, from the master script or from the study’s own file:

python weinhardt2026/analysis/analysis_generative_comparison.py dezfouli2019 [kolff2025 ...] [--n_repeats N] [--cached]
python weinhardt2026/analysis/analysis_generative_comparison.py all
python weinhardt2026/studies/dezfouli2019/analysis_generative.py

--cached re-scores and re-plots the saved generative_metrics.csv without simulating. In a study pipeline the generate_data flag calls run_generative_analysis(study).

Study-specific part — weinhardt2026/studies/<study>/analysis_generative.py defines METRIC_LABELS, N_REPEATS (runs per model: 10, or 1 where the models predict deterministically — bruckner2025, weber2024), load_dataset() (all blocks, built as the models were trained on it), load_models(dataset), generate(model, dataset) and compute_metrics(dataset) (per-session arrays). The master provides vectorized metric helpers (choices_and_rewards, session_mean, stays, choice_entropy, quarter_masks). For the registered studies, data and models come from the STUDIES registry of weinhardt2026/figures/figure_model_evaluation.py (as run.py built them, via load_registered_dataset / load_registered_models); kolff2025 builds its own (group models of kolff2025.py).

Models.

• SPICE: one model per grid-search seed (params/params/): the ensemble member whose SPICE-EQ has the highest training likelihood (select_spice_members; EnsembleMember restricts an estimator’s outputs to that member). The same member is simulated as SPICE-EQ and SPICE-RNN. Every (seed, member) likelihood is written to generative_spice_selection.csv.
• GRU and benchmark: single fits, simulated as many times as there are SPICE models, with other random runs, so every model is scored on the same number of models × runs. A list of fits (e.g. seeds of a benchmark grid search) can take the place of the single fit.

Scores — per model (instance), then mean and SD over instances; all = mean over the metrics:

• distance — Wasserstein distance between the simulated participant values (all runs) and the real ones, divided by the SD of the real values. Lower is better; it asks whether a model gets the population right.
• rho — Spearman correlation across participants between the real value and the model’s participant mean over runs: whether individual differences are recovered.

Outputs in studies/<study>/results/: generative_behavior.png/.pdf (one panel per metric; violins in the order GRU, Benchmark, Real, SPICE-EQ, SPICE-RNN over participants × models × runs; for Benchmark, Real and SPICE-EQ every participant is a dot, joined across the three; d and ρ above each violin and averaged over metrics in the title), generative_metrics.csv (per model, instance, participant, run), generative_scores.csv, generative_scores_instances.csv, generative_spice_selection.csv and behavioral_metrics_real.csv (per-participant real metrics, read by behavioral clustering and the bestbic analyses). With several studies, weinhardt2026/results/generative_scores_all_studies.csv holds the across-metric scores per study and one number per model across studies (study = all, mean over studies).

Coefficient Distributions — analysis_coefficients_distributions.py

Population-level view of the fitted SINDy coefficients: which terms are present, how consistent they are across the ensemble, and how they compare across experiments.

analysis_coefficients_distributions(
    spice_model=None, model_path=None, model_module=None, model_class=None,
    model_config=None, dataset=None, output_dir="analysis_coefficient_distributions",
    max_participants_strip=30, cluster_heatmap=True,
) -> (coeff_df, presence_df, ensemble_consistency_df)

Produces: ensemble-spread plots and CV heatmaps (compute_ensemble_consistency, plot_ensemble_spread, plot_ensemble_cv_heatmap) showing how stable each coefficient is across ensemble members; violin plots of coefficient distributions across participants (plot_coefficient_violins); presence-rate bar charts (plot_presence_rate_bar) — what fraction of participants retain each term after pruning; experiment-comparison plots (plot_experiment_comparison); and a sparsity heatmap (plot_sparsity_heatmap) of which terms are active for which participants.

Individual-Differences Regression — analysis_coefficients_individuals.py

Tests whether individual SINDy coefficients relate to an external criterion (diagnosis, questionnaire score, task performance).

analysis_coefficients_individuals(
    path_data, criterion, analysis,      # analysis: "disc" (discrete/odds-ratio) or "cont" (continuous)
    reference=None, spice_model=None, path_model=None, ...,
    output_dir=None,
)

• Discrete (run_discrete): logistic regression of term presence on group membership (e.g. patient vs. control), reported as odds ratios with forest plots (_plot_forest, _plot_odds_ratios) and per-group presence rates (_plot_presence_rates).
• Continuous (run_continuous): logistic regression of term presence on the standardized criterion, adjusted for log data volume, with profile-likelihood CIs, forest plots and fitted logistic curves; jonckheere_terpstra tests for monotonic trend across ordered groups.
• Result: statistical evidence for which mechanisms (equation terms) differ between groups or scale with a trait — the individual-level counterpart to the morphing analysis above.

Coefficient Clustering — analysis_coefficient_dendrogram.py

Groups the candidate terms within each submodule and draws one dendrogram per module, so a fitted model can be read as a handful of co-occurring groups rather than as dozens of separate terms.

analysis_coefficient_dendrogram(
    output_dir, spice_model=None, model_path=None, spice_class=None, spice_config=None,
    mode="coefficients",          # or "presence"
    metric="abs",                 # coefficients mode: 1-|r| or "signed" (1-r)
    linkage_method="average", distance_threshold=0.5,
    min_active_fraction=0.05, min_overlap=30, pairwise_complete=True,
    n_permutations=0, grouping_alpha=0.05, permutation_seed=0,
    split_signs=True,             # presence mode: split every term by sign
    n_bootstrap=0, bootstrap_seed=0, prefix=None,
)

Two modes, one per axis:

• structural (mode="presence") — do the terms switch on and off together across participants? Distance is 1 − phi, where phi is the correlation between two terms’ 0/1 presence columns. Phi measures co-occurrence beyond what the terms’ frequencies imply and counts shared absence as much as shared presence. (Jaccard is deliberately not used: two terms that are each present in most participants overlap heavily even when entirely independent, so Jaccard merges common terms rather than co-occurring ones.) Terms that exclude each other have phi < 0 and are pushed apart. This is the mode the beta analysis builds on.

  Signs are part of the structure. With split_signs=True (default), every term is first replaced by two pseudo-terms, term (+) and term (-) (split_by_sign): a pseudo-term is present in a participant who has the term with that sign. A term that points one way in some participants and the other way in others describes two different mechanisms — e.g. a choice bias that makes a participant repeat vs one that makes them switch — so a group like {1 (-), choice[t-3] (+)} says which terms switch on together and in which direction. The two signs of one term are never present together, so their phi is negative and they never group. A sign nobody carries is dropped; a rare one (below min_active_fraction) falls to the activity filter. The trade-off: splitting reveals effects that cancel out when both signs are pooled, but costs power for terms whose direction does not matter to the question, and makes rare terms rarer.

• parametric (mode="coefficients") — given both terms are present, do their magnitudes move together? A pruned term is stored as an exact zero, so correlating full columns would mix this with the structural question. pairwise_complete=True (default) correlates each pair only over participants in which both terms are present; pairs with fewer than min_overlap such participants are set to r = 0 and counted as untestable.

Significance-based grouping (n_permutations > 0). Instead of a fixed distance_threshold, every merge gets a permutation p-value (merge_significance):

1. Shuffle each term’s column across participants independently. This keeps how common every term is but destroys any link between terms.
2. Rebuild the tree from the shuffled data and record its tightest merge — the strongest chance pair anywhere in the module.
3. Repeat n_permutations times (10,000 resolves p-values down to 1e-4 in ~13 s per model).

A real merge’s p-value is the share of shuffled trees whose tightest merge was at least as tight, p = (1 + #{null ≤ h}) / (N + 1). Because every merge is compared against the strongest chance pair in the whole module, this controls the family-wise error rate per module — no further correction is needed within a module; it is conservative for looser merges higher up the tree. Merge heights only increase up the tree, so “p < α” is exactly “height below a critical distance”: the cut for α is the k-th smallest null tightest merge, k = ⌈α(N+1) − 1⌉ (the 500th / 100th / 10th of 10,000 for α = .05 / .01 / .001). The groups are the significant merges: the tree is cut just below the grouping_alpha critical distance. Cuts differ per module, since more terms (more pairs for luck to work with) and rarer terms both produce tighter chance merges.

In the figure, each merge — the bracket joining its two children — is coloured by significance tier (dark = p < .001, mid = p < .01, light = p < .05, grey = n.s.), and the three critical distances are drawn as dashed lines in the matching tier colours. Tiers are defined once (STAR_TIERS) and shared by stars, colours and cut lines.

Significant co-occurrence means the terms switch on together more than chance — not that participants usually have all of them. Read the per-merge all (share with every term present) next to any (share with at least one): a large gap marks a loose group.

With n_bootstrap > 0, participants are resampled with replacement and reclustered at the module’s cut, giving per-cluster Hennig clusterwise Jaccard (mean best-match overlap of cluster memberships across draws; ≥0.75 conventionally “stable”) plus the weakest member pair’s co-assignment probability. These are reproducibility measures, not p-values.

• Outputs: dendrogram_*.png (one panel per module), dendrogram_heatmaps_*.png (similarity, in leaf order), dendrogram_clusters_*.csv (term → group, with group all/any activity), with permutations dendrogram_merges_*.csv (one row per merge — every grouping level — with distance, p, stars, all, any), and with bootstrapping dendrogram_stability_*.csv + dendrogram_coassignment_*.png.

Structural Beta Effects — analysis_coefficient_betas.py

Does the presence of a term — or of a group of terms that switch on together — relate to a per-participant behavioral metric? SPICE’s claim is about structural individual differences, so this tests presence only, never coefficient magnitude.

analysis_coefficient_betas(
    data_path, output_dir, spice_model=None, model_path=None, ...,
    criterion_col="reward",
    criterion_type="continuous",  # or "discrete" (e.g. a diagnosis column)
    reference=None,               # discrete: compare every group against this one
    comparisons=None,             # discrete: explicit [(group, reference), ...] pairs
    participant_col="participant", criterion_aggregate="mean",
    min_active_fraction=0.05, n_permutations=10000, grouping_alpha=0.05,
    split_signs=True, invariant_bounds=(0.05, 0.95), alpha=0.05, prefix=None,
)

1. Groups come from the presence dendrogram with significance-based grouping and sign splitting (above). A group counts as present in a participant only when all of its terms are, with their signs. The grouping never looks at the criterion, so it is computed once on all participants and shared by every contrast.
2. Every node inside every group — the group, its subgroups, down to its single terms — plus every ungrouped term is tested: logistic regression of presence on the criterion. Nodes present in fewer than 5% or more than 95% of the contrast’s participants carry no structural variation and are marked invariant instead of tested. Fits that separate perfectly, fail to converge, or — for a group comparison — have a group in which nobody (or everybody) carries the node are marked not estimable: the odds ratio is infinite there, and statsmodels only flags complete separation, so the 2x2 cells are checked directly.
3. Benjamini-Hochberg runs once over all tested nodes. The family is fixed up front: counting only the nodes the rule below visits would make the family shrink exactly when a group looks significant, so groups that look good would face a milder correction.
4. Break-up rule (select_reported): each branch starts at its largest significant group. If that group’s effect is significant it is reported and its subgroups are not; if it is n.s., invariant or not estimable, it is broken into its two children, down to single terms. The reported nodes partition the terms. The rule only decides what is reported, never how strict the correction is.

Continuous vs discrete criteria. A continuous criterion (e.g. reward rate) is one contrast over all participants, predictor = z(criterion); beta is the log-odds of presence per SD. A discrete criterion (e.g. diagnosis) is compared pair by pair — every group against reference, or the pairs in comparisons — each contrast over the participants of its two groups only, predictor = 1 for the group and 0 for the reference; beta is the log odds ratio of presence in the group vs the reference. Each contrast is its own analysis, with its own FDR family and its own break-up, so a group can hold in one contrast and break up in another; as in analysis_coefficients_individuals, there is no correction across contrasts. A discrete criterion is read as each participant’s first value, whatever criterion_aggregate says.

Larger groups tend to break up for two reasons: members whose relation to the metric differs cancel out once all of them are required, and requiring every term shrinks the group toward its rarest member.

Metric choice. Term survival depends on data volume through pruning, so a criterion that scales with how much a participant played (e.g. rewards per block) produces spurious presence effects. The analysis does not adjust for trial count; it checks on every run — Pearson correlation with log trial count for a continuous criterion, Kruskal-Wallis of log trial count across groups for a discrete one — and warns when they are related. Use a metric that does not scale with trial count, such as reward rate (criterion_col="reward", criterion_aggregate="mean").

Every model here is cross-sectional across participants: a term whose presence predicts performance does not establish that adding it would change performance.

• Outputs: per contrast, betas_structural_<prefix>_<criterion>.{csv,png} for a continuous criterion and betas_structural_<prefix>_<criterion>_<group>_vs_<reference>.{csv,png} for each discrete contrast, so separate runs never overwrite each other. The CSV lists every node (contrast, group, parent, children, depth, all/any, beta, OR, p, q, status, reported, broken_up); the forest plot shows the reported, tested leaves (filled markers = significant after FDR).

Deeper Look at One Checkpoint — analysis_bestbic.py

Four analyses that read a single fitted checkpoint (typically the lowest-BIC one from the sparsity scan, via select_lowest_bic_checkpoint) and describe how its discovered equations vary across participants.

analysis_embedding_dynamics(estimator, dataset, output_dir, results_dir,
                            label_map=None, module_symbols=None, n_select=3)
analysis_metric_fingerprint(estimator, output_dir, results_dir, metrics_csv,
                            data_path, metric='avg_reward', label_map=None)
analysis_group_betas(estimator, output_dir, results_dir, data_path,
                     criterion='diag', reference='Control')
analysis_metric_betas(estimator, output_dir, results_dir, data_path,
                      metrics_csv, metric='avg_reward')

1. Embedding extremes — the participants furthest apart in participant-embedding space, with their choice / value-update dynamics and their equations rendered as LaTeX (format_equations_latex; module_symbols maps a module name to its LaTeX symbol, and modules it omits are printed by name).
2. Fingerprint — structural (presence) and parametric (magnitude) coefficient differences ordered by a behavioral metric, adjusted for log data volume.
3. Group betas — presence regression of every coefficient on a discrete grouping column, wrapping analysis_coefficients_individuals and dropping near-constant terms.
4. Metric betas — the same against a continuous participant-level metric, attaching it to the behavioural CSV as a participant-constant column first.

Archived — coefficient compression and mechanism-level differences

analysis_coefficient_compression.py, analysis_mechanism_individuals.py, analysis_concepts.py and analysis_coefficient_ties.py have moved to weinhardt2026/analysis/archive/, together with the core modules they wrapped (sindy_compression.py, sindy_concepts.py, sindy_ties.py, now under archive/core/). They compressed per-participant coefficients into a small number of “mechanisms” via per-module NMF and re-tested group differences at the mechanism level. analysis_coefficient_dendrogram.py above covers the same question — which coefficients belong together — without a factorization step, and is used for interpretation only.

Behavioral Clustering — analysis_behavioral_clustering.py

Tests whether clusters found in raw behavior (e.g. average reward, switch rate, from generate_behavior) align with structural differences in the fitted equations.

analysis_behavioral_clustering(
    spice_model, path_behavioral_metrics, n_clusters=3, output_dir='results',
) -> dict  # {'labels', 'linkage', 'centroids', 'nearest', 'equation_tests', 'alignment_ari', 'df_metrics'}

1. Loads per-participant behavioral metrics CSV (produced by a generative-behavior analysis).
2. Hierarchical clustering (linkage, fcluster, Ward’s method) on standardized behavioral metrics.
3. Extracts equation features (_extract_equation_features: coefficients + presence per participant).
4. Tests whether equation structure differs across behavioral clusters (_test_equation_differences: Kruskal-Wallis / Mann-Whitney) and reports adjusted_rand_score alignment between behavioral clusters and any independently-known grouping.
   • Result: evidence for (or against) the claim that behaviorally-defined subgroups correspond to structurally distinct equations, not just parameter shifts.

Reward-History Kernel — analysis_reward_history_kernel.py

A classic behavioral-analysis baseline, independent of any fitted SPICE model: logistic regression of “stay with previous action” on reward history.

compute_reward_history_kernel(dataset, max_lag=6) -> dict  # 'lags', 'coef_own', 'coef_other', 'se_own', 'se_other', + statsmodels result

For each trial, regresses stay = 1[action[t] == action[t-1]] jointly on reward history for the previously-chosen action (reward_own, lags 1..max_lag) and the other action (reward_other, lags 2..max_lag; lag-1 is structurally collinear with the intercept and dropped). Works on both real and model-generated datasets — run once on human data and once on generate_behavior output to check whether a fitted model reproduces the reward-history kernel shape.

Parameter Recovery — analysis_parameter_recovery.py

Synthetic-data sanity check: simulate a known ground-truth model (e.g. QLearning from weinhardt2026/studies/synthetic/), fit SPICE to the simulated data, and verify the recovered SINDy coefficients match the ground-truth parameters. Handles term collapsing for binary signals (signal^1 == signal^2, e.g. binary reward/choice) when comparing fitted vs. true coefficients.

Sparsity Hyperparameter Scan — analysis_sparsity_hpscan.py

analysis_sparsity_hpscan(
    pkl_pattern, spice_class, spice_config, n_actions, data_path, test_blocks,
    polynomial_degree=2, model_kwargs=None, device=None,
) -> pd.DataFrame

Evaluates a batch of checkpoints from a pruning-threshold × pruning-test hyperparameter sweep (glob-matched .pkl files, e.g. params_array/spice_dezfouli2019_*_*.pkl) on held-out test_blocks, returning a summary table of predictive performance vs. sparsity level — used to pick sindy_threshold_pruning / sindy_ensemble_pruning for a study.

---

Typical Analysis Sequence for a New Study

1. Fit SpiceEstimator on training data, then evaluate on the training and held-out data in one call with analysis_model_evaluation.py.
2. Run analysis_sindy_onestepahead.py if SINDy test performance lags the RNN, to see whether it’s a per-step or accumulation problem.
3. Inspect population-level equation structure with analysis_coefficients_distributions.py.
4. Relate coefficients to external criteria: analysis_coefficients_individuals.py (discrete groups or continuous traits), analysis_coefficient_betas.py for structural effects of co-occurring groups and single terms on a behavioral metric, and/or analysis_morphing.py for a continuous structural trajectory.
5. Write studies/<study>/analysis_generative.py, run python weinhardt2026/analysis/analysis_generative_comparison.py <study> to validate simulated against real behavior, and check the reward-history kernel with compute_reward_history_kernel.
6. Read the fitted equations as a small number of co-occurring groups with analysis_coefficient_dendrogram.py (mode="presence", n_permutations=10000), or check behavioral-cluster/equation alignment with analysis_behavioral_clustering.py.
7. For the selected checkpoint, run the four analysis_bestbic.py analyses (embedding extremes, fingerprint, group and metric betas) to describe how the equations vary across participants.