`confusius.glm`¶

glm ¶

GLM analysis for fUSI data.

This package provides General Linear Model tools for voxel-wise statistical analysis of functional ultrasound imaging data. The code is adapted from Nilearn's GLM module under the BSD-3-Clause license.

Modules:

first_level –

First-level GLM for single-subject fUSI analysis.
second_level –

Second-level GLM for group-level fUSI analysis.

Classes:

Contrast –

Passive container for contrast results.
FirstLevelModel –

First-level GLM estimator for voxel-wise fUSI analysis.
RegressionResults –

Container for regression model results.
SecondLevelModel –

Second-level GLM estimator for group-level fUSI analysis.

Functions:

claron2021_hrf –

Return the Claron et al. (2021) inverse-gamma fUSI HRF preset.
gamma_difference_hrf –

Return an HRF modeled as the difference of two gamma functions.
gamma_hrf –

Return a positive-only gamma HRF on an oversampled time grid.
glover_hrf –

Return the Glover canonical HRF on an oversampled time grid.
inverse_gamma_hrf –

Return an HRF modeled as an inverse gamma distribution.
make_first_level_design_matrix –

Create a first-level design matrix from events and confounds.
make_second_level_design_matrix –

Build a second-level design matrix.
spm_hrf –

Return the SPM canonical HRF on an oversampled time grid.
verhoef2025_hrf –

Return the Verhoef et al. (2025) single-gamma human fUSI HRF preset.

Contrast `dataclass` ¶

Passive container for contrast results.

Holds a contrast's effect/variance estimates and the precomputed test statistics. To construct one from raw (effect, variance, dof), use from_estimate; that is where the t/F statistic, p-value, and z-score are computed.

Supports addition for fixed-effects combination across runs, and scalar multiplication/division for rescaling.

Attributes:

effect ((n_voxels,) numpy.ndarray or (dim, n_voxels) numpy.ndarray) –

Contrast effect estimates.
variance ((n_voxels,) numpy.ndarray) –

Contrast variance estimates.
dim (int) –

Contrast dimension (1 for t, >1 for F).
dof (float) –

Degrees of freedom of the residuals.
stat_type ({'t', 'F'}) –

Type of contrast statistic.
baseline (float) –

Null-hypothesis value used to compute the statistic.
statistic ((n_voxels,) numpy.ndarray) –

Test statistic (t or F).
pvalue ((n_voxels,) numpy.ndarray) –

One-sided p-values (upper tail).
one_minus_pvalue ((n_voxels,) numpy.ndarray) –

1 - pvalue, computed via the CDF for numerical stability in the lower tail.
zscore ((n_voxels,) numpy.ndarray) –

z-scores derived from pvalue/one_minus_pvalue.
tiny (float) –

Small value used to guard against division by zero.
dofmax (float) –

Cap on degrees of freedom used when computing p-values.

Methods:

from_estimate –

Build a Contrast from raw estimates.
from_results –

Build a Contrast from regression results.

from_estimate `classmethod` ¶

from_estimate(
    effect: NDArray[floating],
    variance: NDArray[floating],
    *,
    dim: int | None = None,
    dof: float = 10000000000.0,
    stat_type: Literal["t", "F"] = "t",
    baseline: float = 0.0,
    tiny: float = 1e-50,
    dofmax: float = 10000000000.0,
) -> Contrast

Build a Contrast from raw estimates.

Computes the test statistic, p-values, and z-scores for the supplied baseline and packs everything into a Contrast result object.

Parameters:

effect ¶
((n_voxels,) or (dim, n_voxels) numpy.ndarray) –

Contrast effect estimates.
variance ¶
((n_voxels,) numpy.ndarray) –

Contrast variance estimates.
dim ¶
(int, default: None ) –

Contrast dimension. If not provided, inferred from effect.
dof ¶
(float, default: 1e10 ) –

Degrees of freedom of the residuals.
stat_type ¶
(('t', 'F'), default: "t" ) –

Type of contrast statistic.
baseline ¶
(float, default: 0.0 ) –

Null-hypothesis value tested against. The statistic is computed as (effect - baseline) / sqrt(variance) (t) or sum((effect - baseline)**2) / dim / variance (F). To test a different null, call from_estimate again with a new baseline.
tiny ¶
(float, default: 1e-50 ) –

Small value to guard against division by zero.
dofmax ¶
(float, default: 1e10 ) –

Maximum degrees of freedom used when computing p-values.

Returns:

Contrast –

Contrast result with statistic, pvalue, one_minus_pvalue, and zscore precomputed.

Raises:

ValueError –

