spectralbrain.statistics.bayesian#

Bayesian statistical models for spectral morphometry.

Six models with a scikit-learn-like API: .fit(), .predict(), .score(), .summary(). All models delegate MCMC sampling to the backends (PyMC NUTS, nutpie, or NumPyro).

Models#

HorseshoeRegression — sparse regression for feature selection.
BayesianGroupComparison — BEST (Kruschke 2013) with HDI + ROPE.
HierarchicalLinearModel — multi-site random effects.
GaussianProcessNormative — GP age-trajectory normative.
BayesianSpatialModel — GMRF vertex-wise spatial prior.
BayesianConnectome — hierarchical connectome comparison.

Examples

>>> model = HorseshoeRegression()
>>> model.fit(descriptors, clinical_scores)
>>> model.summary()
>>> predictions = model.predict(new_descriptors)
>>> model.score()  # LOO-CV

Dependencies#

PyMC, ArviZ (optional, lazy-imported).

Classes

`BayesianConnectome`([shrinkage])	Hierarchical Bayesian model for connectome comparison.
`BayesianGroupComparison`([rope])	Bayesian Estimation Supersedes the T-test (BEST).
`BayesianModel`()	Abstract base for all SpectralBrain Bayesian models.
`BayesianSpatialModel`([spatial_strength])	Vertex-wise Bayesian model with spatial GMRF prior.
`GaussianProcessNormative`([kernel, ...])	GP-based normative model for age trajectories.
`HierarchicalLinearModel`([random_effects])	Multi-site hierarchical linear model with random effects.
`HorseshoeRegression`([tau_prior])	Sparse Bayesian regression with horseshoe prior.

class spectralbrain.statistics.bayesian.BayesianConnectome(shrinkage=1.0)[source]#

Bases: BayesianModel

Hierarchical Bayesian model for connectome comparison.

Models each entry of the geometric connectome matrix with a hierarchical prior, testing whether connection strengths differ between groups while sharing information across edges.

Parameters:: shrinkage (float) – Hierarchical shrinkage strength.

Examples

>>> model = BayesianConnectome()
>>> model.fit(connectomes_patients, connectomes_controls)
>>> diff = model.edge_difference_posterior()

edge_difference_matrix()[source]#

Reconstruct the difference as a symmetric matrix.

Returns:: ndarray, shape (R, R)
Return type:: ndarray[tuple[Any, …], dtype[floating]]

edge_difference_posterior()[source]#

Posterior mean of per-edge group difference.

Returns:: ndarray, shape (n_edges,)
Return type:: ndarray

fit(group_a_connectomes, group_b_connectomes, **kwargs)[source]#

Fit connectome comparison model.

Parameters:

group_a_connectomes (ndarray, shape (n_a, R, R))
group_b_connectomes (ndarray, shape (n_b, R, R))

Return type:

BayesianConnectome

predict(X_new=None, **kwargs)[source]#: Not applicable for connectome comparison; use edge_difference_posterior().

class spectralbrain.statistics.bayesian.BayesianGroupComparison(rope=(-0.1, 0.1))[source]#

Bases: BayesianModel

Bayesian Estimation Supersedes the T-test (BEST).

Estimates the full posterior distribution of group means, standard deviations, effect size, and their differences. Reports HDI (Highest Density Interval) and probability of the difference exceeding a ROPE.

Parameters:: rope (tuple of float) – Region Of Practical Equivalence. Default (-0.1, 0.1) in units of the pooled standard deviation.

Examples

>>> model = BayesianGroupComparison(rope=(-0.1, 0.1))
>>> model.fit(group_a_values, group_b_values)
>>> model.summary()
>>> model.effect_size_posterior()

effect_size_posterior()[source]#

Posterior samples of Cohen’s d.

Returns:: ndarray, shape (n_samples,)
Return type:: ndarray

fit(group_a, group_b, **kwargs)[source]#

Fit BEST model.

Parameters:

group_a (ndarray, shape (n_a,) and (n_b,))
group_b (ndarray, shape (n_a,) and (n_b,))

Return type:

BayesianGroupComparison

predict(X_new=None, **kwargs)[source]#: Not applicable for group comparison; use effect_size_posterior().

rope_probability()[source]#

Probability of effect size in, below, and above ROPE.

Returns:: dict – Keys: "p_rope" (inside), "p_below" (below), "p_above" (above ROPE).
Return type:: dict[str, float]

class spectralbrain.statistics.bayesian.BayesianModel[source]#

Bases: ABC

Abstract base for all SpectralBrain Bayesian models.

Subclasses implement _build_model() to construct a PyMC model, and optionally override _build_posterior_predictive() for custom prediction logic.

trace_#

Posterior samples (populated after .fit()).

Type:: arviz.InferenceData or None

model_#

The PyMC model object.

Type:: pymc.Model or None

classmethod load_trace(path)[source]#

Load a saved trace.

Parameters:: path (str | PathLike)
Return type:: Any

fit(X, y, *, sampler='auto', draws=2000, tune=1000, chains=4, cores=4, target_accept=0.95, seed=42, **kwargs)[source]#

Fit the model via MCMC sampling.

Parameters:

X (ndarray, shape (n, d)) – Feature matrix.
y (ndarray, shape (n,)) – Target variable.
sampler (str) – "auto" tries nutpie → numpyro → nuts.
draws (int) – MCMC configuration.
tune (int) – MCMC configuration.
chains (int) – MCMC configuration.
cores (int) – MCMC configuration.
target_accept (float)
seed (int)
**kwargs – Extra arguments passed to the sampler.

