Across-Trials Testing with GLHMM toolbox

In this tutorial, we’ll explore how to perform across-trials testing using the GLHMM toolbox, based on the paper The Gaussian-Linear Hidden Markov Model: a Python Package. This test is designed to study differences between trials within one or more experimental sessions. It’s useful for analysing how brain dynamics vary across different stimuli or task conditions.

To focus on the statistical testing, we use synthetic data generated using the Genephys toolbox by Vidaurre (2023). Genephys simulates EEG or MEG data in psychometric experiments, where a subject is repeatedly exposed to stimuli,

For this tutorial:

• \(D\): Gamma values — the decoded state probabilities from an HMM. These were computed from synthetic EEG/MEG-like data across 1421 trials (250 timepoints per trial), grouped into 10 sessions.

• \(R\): categorical labels representing the stimulus shown on each trial (e.g., animate vs inanimate objects).

Note: Due to rendering issues when viewing this notebook through github, internal links, like the table of contents, may not work correctly. To ensure that the notebook renders correctly, you can view it through this link.

Authors: Nick Yao Larsen nylarsen@cfin.au.dk

Table of Contents

1. Preparation

2. Load data

   • Data restructuring

3. Across-trials testing

   • Across-trials testing - Multivariate

   • Across-trials testing - Univariate

1. Preparation

If you don’t have the GLHMM package installed, this notebook will help you install it automatically.

We will also download example data from the Open Science Framework (OSF).

---

# Install packages if needed
try:
    import google.colab
    IN_COLAB = True
except ImportError:
    IN_COLAB = False

if IN_COLAB:
    print("Running in Google Colab. Installing GLHMM...")
    !pip install git+https://github.com/vidaurre/glhmm

# Install osfclient if missing
try:
    import osfclient
except ImportError:
    print("Installing osfclient...")
    !pip install osfclient

# Now import everything we need
import numpy as np
from pathlib import Path
from glhmm import graphics, statistics
np.random.seed(0)  # For reproducibility

2. Load data

To get started, we’ll download the synthetic data needed use it to run the across-sessions test. If they already exist, we will skip downloading.

We’ll use the osfclient package to fetch the files directly from the Open Science Framework (OSF). If you prefer, you can also download them manually from the OSF project page.

# Set up data directory
data_dir = Path.cwd() / "files" / "data_statistical_testing"

if not data_dir.exists():
    print(f"Creating {data_dir}...")
    data_dir.mkdir(parents=True, exist_ok=True)
else:
    print(f"Data directory {data_dir} already exists.")

# Files to download
files = [
    "Gamma_trials.npy",
    "R_trials.npy",
    "idx_trials_session.npy",
]

# Download the files from OSF if they don't exist locally
for fname in files:
    local_path = data_dir / fname
    remote_path = f"data_statistical_testing/{fname}"

    if local_path.exists():
        print(f"✓ {fname} already exists — skipping.")
    else:
        print(f"Downloading {fname}...")
        # as_posix() ensures forward slashes on Windows for shell compatibility
        !osf -p 8qcyj fetch {remote_path} {local_path.as_posix()}

Data directory c:\Users\au323479\Github\glhmm_28_04\docs\notebooks\files\data_statistical_testing already exists.
✓ Gamma_trials.npy already exists — skipping.
✓ R_trials.npy already exists — skipping.
✓ idx_trials_session.npy already exists — skipping.

Load and check the data The folder data_statistical_testing contains:

• ``Gamma_trials``: decoded HMM state probabilities (Gamma values)

  Shape (n_timepoints × n_trials, n_states) : 355,250 × 6 — where 355,250 = 250 timepoints × 1421 subjects

  Each row is a timepoint; each column corresponds to one of 6 HMM states

• ``R_trials``: stimulus condition label for each trial (0 or 1) Shape (n_trials,): 1D array with one value for each of 1421 trials

• ``idx_trials_session``: session boundaries

  Shape (n_sessions, 2): 2D array with 10 rows — each row defines the start and end index for a session