If effect is not 1D or 2D, variance is not 1D, or stat_type is not "t" or "F".

from_results `classmethod` ¶

from_results(
    results: RegressionResults,
    contrast_vec: NDArray[floating],
    *,
    stat_type: Literal["t", "F"] | None = None,
    baseline: float = 0.0,
) -> Contrast

Build a Contrast from regression results.

Dispatches to a t or F contrast based on the shape of contrast_vec (or stat_type when provided), then packages the per-voxel effect, variance, and degrees of freedom into a Contrast. For F-contrasts, the whitened-effect representation is used so that the standard sum(effect²) / dim / variance formula in from_estimate recovers the proper quadratic-form F statistic — see compute_f_contrast for the whitening derivation.

Parameters:

results ¶
(RegressionResults) –

Fitted GLM results.
contrast_vec ¶
((n_columns,) or (q, n_columns) numpy.ndarray) –

Numeric contrast vector (1D → t) or matrix (2D → F).
stat_type ¶
(('t', 'F'), default: "t" ) –

Force the contrast type. By default, inferred from the shape of contrast_vec.
baseline ¶
(float, default: 0.0 ) –

Null-hypothesis value tested against.

Returns:

Contrast –

Result with statistic, pvalue, one_minus_pvalue, and zscore precomputed.

FirstLevelModel ¶

Bases: BaseEstimator

First-level GLM estimator for voxel-wise fUSI analysis.

Fits a General Linear Model to fUSI DataArrays and computes statistical contrasts. Supports multiple runs (fixed-effects combination) and autoregressive noise modelling.

This implementation is adapted from nilearn.glm.first_level.FirstLevelModel.

Parameters:

hrf_model ¶
(('glover', 'spm', 'verhoef2025', 'claron2021', 'fir'), default: "glover" ) –

Hemodynamic response function model. A callable matching the [HRFModel][confusius.glm._hrf_models.HRFModel] protocol (a function taking dt and oversampling and returning a 1D array) is invoked to produce a custom HRF kernel. If not provided, skips HRF convolution and uses the raw block regressors, matching make_first_level_design_matrix.
drift_model ¶
(('cosine', 'polynomial'), default: "cosine" ) –

Drift model for low-frequency confounds.
low_cutoff ¶
(float, default: 0.01 ) –

Low cutoff frequency in hertz (used with drift_model="cosine").
drift_order ¶
(int, default: 1 ) –

Polynomial order when drift_model="polynomial".
fir_delays ¶
(list of int, default: None ) –

FIR delays in volumes (required when hrf_model="fir").
noise_model ¶
('ols', default: "ols" ) –

Noise model. "arN" estimates AR(N) coefficients from OLS residuals and refits with a whitened model.
minimize_memory ¶
(bool, default: True ) –

Whether to keep only the statistics needed for contrast computation (per-run RegressionResults are discarded after extracting the contrast-relevant fields).
oversampling ¶
(int, default: 50 ) –

Oversampling factor for HRF convolution.
min_onset ¶
(float, default: -24.0 ) –

Minimum onset time in seconds for event regressors.
uniformity_tolerance ¶
(float, default: 1e-2 ) –

Maximum allowed per-interval relative deviation from the median consecutive interval in the run time coordinate (see [get_representative_step][confusius._utils.get_representative_step]). Increase this value to tolerate slight timestamp jitter.
smoothing_fwhm ¶
(float or dict[str, float], default: None ) –

Full width at half maximum of the Gaussian smoothing kernel, in the physical units of the spatial coordinates, applied to each run before model fitting. A scalar smooths all non-time dimensions (caution if the data has a pose dimension); a dict maps dimension names to per-dimension FWHM. Smoothing is delegated to smooth_volume, which requires each smoothed dimension to have uniform coordinate spacing.

Warning

If the data has any other dimension besides time, a scalar smoothing_fwhm will smooth these dimensions as well.

Attention

The dict may also include the "time" key to smooth along the time axis as well (this is intended behavior).

If not provided, no smoothing is applied.
mask ¶
(DataArray, default: None ) –

Boolean spatial mask selecting voxels to include in model fitting. Must match the spatial dimensions and coordinates of each run.

Attributes:

design_matrices_ (list of pandas.DataFrame) –

Design matrix for each run (set after fit).
results_ (list of RegressionResults) –

Per-run regression results (set after fit).

Examples:

>>> import numpy as np, pandas as pd, xarray as xr
>>> data = xr.DataArray(
...     np.random.randn(200, 2, 3, 4),
...     dims=["time", "z", "y", "x"],
...     coords={"time": np.arange(200) * 0.5},
... )
>>> events = pd.DataFrame(
...     {"trial_type": ["A", "B"] * 5,
...      "onset": np.arange(10) * 10.0,
...      "duration": [2.0] * 10}
... )
>>> model = FirstLevelModel(noise_model="ols")
>>> model.fit(data, events=events)
>>> z_map = model.compute_contrast("A - B")

Methods:

compute_contrast –

Compute a contrast and return a statistical map.
fit –

Fit the GLM to fUSI data.
get_metadata_routing –

Get metadata routing of this object.
get_params –

Get parameters for this estimator.
set_params –

Set the parameters of this estimator.

compute_contrast ¶

compute_contrast(
    contrast_def: str | NDArray[floating],
    stat_type: Literal["t", "F"] | None = None,
    output_type: Literal[
        "zscore",
        "statistic",
        "pvalue",
        "effect",
        "variance",
    ] = "zscore",
    baseline: float = 0.0,
) -> DataArray

Compute a contrast and return a statistical map.

Parameters:

contrast_def ¶
(str or ndarray) –

Contrast definition. A string is parsed as an expression over design-matrix column names (e.g. "A - B"). A 1D array specifies a t-contrast vector; a 2D array specifies an F-contrast matrix.
stat_type ¶
(('t', 'F'), default: "t" ) –

Force the contrast type. By default inferred from the shape of the contrast vector (1D → t, 2D → F).
output_type ¶
(('zscore', 'statistic', 'pvalue', 'effect', 'variance'), default: "zscore" ) –

Which statistical map to return.
baseline ¶
(float, default: 0.0 ) –

Null-hypothesis value tested against. The statistic is (effect - baseline) / sqrt(variance) for t-contrasts and sum((effect - baseline)**2) / dim / variance for F-contrasts.

Returns:

DataArray –

Statistical map with the spatial dimensions of the input data.

Raises:

ValueError –

If the model has not been fitted or the contrast definition is invalid.

fit ¶

fit(
    run_data: DataArray | list[DataArray],
    events: DataFrame | list[DataFrame] | None = None,
    confounds: DataFrame
    | NDArray[floating]
    | list[DataFrame | NDArray[floating] | None]
    | None = None,
    design_matrices: DataFrame
    | list[DataFrame]
    | None = None,
) -> FirstLevelModel

Fit the GLM to fUSI data.

Either events or design_matrices must be provided. When design_matrices is given, events, confounds, and the design-related constructor parameters (hrf_model, etc.) are ignored.

Parameters:

run_data ¶
(xarray.DataArray or list of xarray.DataArray) –

Single-run or multi-run fUSI data. Must have a time dimension; all other dimensions are treated as spatial (e.g. (time, z, y, x) or (time, pose, z, y, x)).
events ¶
(pandas.DataFrame or list of pandas.DataFrame, default: None ) –

Events table(s) with onset, duration, and trial_type columns. Onsets are in the same physical time units as the time coordinate of run_data.
confounds ¶
(pandas.DataFrame, numpy.ndarray, or list, default: None ) –

Confound regressors per run.
design_matrices ¶
(pandas.DataFrame or list of pandas.DataFrame, default: None ) –

Pre-built design matrices. Overrides events / confounds.

Returns:

FirstLevelModel –

Fitted estimator (self).

Raises:

ValueError –

If neither events nor design_matrices is provided, if list lengths are inconsistent across arguments, or if runs have different spatial shapes or mismatched spatial coordinates.

get_metadata_routing ¶

get_metadata_routing()

Get metadata routing of this object.

Please check :ref:User Guide <metadata_routing> on how the routing mechanism works.

Returns:

routing ( MetadataRequest ) –

A :class:~sklearn.utils.metadata_routing.MetadataRequest encapsulating routing information.

get_params ¶

get_params(deep=True)

Get parameters for this estimator.

Parameters:

deep ¶
(bool, default: True ) –

If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns:

params ( dict ) –

Parameter names mapped to their values.

set_params ¶

set_params(**params)

Set the parameters of this estimator.

The method works on simple estimators as well as on nested objects (such as :class:~sklearn.pipeline.Pipeline). The latter have parameters of the form <component>__<parameter> so that it's possible to update each component of a nested object.

Parameters:

**params ¶
(dict, default: {} ) –

Estimator parameters.

Returns:

self ( estimator instance ) –

Estimator instance.

RegressionResults ¶

Container for regression model results.

Stores fitted parameters, residuals, and provides methods for computing contrasts and statistical tests.