Returns:

self

Return type:

BayesianModel

predict(X_new, *, n_samples=500, seed=None)[source]#

Generate posterior predictive samples.

Parameters:

X_new (ndarray, shape (m, d)) – New feature values.
n_samples (int) – Number of posterior predictive draws.
seed (int)

Returns:

ndarray, shape (n_samples, m) – Posterior predictive samples.

Return type:

ndarray

save(path)[source]#

Save trace to NetCDF.

Parameters:: path (str | PathLike)
Return type:: Path

score(method='loo')[source]#

Model comparison score via LOO-CV or WAIC.

Parameters:: method (str) – "loo" — Leave-One-Out via PSIS. "waic" — Widely Applicable Information Criterion.
Returns:: float – Expected log pointwise predictive density (elpd).
Return type:: float

summary(var_names=None, hdi_prob=0.94)[source]#

ArviZ summary table.

Parameters:

var_names (list of str, optional)
hdi_prob (float)

Returns:

pandas.DataFrame

Return type:

Any

class spectralbrain.statistics.bayesian.BayesianSpatialModel(spatial_strength=10.0)[source]#

Bases: BayesianModel

Vertex-wise Bayesian model with spatial GMRF prior.

Places a Gaussian Markov Random Field prior on the vertex-level effects, so neighbouring vertices share information. This is Bayesian spatial smoothing — more principled than Gaussian kernel pre-smoothing.

Parameters:: spatial_strength (float) – Precision of the GMRF prior (higher = more spatial smoothing).

Examples

>>> model = BayesianSpatialModel(spatial_strength=10.0)
>>> model.fit(group_labels, vertex_descriptors,
...           adjacency=mesh_adjacency)

fit(group_labels, vertex_data, *, adjacency, **kwargs)[source]#

Fit spatial model.

Parameters:

group_labels (ndarray, shape (S,)) – Group assignment (0/1) per subject.
vertex_data (ndarray, shape (S, N)) – Per-subject vertex-wise descriptor values.
adjacency (sparse matrix, shape (N, N)) – Mesh or kNN adjacency.

vertex_effect_map()[source]#

Posterior mean of vertex-level group effect.

Returns:: ndarray, shape (N,)
Return type:: ndarray

class spectralbrain.statistics.bayesian.GaussianProcessNormative(kernel='matern52', lengthscale_prior=10.0)[source]#

Bases: BayesianModel

GP-based normative model for age trajectories.

Fits a Gaussian Process over age (or any continuous covariate) to model the normative distribution of a spectral descriptor. Individual deviations are computed as z-scores from the posterior predictive.

Parameters:

kernel (str) – "matern32" or "matern52" or "rbf".
lengthscale_prior (float) – Prior mean for the GP lengthscale (in years for age).

Examples

>>> gp = GaussianProcessNormative(kernel="matern52")
>>> gp.fit(ages_controls[:, None], descriptor_controls)
>>> z_patient = gp.deviation(age_patient, descriptor_patient)

deviation(age, observed_value)[source]#

Z-score deviation of an individual from the normative.

Parameters:

age (float)
observed_value (float)

Returns:

float – Z-score (positive = above normative).

Return type:

float

predict(X_new, **kwargs)[source]#

Posterior predictive mean and std at new points.

Parameters:

X_new (ndarray, shape (m, d))

Returns:

mean (ndarray, shape (m,))
std (ndarray, shape (m,))

Return type:

tuple[ndarray, ndarray]

class spectralbrain.statistics.bayesian.HierarchicalLinearModel(random_effects='intercept')[source]#

Bases: BayesianModel

Multi-site hierarchical linear model with random effects.

Models spectral descriptors with fixed effects (group, age, sex) and random intercepts/slopes per site, handling batch effects within the model rather than post-hoc harmonisation.

y ~ α + β·X + u_site + ε

Parameters:: random_effects (str) – "intercept" — random intercept per site. "slope" — random intercept + slope per site.

Examples

>>> model = HierarchicalLinearModel(random_effects="intercept")
>>> model.fit(X, y, site_labels=sites)

fit(X, y, *, site_labels, **kwargs)[source]#

Fit with site labels for random effects.

Parameters:: site_labels (ndarray)

site_effects()[source]#

Posterior mean of random intercepts per site.

Returns:: ndarray, shape (n_sites,)
Return type:: ndarray

class spectralbrain.statistics.bayesian.HorseshoeRegression(tau_prior=0.1)[source]#

Bases: BayesianModel

Sparse Bayesian regression with horseshoe prior.

The horseshoe prior (Carvalho, Polson & Scott 2009) provides aggressive shrinkage of irrelevant coefficients toward zero while leaving large effects unshrunk — ideal for selecting which of 20+ spectral descriptors predict a clinical outcome.

Parameters:: tau_prior (float) – Global shrinkage scale (smaller = more sparse). Rule of thumb: p_eff / (d - p_eff) / sqrt(n), where p_eff is expected number of relevant features.

Examples

>>> model = HorseshoeRegression(tau_prior=0.1)
>>> model.fit(descriptors, clinical_scores)
>>> model.summary()
>>> model.predict(new_descriptors)

feature_importance()[source]#

Posterior mean of |β| — higher = more important.

Returns:: ndarray, shape (d,)
Return type:: ndarray