Note: For details on training from scratch, follow the tutorials Standard Gaussian Hidden Markov Model or Gaussian-Linear Hidden Markov Model. In this notebook, we use precomputed Gamma values to focus on the across-sessions statistical test.

Let’s load the files into memory:

# Define the data folder path
PATH_PARENT = Path.cwd()
PATH_DATA = PATH_PARENT / "files" / "data_statistical_testing"

# Load brain and behavioral data
Gamma = np.load(PATH_DATA / "Gamma_trials.npy")
R_trials = np.load(PATH_DATA / "R_trials.npy")
idx_trials_session = np.load(PATH_DATA / "idx_trials_session.npy")

Check data dimensions

# Check data dimensions
print("Gamma shape:", Gamma.shape)
print("R_trials shape:", R_trials.shape)
print("idx_trials_session shape:", idx_trials_session.shape)

Gamma shape: (355250, 6)
R_trials shape: (1421,)
idx_trials_session shape: (10, 2)

Now that we have our decoded Gamma values, we need to reshape the data back into a trial-wise format. This step is essential because we will run statistical tests at each time point, rather than using HMM-aggregated features like fractional occupancy.

To do this, we’ll use the function reconstruct_concatenated_to_3D. It converts a concatenated 2D array into a structured 3D array — restoring the original format of (timepoints × trials × features).

Convert Gamma to trial-wise format

# Define the dimensions
n_timepoints = 250
n_trials = len(R_trials)
n_features = Gamma.shape[1]

# Reshape Gamma into a 3D array: timepoints × trials × features
Gamma_epoch = statistics.reconstruct_concatenated_to_3D(
    Gamma,
    n_timepoints=n_timepoints,
    n_entities=n_trials,
    n_features=n_features
)

# Check the new shape
print("Gamma_epoch shape:", Gamma_epoch.shape)  # Expected: (250, 1500, 6)

Gamma_epoch shape: (250, 1421, 6)

Sanity check: does the reshaping work?

To confirm the data was reshaped correctly, let’s compare a slice of the new 3D matrix with the original 2D version.

We’ll compare:

Gamma_epoch[:, 0, :]: all timepoints from the first trial in the reshaped 3D matrix

Gamma_sessions[0:250, :]: the first 250 timepoints in the original 2D matrix

np.array_equal(Gamma_epoch[:, 0, :], Gamma[0:250, :])

True

If this returns True, the reshaping was successful — every timepoint from the first trial is correctly aligned in the new format.

3. Across-trials testing

Now that our data is organized into trials and sessions, we’re ready to test whether changes in brain activity (Gamma) are related to behavioral variation across trials using the test_across_trials function.

Before we dive into the test, let’s quickly explain what’s going on.

In this type of test, we want to find out if the brain behaves differently depending on the type of stimulus presented — and whether these differences are consistent across trials within a session. For example, does the brain respond differently when the subject sees an animate object versus an inanimate object?

To answer that, we use permutation testing. This involves shuffling the trial labels within each session and recomputing the statistical test for each shuffled version of the data. By repeating this process many times, we build a null distribution of results. If the result from the original (unshuffled) data stands out from this distribution, it suggests that the brain activity is systematically influenced by the type of stimulus.

Figure 5C in Vidaurre et al. 2023: A 9 x 4 matrix representing permutation testing across sessions. Each number corresponds to a trial within a session and permutations are performed between sessions, with each session containing the same number of trials.

Across-trials testing - Multivariate

We’ll begin by running a multivariate test. This analysis checks whether patterns of brain state activity can explain differences in stimulus condition (e.g., animate vs inanimate) across trials.

If brain responses do not vary meaningfully with trial type, we expect little to no systematic pattern. However, if specific states consistently differ between stimulus categories, this analysis can reveal those relationships.

Inputs:

• D_data: the Gamma matrix reshaped into trials (Gamma_epoch)

• R_data: stimulus labels for each trial (R_trials, values 0 or 1)

• indices_blocks: trial indices for each session (idx_trials_session)

