Across-State-Visits Testing with GLHMM toolbox

This tutorial explains how to use the GLHMM toolbox, part of the paper The Gaussian-Linear Hidden Markov Model: a Python Package, to explore whether brain states identified by a Hidden Markov Model (HMM) are related to a continuous signal, like a behavioral measurement.

For simplicity, we use synthetic data instead of real data. The brain data (\(D\)) comes from the Viterbi path, which represents the sequence of brain states identified by the HMM. These states capture patterns in brain activity over time. The behavioral data (\(R\)) is a continuous signal recorded at the same time as the brain data.

The purpose of the test is to determine if the brain states decoded by the HMM are associated with patterns in the behavioral signal. While data preparation requires some explanation, running the test_across_state_visits function is straightforward. Simply provide the input variables (\(D\) and \(R\)) and specify the method for the analysis.

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. Load and prepare data

2. Permutation testing tutorial - Across-Visits

   • Multivariate

   • Univariate

   • One-state-vs-the-rest (OSR)

   • One-state-vs-another-state (OSA)

Install necessary packages

If you dont have the GLHMM-package installed, then run the following command in your terminal:

pip install git+https://github.com/vidaurre/glhmm

Import libraries

Let’s start by importing the required libraries and modules.

# Importing libraries
import numpy as np
from pathlib import Path
from glhmm import glhmm, graphics, statistics

1.Load and prepare data

Example data is available in the folder data_statistical_testing. The file vpath.npy contains the Viterbi path, and sig_data.npy is synthetic timeseries.

The Viterbi path, derived from synthetic brain data, was decoded using a Hidden Markov Model (HMM). Detailed training steps for the HMM model can be found in this tutorial.

To generate the Viterbi path yourself, use the following code after training the model:

vpath = hmm.decode(X=None, Y=D_data, indices=idx, viterbi=True)

For this example, the Viterbi path is precomputed to focus on statistical analysis between the brain data (\(D\)) and synthetic timeseries (\(R\)).

We will now load the data.

# Get the current directory
PATH_PARENT = Path.cwd()
# Path of the location of the data
PATH_DATA = PATH_PARENT / "data_statistical_testing"

# Load vpath (D)
vpath = np.load(PATH_DATA/"vpath.npy")

# Load sig_data (R)
sig_data = np.load(PATH_DATA/"sig_data.npy")

print(f"Data dimension of vpath: {vpath.shape}")
print(f"Data dimension of sig_data: {sig_data.shape}")

Data dimension of vpath: (5000, 5)
Data dimension of sig_data: (5000, 1)

Look at data

Now we can look at the data structure.

• vpath: 2D array of shape (n_timepoints, n_states)

• sig_data: 1D array of shape (n_timepoints,)

vpath has 5000 rows and 5 columns. Each row represents a moment in time, and the columns show which brain state the HMM identified at that moment.

sig_data has 5000 rows and 1 column. Each row contains the value of a continuous signal measured at the same time as the brain data.

Visualize Viterbi path

Now, let’s visualize the distinct states in the Viterbi path from our trained Hidden Markov Model.

The plot provides a clear depiction of each time point assigned to a specific HMM state, with each state represented by a distinct color. This visualization allows us to easily discern the temporal distribution and transitions between different states.

graphics.plot_vpath(vpath, figsize=(6,4), ylabel="States")

Plot Viterbi path and discrete states

This plot shows the brain states identified by the HMM over time, with each state represented by a different color. It makes it easy to see which state the model assigned to each moment.

At the top of the plot, there’s an extra layer that shows the state assignments as numbers on the y-axis. For example, if you see the number 5 on the y-axis in blue, it means the data at that moment belongs to state 5. This format is similar to Figure 3D in our paper.

# Create a 1D array from the one-hot encoded Viterbi path data
vpath_1D=statistics.generate_vpath_1D(vpath)
# Convert the signal to be in a range from 0 to 1,
sig_state = vpath_1D/np.max(vpath_1D)
# Plot discrete states with the Viterbi path
graphics.plot_vpath(vpath,sig_state, figsize=(6,4), ylabel="States", yticks=True)

Plot Viterbi path and signal

This plot shows the Viterbi path (D) together with the signal (R). Each time point is color-coded based on the brain state assigned by the HMM, making it easy to compare the signal and the state changes.