This implementation is adapted from [nilearn.glm.regression.RegressionResults][ nilearn.glm.regression.RegressionResults] and [nilearn.glm.model.LikelihoodModelResults][ nilearn.glm.model.LikelihoodModelResults].

Parameters:

theta ¶
((n_regressors, n_voxels) numpy.ndarray) –

Parameter estimates (beta coefficients).
Y ¶
((n_timepoints, n_voxels) numpy.ndarray) –

Original data.
model ¶
(OLSModel or ARModel) –

Fitted model instance.
whitened_Y ¶
((n_timepoints, n_voxels) numpy.ndarray) –

Whitened data.
whitened_residuals ¶
((n_timepoints, n_voxels) numpy.ndarray) –

Whitened residuals.
dispersion ¶
((n_voxels,) numpy.ndarray) –

Residual variance estimate per voxel.
cov ¶
((n_regressors, n_regressors) numpy.ndarray) –

Normalized covariance matrix.

Attributes:

theta ((n_regressors, n_voxels) numpy.ndarray) –

Parameter estimates.
cov ((n_regressors, n_regressors) numpy.ndarray) –

Covariance matrix.
dispersion ((n_voxels,) numpy.ndarray) –

Variance estimates.
df_total (int) –

Total degrees of freedom.
df_model (int) –

Model degrees of freedom.
df_residuals (int) –

Residual degrees of freedom.

Methods:

compute_f_contrast –

Compute F-statistic for a multi-dimensional contrast.
compute_t_contrast –

Compute t-statistic for a contrast.

mse `property` ¶

mse: NDArray[floating]

Mean squared error: SSE / df_residuals.

Returns:

(n_voxels,) numpy.ndarray –

Mean squared error per voxel.

predicted `property` ¶

predicted: NDArray[floating]

Fitted values: X @ theta.

Returns:

(n_timepoints, n_voxels) numpy.ndarray –

Fitted values.

Raises:

RuntimeError –

If the model reference was dropped via minimize_memory=True.

residuals `property` ¶

residuals: NDArray[floating]

Raw residuals: Y - X @ theta.

Returns:

(n_timepoints, n_voxels) numpy.ndarray –

Raw residuals.

Raises:

RuntimeError –

If the heavy fields were stripped via minimize_memory=True.

sse `property` ¶

sse: NDArray[float64]

Sum of squared errors (RSS).

Returns:

(n_voxels,) numpy.ndarray –

Sum of squared whitened residuals per voxel.

Raises:

RuntimeError –

If whitened_residuals was dropped via minimize_memory=True.

compute_f_contrast ¶

compute_f_contrast(
    contrast: NDArray[floating],
) -> dict[str, NDArray[floating]]

Compute F-statistic for a multi-dimensional contrast.

Parameters:

contrast ¶
((q, n_regressors) numpy.ndarray) –

Contrast matrix with q rows (q > 1 for F-test).

Returns:

dict –
Dictionary with keys:
- "effect": Raw contrast effect c·θ, shape (q, n_voxels).
- "whitened_effect": Effect whitened by cholesky(invcov), shape (q, n_voxels). Designed so that sum(whitened_effect**2, axis=0) / q / dispersion == F, which is what Contrast consumes downstream.
- "covariance": Per-voxel contrast covariance with dispersion, shape (n_voxels, q, q).
- "dispersion": Per-voxel residual variance, shape (n_voxels,).
- "F": F-statistic (n_voxels,).
- "df_num": Numerator degrees of freedom (q).
- "df_den": Denominator degrees of freedom.

compute_t_contrast ¶

compute_t_contrast(
    contrast: NDArray[floating],
) -> dict[str, NDArray[floating]]

Compute t-statistic for a contrast.

Parameters:

contrast ¶
((n_regressors,) or (1, n_regressors) numpy.ndarray) –

Contrast vector defining the linear combination of parameters.

Returns:

dict –

Dictionary with keys: - "effect": Contrast effect size (scalar or (n_voxels,)). - "sd": Standard error of contrast (scalar or (n_voxels,)). - "t": t-statistic (scalar or (n_voxels,)). - "df_den": Denominator degrees of freedom.

SecondLevelModel ¶

Bases: BaseEstimator

Second-level GLM estimator for group-level fUSI analysis.

Fits a voxel-wise OLS model to a collection of first-level contrast maps (one per subject or session). This is equivalent to a mass-univariate one-sample t-test by default, but any linear group-level model can be specified via a design matrix.

second_level_input accepts either a list of fitted FirstLevelModel objects or a list of spatial xarray.DataArray maps (e.g. output of compute_contrast with output_type="effect"). When FirstLevelModel objects are passed, first_level_contrast must be provided so the effect map can be extracted automatically from each model.