Settings:

• method = "multivariate": perform multivariate regression

• Nnull_samples = 10_000: number of permutations to build the null distribution

# Set the parameters for across trials testing
method = "multivariate"
Nnull_samples = 10_000 # Number of permutations (default = 0)

# Perform across-subject testing
result_multivariate_session  =statistics.test_across_trials(D_data=Gamma_epoch,
                                                            R_data=R_trials,
                                                            indices_blocks=idx_trials_session,
                                                            method=method,
                                                            Nnull_samples=Nnull_samples)

performing permutation testing per timepoint
Total possible permutations: 1.29e+2459
Running number of permutations: 10000

100%|██████████| 250/250 [17:58<00:00,  4.31s/it]

Understanding the output

The result is stored in a dictionary called result_multivariate. Here’s what it contains:

• pval: array of p-values with shape (1, q). Each value corresponds to a behavioral variable. See the GLHMM paper for details.

• base_statistics: test statistics for the unpermuted data. For multivariate tests, this is the F-statistic

• null_stat_distribution: array of test statistics computed from each permutation

• statistical_measures: dictionary showing which statistic was used (e.g. F-stat)

• test_type: Indicates the type of permutation test performed. In this case, it is across_trials.

• method: test method used, here it is "multivariate"

• max_correction: whether MaxT correction was applied during permutation

• Nnull_samples: number of permutations used to generate the null

• test_summary: structured summary of results including F-test and model coefficients

Display the test summary

We can print a clean summary of the result using the helper function below.

statistics.display_test_summary(result_multivariate_session)

Model Summary (timepoint 0):
    Outcome  F-stat  df1  df2  p-value (F-stat)
Regressor 1  0.3603    6 1415            0.8279

Coefficients Table (timepoint 0):
Predictor     Outcome    T-stat  p-value      LLCI     ULCI
  State 1 Regressor 1  0.378125   0.7050 -1.922703 1.873748
  State 2 Regressor 1 -0.571448   0.5629 -1.854683 2.012392
  State 3 Regressor 1  0.380939   0.7068 -1.862957 2.085601
  State 4 Regressor 1  0.734707   0.4692 -2.148568 1.750658
  State 5 Regressor 1  0.083531   0.9339 -1.798056 2.101710
  State 6 Regressor 1 -1.100686   0.2724 -1.998593 1.932921

Visualising the results

We can now plot the p-values over time using the plot_p_values_over_time function from the graphics module.

Lower p-values (below α = 0.05) are shown in warm colours, while higher p-values are shown in cool colours.

graphics.plot_p_values_over_time(
    result_multivariate_session["pval"],
    title_text="Multivariate - Uncorrected",
    xlabel="Time points",
    num_x_ticks=11
)

Multiple testing correction

To reduce the chance of false positives, we apply Benjamini-Hochberg correction to control the false discovery rate:

Visualisation of results

Now that we have performed our test, we can then visualise the p-value array.

We will use the function plot_heatmap from the graphics module.

Note: Elements marked in red indicate a p-value below 0.05, signifying statistical significance.

Multiple Comparison

To take into account for type 1 error, we can apply the Benjamini/Hochberg correction.

pval_corrected, rejected_corrected =statistics.pval_correction(result_multivariate_session,
                                                               method='fdr_bh')
# Plot p-values
graphics.plot_p_values_over_time(pval_corrected,
                                 title_text ="Multivariate - Benjamini/Hochberg",
                                 xlabel="HMM States",
                                 num_x_ticks=11)

Cluster-based correction

To control for family-wise error rate more strictly, we can apply cluster-based permutation testing. This is useful when you expect effects to be temporally contiguous.

Make sure test_statistics_option=True was set earlier, as this is required to access the full permutation distribution:

pval_cluster = statistics.pval_cluster_based_correction(result_multivariate_session)

graphics.plot_p_values_over_time(
    pval_cluster,
    title_text="Multivariate - Cluster correction",
    xlabel="Time points",
    num_x_ticks=11
)