To create this plot, the signal data is adjusted to fit on the same scale as the Viterbi path (from 0 to 1). This adjustment is just for visualization, helping us overlay the two and see how well the model reflects the patterns in the signal.

# Normalize the sig_data to the range [0, 1]
min_value = np.min(sig_data)
max_value = np.max(sig_data)
normalized_sig_data = ((sig_data - min_value) / (max_value - min_value))

# Plot vpath and sig_data
graphics.plot_vpath(vpath,normalized_sig_data, figsize=(6,4))

Let’s delve into the values assigned to each state. To achieve this, we calculate the average of the values associated with each state.

This brief analysis provides insights into the typical or central values characterizing each state. By examining these averages, we gain a clearer understanding of the distinctive characteristics and patterns represented by the different states in our HMM.

vpath_1D=statistics.generate_vpath_1D(vpath)
val_state =[np.mean(sig_data[vpath_1D == i+1]) for i in range(vpath.shape[1]) ]
val_state

[np.float64(4.912042797300882),
 np.float64(-0.05858422375450776),
 np.float64(-3.795696791047377),
 np.float64(2.352225972831181),
 np.float64(-2.5030418147224505)]

The above figure shows an relation between the simulated measurements (sig_data) and the Viterbi path (vpath). Notably, a closer examination of the variable val_state shows distinct values assigned to each specific state:

• State 1 (green) aligns with highest value

• State 2 (yellow) corresponds to values in the middle, close to 0.

• State 3 (purple) corresponds to the lowest negative values.

• States 4 (red) and State 5 (blue) are values that fall between the extremes values. This gives us a perspective on the behavioral dynamics across different states on our signal.

2. Permutation testing tutorial - Across-Visits

As we move on to the next part of this tutorial, let’s dive into how we can use the test_across_state_visits function. This function helps us to find connections between HMM state time course (D) and behavioral variables or individual traits (R) using permutation testing.

Permutation testing

Permutation testing does not assume any particular data distribution and the procedure shuffles the data around to create a null distribution. This null distribution comes in handy for testing our hypotheses without making any assumptions about the data. This null distribution becomes our benchmark to test the question: is there any real difference or relationship between the variables we’re interested in?

Across visits - Multivariate

The multivariate analysis aims to determine how the Viterbi path (D_data), characterized by different states, contributes to explaining the observed variability in the simulated signal (R_data). By quantifying explained variance, this analysis evaluates whether state transitions significantly influence changes in signal values over time.

A permutation test for explained variance assesses the statistical significance of this relationship. A significant result indicates that the Viterbi path meaningfully explains the signal’s variability, while a non-significant result suggests the observed relationship is likely due to random chance.

To run the test, use the across_state_visits function by providing the inputs:

Inputs:

• D_data: The Viterbi path.

• R_data: The simulated signal.

Settings:

• method = "multivariate": Specifies that the test should perform multivariate analysis.

• Nnull_samples: Number of samples (optional, default is 1000).

For details on additional settings, refer to the function documentation.

# Set the parameters for across_visits testing
method = "multivariate"
Nnull_samples = 10_000 # Number of permutations (default = 1000)
result_multivariate  =statistics.test_across_state_visits(vpath,
                                                          sig_data,
                                                          method=method,
                                                          Nnull_samples=Nnull_samples)

100%|██████████| 10001/10001 [00:11<00:00, 907.75it/s]

What we can see here is that result_multivariate is a dictionary containing the outcomes of a statistical analysis conducted using the specified method and test_type.

Let us break it down:

• pval: This array holds the p-values resulting from the permutation test. Each value corresponds to a behavioral variable and will have shape of 1 by q, see GLHMM paper.

• base_statistics: Stores the base statistics of the tests. In this case it is the explained variance \(R^2\)

• test_statistic: Will by default always return a list of the base (unpermuted) statistics when test_statistic_option=False. This list can store the test statistics associated with the permutation test. It provides information about the permutation distribution that is used to calculate the p-values. The output will exported if we set test_statistic_option=True

• statistical_measures: A dictionary that marks the units used as test_statistics

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

• method: Specifies the analytical method employed, which is 'multivariate', which means that the analysis is carried out using regression-based permutation testing.

• max_correction: Boolean value that indicates whether Max correction has been applied when performing permutation testing

• Nnull_samples: Is the number of samples that has been performed.

• test_summary: A dictionary summarizing the test results based on the applied method.

Visualization of results

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

We will import the function plot_p_values_bar from module graphics.py