This implementation is adapted from [nilearn.glm.second_level.SecondLevelModel][ nilearn.glm.second_level.SecondLevelModel].

Attributes:

design_matrix_ (DataFrame) –

Design matrix used for fitting (set after fit).
results_ (RegressionResults) –

Regression results (set after fit).

Examples:

>>> import numpy as np, xarray as xr
>>> maps = [
...     xr.DataArray(np.random.randn(2, 3, 4), dims=["z", "y", "x"])
...     for _ in range(10)
... ]
>>> model = SecondLevelModel()
>>> model.fit(maps)
>>> z_map = model.compute_contrast("intercept")

Methods:

compute_contrast –

Compute a group-level contrast and return a statistical map.
fit –

Fit the group-level GLM.
get_metadata_routing –

Get metadata routing of this object.
get_params –

Get parameters for this estimator.
set_params –

Set the parameters of this estimator.

compute_contrast ¶

compute_contrast(
    second_level_contrast: str
    | NDArray[floating] = "intercept",
    stat_type: Literal["t", "F"] | None = None,
    output_type: Literal[
        "zscore",
        "statistic",
        "pvalue",
        "effect",
        "variance",
    ] = "zscore",
    baseline: float = 0.0,
) -> DataArray

Compute a group-level contrast and return a statistical map.

Parameters:

second_level_contrast ¶
(str or ndarray, default: "intercept" ) –

Contrast definition. A string is parsed as an expression over design-matrix column names (e.g. "group_A - group_B"). A 1D array specifies a t-contrast vector; a 2D array specifies an F-contrast matrix. Defaults to "intercept" for a one-sample t-test.
stat_type ¶
(('t', 'F'), default: "t" ) –

Force the contrast type. By default inferred from the shape of the contrast (1D → t, 2D → F).
output_type ¶
(('zscore', 'statistic', 'pvalue', 'effect', 'variance'), default: "zscore" ) –

Which statistical map to return.
baseline ¶
(float, default: 0.0 ) –

Null-hypothesis value tested against. The statistic is (effect - baseline) / sqrt(variance) for t-contrasts and sum((effect - baseline)**2) / dim / variance for F-contrasts.

Returns:

DataArray –

Statistical map with the spatial dimensions of the input maps.

Raises:

ValueError –

If the model has not been fitted or the contrast definition is invalid.

fit ¶

fit(
    second_level_input: list[DataArray]
    | list[FirstLevelModel],
    first_level_contrast: str
    | NDArray[floating]
    | None = None,
    confounds: DataFrame | None = None,
    design_matrix: DataFrame | None = None,
) -> SecondLevelModel

Fit the group-level GLM.

Accepts either a list of fitted FirstLevelModel objects or a list of spatial contrast maps (xarray.DataArray with no time dimension). All maps must share the same spatial dimensions and shape.

Parameters:

second_level_input ¶
(list of FirstLevelModel or list of xarray.DataArray) –

Input data. When a list of FirstLevelModel is provided, first_level_contrast must be given and effect maps are extracted automatically via compute_contrast. When a list of xarray.DataArray is provided, maps are used directly.
first_level_contrast ¶
(str or ndarray, default: None ) –

Contrast passed to each FirstLevelModel.compute_contrast call when second_level_input is a list of FirstLevelModel. Ignored otherwise.
confounds ¶
((n_subjects, n_confounds) pandas.DataFrame, default: None ) –

Subject-level confound regressors in the same order as second_level_input. Ignored when design_matrix is provided.
design_matrix ¶
((n_subjects, n_regressors) pandas.DataFrame, default: None ) –

Pre-built group-level design matrix. When provided, confounds is ignored.

Returns:

SecondLevelModel –

Fitted estimator (self).

Raises:

ValueError –

If second_level_input is empty, if FirstLevelModel objects are passed without first_level_contrast, if maps have inconsistent spatial shapes or mismatched spatial coordinates, or if design_matrix row count does not match the number of inputs.
TypeError –

If second_level_input is not a list of FirstLevelModel or xarray.DataArray.

get_metadata_routing ¶

get_metadata_routing()

Get metadata routing of this object.

Please check :ref:User Guide <metadata_routing> on how the routing mechanism works.

Returns:

routing ( MetadataRequest ) –

A :class:~sklearn.utils.metadata_routing.MetadataRequest encapsulating routing information.

get_params ¶

get_params(deep=True)

Get parameters for this estimator.

Parameters:

deep ¶
(bool, default: True ) –