Visualise average probabilities and differences

To help interpret the statistical results from result_multivariate_trials[“pval”], we can compare the average state time courses between the two conditions. This can be done using the function plot_condition_difference, which shows the average probability of each state over time for each condition (0 or 1), along with their difference. We can then highlight time points where the difference is statistically significant based on our multiple testing correction.

# Detect the intervals of when there is a significant difference, will be highlighed
alpha = 0.05
intervals =statistics.detect_significant_intervals(pval_cluster, alpha)
# Plot the average probability for each state over time for two conditions and their difference.
graphics.plot_condition_difference(Gamma_epoch,
                                   R_trials,
                                   figsize=(12,3),
                                   vertical_lines=intervals,
                                   highlight_boxes=True,
                                   title="Average Probability and Differences (Cluster corrected)",)

Across-trials testing - Univariate

Now let’s run the univariate test, which looks at each brain state individually. This helps us identify whether a specific state shows a relationship with the stimulus condition (e.g., animate vs inanimate).

This test runs separately for each state and each timepoint, giving a more detailed view of where and when differences occur.

Inputs:

• D_data: the Gamma matrix reshaped into trials (Gamma_epoch)

• R_data: stimulus labels (R_trials)

• indices_blocks: session boundaries (idx_trials_session)

Settings:

• method = "univariate": run one test per state

• Nnull_samples = 10_000: number of permutations

# Set the parameters for between-subject testing
method = "univariate"
Nnull_samples = 10_000 # Number of permutations (default = 1000)
# Perform across-subject testing
result_univariate =statistics.test_across_trials(D_data=Gamma_epoch,
                                                R_data=R_trials,
                                                indices_blocks=idx_trials_session,
                                                method=method,
                                                Nnull_samples=Nnull_samples)

performing permutation testing per timepoint
Total possible permutations: 1.29e+2459
Running number of permutations: 10000

100%|██████████| 250/250 [09:02<00:00,  2.17s/it]

Visualising univariate results

We can now plot the p-values as a matrix where each row is a state and each column is a timepoint.

This is done by using the plot_p_value_matrix function from the graphics module. > Lower p-values (below α = 0.05) are shown in warm colours, while higher p-values are shown in cool colours.

# Plot p-values
# P-values between reaction time and each state as function of time
graphics.plot_p_value_matrix(result_univariate["pval"].T,
                             title_text ="Univariate - Uncorrected",
                             figsize=(8, 5),
                             xlabel="Time points",
                             ylabel="HMM States",
                             annot=False,
                             num_x_ticks=11)

Multiple Comparison

To be sure there is no type 1 error, we can apply the Benjamini/Hochberg to control the False Discovery Rate

pval_corrected, rejected_corrected =statistics.pval_correction(result_univariate,
                                                               method='fdr_bh')
# Plot p-values
graphics.plot_p_value_matrix(pval_corrected.T,
                             title_text ="Univaraite - Benjamini/Hochberg",
                             figsize=(8, 5),
                             xlabel="Time points",
                             ylabel="HMM States",
                             annot=False,
                             num_x_ticks=11)

Cluster-based correction for univariate test

We can also apply cluster-based correction here to account for contiguous significant effects across time:

pval_cluster =statistics.pval_cluster_based_correction(result_univariate, individual_feature=True,)
# Plot p-values
graphics.plot_p_value_matrix(pval_cluster.T,
                             title_text ="Univariate - Cluster corrected",
                             figsize=(8, 5),
                             xlabel="Time points",
                             ylabel="HMM States",
                             annot=False,
                             num_x_ticks=11)

Family-Wise Error Rate (FWER) correction

Applying FWER correction

# Correct for FWER
pval_FWER = statistics.pval_FWER_correction(result_univariate)

# Correct p-values using FWER
graphics.plot_p_value_matrix(pval_FWER.T,
                             title_text ="Univaraite - FWER",figsize=(8, 5),
                             xlabel="Time points",
                             ylabel="HMM States",
                             annot=False,
                             num_x_ticks=11)