# Plot p-values
graphics.plot_p_values_bar(result_multivariate["pval"],
                           title_text ="Multivariate - uncorrected",
                           figsize=(3, 3),
                           alpha=0.05,
                           xticklabels=["Signal"])

Conclusion - Multivariate

The multivariate test between the Viterbi path and the signal resulted in a statistically significant p-value. This indicates a meaningful relationship between the brain states identified by the HMM and the signal, suggesting that the states capture important dynamics related to the observed data.

Across visits - Univariate test

The univariate analysis evaluates pairwise relationships between discrete states in the Viterbi path (D_data) and the variability in the signal (R_data). Each state is tested individually against the signal, resulting in one p-value per state. For example, with 5 states and 1 signal, the test will produce 5 p-values.

The goal is to assess whether specific states in the Viterbi path significantly contribute to changes in the signal values over time. A permutation test determines the statistical significance of these relationships:

To run the test, we set the following parameters:

Inputs:

• D_data: The Viterbi path.

• R_data: The simulated signal.

Settings:

• method = "univariate": Specifies that the test should perform multivariate analysis.

• Nnull_samples: Number of samples (optional, default is 1000).

The function test_across_state_visits is used as follows:

# Set the parameters for across_visits testing
method = "univariate"
Nnull_samples = 10_000
result_univariate  =statistics.test_across_state_visits(vpath,
                                                        sig_data,
                                                        method=method,
                                                        Nnull_samples=Nnull_samples)

100%|██████████| 10001/10001 [00:12<00:00, 824.20it/s]

Now that we have the permutation distribution [‘test_statistic] as we can see, it is because we set test_statistic_option=True

Visualization of results

Now that we have performed our test, we can then visualize the p-value array. We will import the function plot_heatmap, plot_scatter_with_labels and plot_histograms from module helperfunctions.py

# Plot p-values
graphics.plot_p_value_matrix(result_univariate["pval"].T,
                             title_text ="Univariate - Uncorrected",
                             figsize=(7, 2.5),
                             xlabel="HMM States",
                             alpha=0.05,
                             x_tick_min=1)

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 ="Univariate - Benjamini/Hochberg",
                             figsize=(7, 2.5),
                             xlabel="HMM States",
                             x_tick_min=1)

Instead of using a heatmap, we can also visualize the results with a bar plot

# Set the threshold of alpha to be 0.05
alpha = 0.05
variables = [f"State {i+1}" for i in range(len(pval_corrected))] # construct the variable names
graphics.plot_p_values_bar(pval_corrected,
                           alpha = alpha,
                           xticklabels=variables,
                           title_text="Univariate - Benjamini/Hochberg")

Plot permutation distribution

Presented here are the test statistics(result["null_stat_distribution"]) of our permutation distributions for different states.

The red line shows the observed statistic, while the datapoints of the histogram represent the permutation distribution

# Plot test statistics for pvals
significant_timestamp_position = np.where(pval_corrected < alpha)
for i in significant_timestamp_position[0]:
    graphics.plot_permutation_distribution(result_univariate["null_stat_distribution"][:,i],
                                           title_text=f"Permutation distribution of State Nr.{i+1} ")

Conclusion - Univariate

The results show that states 1, 3, and 4 have a significant relationship with the signal, while states 2 and 5 do not. The signal’s highest peak aligns with state 1, the lowest with state 3, and the second-highest with state 4, supporting these findings.

This suggests that states with extreme signal values have a stronger influence on the signal’s changes, while states 2 and 5 may have a weaker or less clear impact.

Across visits - One-state-vs-the-rest (OSR)

The OSR test assesses whether a specific state in the Viterbi path (D_data) significantly differs from the combined influence of all other states on the signal (R_data). This test generates five p-values, one for each state.

By default, the test checks if a state is larger than the rest using state_com='larger'. To test if a state is smaller, set state_com='smaller'.

Inputs:

• D_data: The Viterbi path.

• R_data: The simulated signal.

Settings:

• method = "ors": Specifies that the test should perform multivariate analysis.

• Nnull_samples: Number of permutations (optional, default is 1000).

# Set the parameters for across_visits testing
method = "osr"
Nnull_samples = 10_000

result_one_vs_rest  =statistics.test_across_state_visits(vpath,
                                                         sig_data,
                                                         method=method,
                                                         Nnull_samples=Nnull_samples)