If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns:

params ( dict ) –

Parameter names mapped to their values.

set_params ¶

set_params(**params)

Set the parameters of this estimator.

The method works on simple estimators as well as on nested objects (such as :class:~sklearn.pipeline.Pipeline). The latter have parameters of the form <component>__<parameter> so that it's possible to update each component of a nested object.

Parameters:

**params ¶
(dict, default: {} ) –

Estimator parameters.

Returns:

self ( estimator instance ) –

Estimator instance.

claron2021_hrf ¶

claron2021_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    alpha: float = 2.5,
    beta: float = 12.7,
    onset: float = 0.0,
) -> NDArray[floating]

Return the Claron et al. (2021) inverse-gamma fUSI HRF preset.

HRF model proposed for functional ultrasound imaging of the rodent spinal cord in Claron et al. (2021).

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
alpha ¶
(float, default: 2.5 ) –

Shape parameter of the inverse gamma distribution.
beta ¶
(float, default: 12.7 ) –

Scale parameter of the inverse gamma distribution.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

References

Claron, J., Hingot, V., Rivals, I., Rahal, L., Couture, O., Deffieux, T., Tanter, M., & Pezet, S. (2021). Large-scale functional ultrasound imaging of the spinal cord reveals in-depth spatiotemporal responses of spinal nociceptive circuits in both normal and inflammatory states. Pain, 162(4), 1047-1059. https://doi.org/10.1097/j.pain.0000000000002078 ↩

gamma_difference_hrf ¶

gamma_difference_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    onset: float = 0.0,
    delay: float = 6.0,
    undershoot: float = 16.0,
    dispersion: float = 1.0,
    undershoot_dispersion: float = 1.0,
    ratio: float = 1.0 / 6.0,
) -> NDArray[floating]

Return an HRF modeled as the difference of two gamma functions.

The general parameterization underlying glover_hrf and spm_hrf; expose the parameters directly to build custom double-gamma HRFs.

Parameters:

dt ¶
(float) –

Sampling interval in seconds.
oversampling ¶
(int, default: 50 ) –

Temporal oversampling factor relative to the acquisition grid.
time_length ¶
(float, default: 32.0 ) –

Duration of the HRF in seconds.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.
delay ¶
(float, default: 6.0 ) –

Peak delay of the first gamma in seconds.
undershoot ¶
(float, default: 16.0 ) –

Peak delay of the undershoot gamma in seconds.
dispersion ¶
(float, default: 1.0 ) –

Dispersion of the peak gamma.
undershoot_dispersion ¶
(float, default: 1.0 ) –

Dispersion of the undershoot gamma.
ratio ¶
(float, default: 1.0/6.0 ) –

Ratio of undershoot to peak amplitude.

Returns:

(n_timepoints,) numpy.ndarray –

Normalized HRF sampled on an oversampled time grid.

gamma_hrf ¶

gamma_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    peak_delay: float = 5.0,
    dispersion: float = 1.0,
    onset: float = 0.0,
) -> NDArray[floating]

Return a positive-only gamma HRF on an oversampled time grid.

HRF model for functional ultrasound signals that show a single positive lobe without a BOLD-like post-stimulus undershoot. The gamma distribution is parameterized so that its mode occurs peak_delay seconds after onset and its scale is dispersion seconds.

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
peak_delay ¶
(float, default: 5.0 ) –

Time in seconds at which the gamma HRF peaks.
dispersion ¶
(float, default: 1.0 ) –

Scale parameter of the gamma distribution in seconds.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

glover_hrf ¶

glover_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    onset: float = 0.0,
) -> NDArray[floating]

Return the Glover canonical HRF on an oversampled time grid.

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

inverse_gamma_hrf ¶

inverse_gamma_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    alpha: float = 2.5,
    beta: float = 12.7,
    onset: float = 0.0,
) -> NDArray[floating]

Return an HRF modeled as an inverse gamma distribution.

Continuous-time inverse-gamma HRF whose sampled kernel is normalized to unit mass for GLM convolution. This corresponds to the shape beta**alpha / gamma(alpha) * t**(-(alpha + 1)) * exp(-beta / t) for t > 0, omitting any overall amplitude prefactor.

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
alpha ¶
(float, default: 2.5 ) –

Shape parameter of the inverse gamma distribution.
beta ¶
(float, default: 12.7 ) –

Scale parameter of the inverse gamma distribution.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

make_first_level_design_matrix ¶

make_first_level_design_matrix(
    volume_times: NDArray[floating],
    events: DataFrame | None = None,
    hrf_model: HrfModelSpec = None,
    drift_model: DriftModelSpec = "cosine",
    low_cutoff: float = 0.01,
    drift_order: int = 1,
    fir_delays: list[int] | None = None,
    confounds: NDArray[floating] | DataFrame | None = None,
    confound_names: list[str] | None = None,
    oversampling: int = 50,
    min_onset: float = -24.0,
    uniformity_tolerance: float = 0.01,
) -> DataFrame

Create a first-level design matrix from events and confounds.

Parameters:

volume_times ¶
((n_volumes,) numpy.ndarray) –

Acquisition time of each volume in seconds.
events ¶
(DataFrame, default: None ) –

Events table with onset, duration, and trial_type columns.
hrf_model ¶
(('glover', 'spm', 'verhoef2025', 'claron2021', 'fir'), default: "glover" ) –

Hemodynamic response function model. A callable matching the [HRFModel][confusius.glm._hrf_models.HRFModel] protocol (i.e., a function taking dt and oversampling and returning a 1D array) is invoked to produce a custom HRF kernel. If not provided, no convolution is performed and the design matrix will contain boxcar regressors sampled at the acquisition times.
drift_model ¶
(('cosine', 'polynomial'), default: "cosine" ) –

Drift model for low-frequency confounds.
low_cutoff ¶
(float, default: 0.01 ) –

Low cutoff frequency in hertz (used with drift_model="cosine").
drift_order ¶
(int, default: 1 ) –

Polynomial order when drift_model="polynomial".
fir_delays ¶
(list of int, default: None ) –

FIR delays in volumes (required when hrf_model="fir").
confounds ¶
((n_volumes, n_confounds) numpy.ndarray or pandas.DataFrame, default: None ) –

Confound regressors.
confound_names ¶
(list of str, default: None ) –

Names for confound regressors. Inferred from DataFrame columns if not given.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for HRF convolution.
min_onset ¶
(float, default: -24.0 ) –

Minimum onset time in seconds for event regressors.
uniformity_tolerance ¶
(float, default: 1e-2 ) –

Maximum allowed per-interval relative deviation from the median consecutive interval in volume_times (see [get_representative_step][confusius._utils.coordinates.get_representative_step]). Raise a ValueError if any interval exceeds this threshold. Increase this value to tolerate slight timestamp jitter (e.g. from acquisition clocks).

Returns:

(n_volumes, n_columns) pandas.DataFrame –

Design matrix indexed by volume_times. Columns are, in order: one regressor per condition in events (or per FIR delay when hrf_model="fir"), then any confounds, then the drift regressors and the constant column.

make_second_level_design_matrix ¶

make_second_level_design_matrix(
    n_subjects: int, confounds: DataFrame | None = None
) -> DataFrame

Build a second-level design matrix.

Creates a design matrix for group-level analysis. By default, includes only an intercept column (one-sample t-test / group mean). If confounds are provided, they are prepended as additional regressors.

Parameters:

n_subjects ¶
(int) –

Number of subjects (rows in the design matrix).
confounds ¶
((n_subjects, n_confounds) pandas.DataFrame, default: None ) –

Subject-level confound regressors. Columns are added before the intercept.

Returns:

(n_subjects, n_confounds + 1) pandas.DataFrame –

Design matrix.

Raises:

ValueError –

If confounds has a number of rows different from n_subjects.

Examples:

>>> import pandas as pd
>>> dm = make_second_level_design_matrix(5)
>>> list(dm.columns)
['intercept']
>>> dm.shape
(5, 1)

>>> confounds = pd.DataFrame({"age": [25, 30, 35, 40, 45]})
>>> dm = make_second_level_design_matrix(5, confounds=confounds)
>>> list(dm.columns)
['age', 'intercept']

spm_hrf ¶

spm_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    onset: float = 0.0,
) -> NDArray[floating]

Return the SPM canonical HRF on an oversampled time grid.

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

verhoef2025_hrf ¶

verhoef2025_hrf(
    dt: float,
    oversampling: int = 50,
    time_length: float = 32.0,
    peak_delay: float = 5.0,
    dispersion: float = 1.0,
    onset: float = 0.0,
) -> NDArray[floating]

Return the Verhoef et al. (2025) single-gamma human fUSI HRF preset.

HRF model proposed for human 4D functional ultrasound imaging in Verhoef et al. (2025). The preset follows the reported single-gamma HRF, excluding any post-stimulus undershoot.

Parameters:

dt ¶
(float) –

Sampling interval of the original data in seconds.
oversampling ¶
(int, default: 50 ) –