100%|██████████| 10001/10001 [00:13<00:00, 722.90it/s]

Visualization of results

Now that we have performed our test, we can then visualize the p-value array. We will plot the p-values using the function plot_p_value_matrix

# Set the threshold of alpha to be 0.05
alpha = 0.05
variables = [f"State {i+1}" for i in range(len(result_one_vs_rest["pval"]))] # construct the variable names
title_text = "One vs rest (Larger) - Uncorrected"
graphics.plot_p_values_bar(result_one_vs_rest["pval"],
                           alpha = alpha,
                           xticklabels=variables,
                           title_text=title_text)

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, _ =statistics.pval_correction(result_one_vs_rest,
                                              method='fdr_bh')
# Plot p-values
graphics.plot_p_values_bar(pval_corrected,alpha = alpha,
                           xticklabels=variables,
                           title_text="One vs rest (Larger) - Benjamini/Hochberg")

val_state

[np.float64(4.912042797300882),
 np.float64(-0.05858422375450776),
 np.float64(-3.795696791047377),
 np.float64(2.352225972831181),
 np.float64(-2.5030418147224505)]

After applying the Benjamini/Hochberg method to control the False Discovery Rate (FDR), the adjusted p-values for the one vs. rest test are:

• State 1: 0.005

• State 2: 0.724

• State 3: 1

• State 4: 0.048

• State 5: 1

States 1 and 4, which have the lowest adjusted p-values, correspond to the highest values in sig_data. A significant p-value in the one vs. rest test indicates that the mean signal for a specific state is significantly different from the combined influence of all other states.

Across visits - One-state-vs-another-state (OSA)

The OSA test compares the mean signal differences between pairs of states (e.g., state 1 vs. state 2, state 1 vs. state 3, etc.). This test evaluates whether the observed mean difference between specific state pairs is statistically significant.

Inputs:

• D_data: The Viterbi path.

• R_data: The simulated signal.

Settings:

• method = "osa": Specifies that the test should perform multivariate analysis.

• Nnull_samples: Number of permutations (optional, default is 1000).

# Set the parameters for across_visits testing
method = "osa"
Nnull_samples = 10_000
test_statistics_option=True

result_state_pairs  =statistics.test_across_state_visits(vpath,
                                                         sig_data,
                                                         method=method,
                                                         Nnull_samples=Nnull_samples,)

Pairwise comparisons: 100%|██████████| 10/10 [01:16<00:00,  7.64s/it]

Visualization of results

Now that we have performed our test, we can then visualize the p-value array. We will plot the p-values using the function plot_heatmap. Notably, in this instance, we designate the values on the diagonal as NaN (Not a Number) since these values are expected to be zeros and can be safely ignored in the visualization.

# Plot p-values# Plot p-values
graphics.plot_p_value_matrix(result_state_pairs["pval"],
                             title_text ="State pairs - Uncorrected",
                             figsize=(7, 3),
                             xlabel="State X",
                             ylabel="State Y",
                             alpha=0.05,
                             none_diagonal=True,
                             annot=True,
                             x_tick_min=1)

Interpret the p-value matrix

The matrix has a shape of (number of states by number of states), and the values below the diagonal (lower-left) represent comparisons where the lower state number is less than the higher state number. Values above the diagonal (upper-right) represent comparisons where the lower state number is greater than the higher state number.

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_state_pairs, method='fdr_bh')
# Plot p-values# Plot p-values
graphics.plot_p_value_matrix(pval_corrected,
                             title_text ="State pairs - Benjamini/Hochberg" ,
                             figsize=(7, 3),
                             xlabel="State X",
                             ylabel="State Y",
                             alpha=0.05,
                             none_diagonal=True,
                             annot=True,
                             x_tick_min=1)

# Look at the val_sates values
val_state

[np.float64(4.912042797300882),
 np.float64(-0.05858422375450776),
 np.float64(-3.795696791047377),
 np.float64(2.352225972831181),
 np.float64(-2.5030418147224505)]

The results show clear differences in the signal between certain pairs of brain states:

• State 1 (the highest signal value) is significantly different from states 2, 3, and 5.

• State 3 (the lowest signal value) is significantly different from states 1 and 4.

• States 4 and 5, which are in the middle, also show a clear difference from each other and from state 1.

These differences suggest that specific pairs of brain states are linked to noticeable changes in the signal. This could help us understand how the signal shifts during transitions between states, such as when studying responses to pain.