Oversampling factor for the HRF time grid.
time_length ¶
(float, default: 32.0 ) –

Total length of the HRF in seconds.
peak_delay ¶
(float, default: 5.0 ) –

Time in seconds at which the HRF peaks.
dispersion ¶
(float, default: 1.0 ) –

Scale parameter of the gamma distribution in seconds.
onset ¶
(float, default: 0.0 ) –

Onset of the HRF in seconds.

Returns:

(n_samples,) numpy.ndarray –

HRF values on the oversampled time grid, normalized to sum to 1.

References

Verhoef, L., Soloukey, S., Springeling, G., Flikweert, A. J., Lippe, B., de Jong, A. J., Radeljic-Jakic, N., Baas, M., Voorneveld, J., Vincent, A. J. P. E., & Kruizinga, P. (2025). Miniaturized Four-Dimensional Functional Ultrasound for Mapping Human Brain Activity. medRxiv. https://doi.org/10.1101/2025.08.19.25332261 ↩

confusius.glm¶

glm ¶

Contrast dataclass ¶

from_estimate classmethod ¶

effect ¶

variance ¶

dim ¶

dof ¶

stat_type ¶

baseline ¶

tiny ¶

dofmax ¶

from_results classmethod ¶

results ¶

contrast_vec ¶

stat_type ¶

baseline ¶

FirstLevelModel ¶

hrf_model ¶

drift_model ¶

low_cutoff ¶

drift_order ¶

fir_delays ¶

noise_model ¶

minimize_memory ¶

oversampling ¶

min_onset ¶

uniformity_tolerance ¶

smoothing_fwhm ¶

mask ¶

compute_contrast ¶

contrast_def ¶

stat_type ¶

output_type ¶

baseline ¶

fit ¶

run_data ¶

events ¶

confounds ¶

design_matrices ¶

get_metadata_routing ¶

get_params ¶

deep ¶

set_params ¶

**params ¶

RegressionResults ¶

theta ¶

Y ¶

model ¶

whitened_Y ¶

whitened_residuals ¶

dispersion ¶

cov ¶

mse property ¶

predicted property ¶

residuals property ¶

sse property ¶

compute_f_contrast ¶

contrast ¶

compute_t_contrast ¶

contrast ¶

SecondLevelModel ¶

compute_contrast ¶

second_level_contrast ¶

stat_type ¶

output_type ¶

baseline ¶

fit ¶

second_level_input ¶

first_level_contrast ¶

confounds ¶

design_matrix ¶

get_metadata_routing ¶

get_params ¶

deep ¶

set_params ¶

**params ¶

claron2021_hrf ¶

dt ¶

oversampling ¶

`confusius.glm`¶

Contrast `dataclass` ¶

from_estimate `classmethod` ¶

`effect` ¶

`variance` ¶

`dim` ¶

`dof` ¶

`stat_type` ¶

`baseline` ¶

`tiny` ¶

`dofmax` ¶

from_results `classmethod` ¶

`results` ¶

`contrast_vec` ¶

`stat_type` ¶

`baseline` ¶

`hrf_model` ¶

`drift_model` ¶

`low_cutoff` ¶

`drift_order` ¶

`fir_delays` ¶

`noise_model` ¶

`minimize_memory` ¶

`oversampling` ¶

`min_onset` ¶

`uniformity_tolerance` ¶

`smoothing_fwhm` ¶

`mask` ¶

`contrast_def` ¶

`stat_type` ¶

`output_type` ¶

`baseline` ¶

`run_data` ¶

`events` ¶

`confounds` ¶

`design_matrices` ¶

`deep` ¶

`params`** ¶

`theta` ¶

`Y` ¶

`model` ¶

`whitened_Y` ¶

`whitened_residuals` ¶

`dispersion` ¶

`cov` ¶

mse `property` ¶

predicted `property` ¶

residuals `property` ¶

sse `property` ¶

`contrast` ¶

`contrast` ¶

`second_level_contrast` ¶

`stat_type` ¶

`output_type` ¶

`baseline` ¶

`second_level_input` ¶

`first_level_contrast` ¶

`confounds` ¶

`design_matrix` ¶

`deep` ¶

`params`** ¶

`dt` ¶

`oversampling` ¶

`time_length` ¶

`alpha` ¶

`beta` ¶

`onset` ¶

`dt` ¶

`oversampling` ¶

`time_length` ¶

`onset` ¶

`delay` ¶

`undershoot` ¶

`dispersion` ¶

`undershoot_dispersion` ¶

`ratio` ¶

`dt` ¶

`oversampling` ¶

`time_length` ¶

`peak_delay` ¶