Two sample permutation test of respiratory phase angles at the group level#

This notebook runs permutation tests to compare two samples of respiratory phase angles. This is useful when you want to compare the phase angles of two different conditions or groups.

[2]:

import pickle as pkl
import numpy as np
from pathlib import Path

from pyriodic.preproc import RawSignal
from pyriodic.permutation import permutation_test_between_samples
from pyriodic.phase_events import create_phase_events

/Users/au661930/Library/CloudStorage/OneDrive-Aarhusuniversitet/Dokumenter/projects/_BehaviouralBreathing/code/AnalysisBreathingBehaviour/BreathingBehaviourVenv/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

[3]:

# preprocess the data and extract phase angles for each subject for correct and incorrect trials
data_path = Path("../../data/raw/intermediate")
circs = {}
for file_path in data_path.iterdir():
    if file_path.suffix != ".pkl":
        continue
    subj_id = file_path.name.split("_")[1]
    with open(file_path, 'rb') as f:
            data = pkl.load(f)


    behav_data = data["behav_data"]
    resp_ts, sfreq, event_samples, event_ids, labels, correct = data["raw"], data["sfreq"], behav_data["event_samples"], behav_data["event_ids"], behav_data["event_type"], behav_data["correct"]

    start_sample = event_samples[0] - 1000  # start 1000 samples before the first event

    # subtract start sample from events (the first column)
    event_samples -= start_sample


    # initialise RawSignal object
    raw = RawSignal(resp_ts, fs=sfreq)
    raw.filter_bandpass(low=0.1, high=1.0)
    raw.smoothing(window_size=500)
    raw.zscore()

    PA, peaks, troughs = raw.phase_twopoint(prominence=0.1)
    correct_event_samples = [
        samp for samp, label, corr in zip(event_samples, labels, correct) if "target" in label and corr == 1.0
    ]

    correct_targets = create_phase_events(
        PA, events=correct_event_samples
    )

    incorrect_event_samples = [
        samp for samp, label, corr in zip(event_samples, labels, correct) if "target" in label and corr == 0.0
    ]
    incorrect_targets = create_phase_events(
        PA, events=incorrect_event_samples
    )

    circs[subj_id] = {
        "correct": correct_targets,
        "incorrect": incorrect_targets
    }

Rejected 0 out of 273 events (0.0%)
Rejected 0 out of 54 events (0.0%)
Rejected 0 out of 305 events (0.0%)
Rejected 0 out of 73 events (0.0%)
Rejected 0 out of 320 events (0.0%)
Rejected 0 out of 70 events (0.0%)
Rejected 0 out of 297 events (0.0%)
Rejected 0 out of 67 events (0.0%)
Rejected 0 out of 303 events (0.0%)
Rejected 0 out of 71 events (0.0%)
Rejected 0 out of 226 events (0.0%)
Rejected 0 out of 77 events (0.0%)
Rejected 0 out of 242 events (0.0%)
Rejected 0 out of 87 events (0.0%)
Rejected 0 out of 267 events (0.0%)
Rejected 0 out of 88 events (0.0%)
Rejected 0 out of 261 events (0.0%)
Rejected 0 out of 114 events (0.0%)
Rejected 0 out of 212 events (0.0%)
Rejected 0 out of 38 events (0.0%)
Rejected 0 out of 243 events (0.0%)
Rejected 0 out of 79 events (0.0%)
Rejected 0 out of 312 events (0.0%)
Rejected 0 out of 56 events (0.0%)
Rejected 0 out of 293 events (0.0%)
Rejected 0 out of 71 events (0.0%)
Rejected 0 out of 253 events (0.0%)
Rejected 0 out of 78 events (0.0%)
Rejected 0 out of 314 events (0.0%)
Rejected 0 out of 70 events (0.0%)
Rejected 0 out of 316 events (0.0%)
Rejected 0 out of 44 events (0.0%)
Rejected 0 out of 332 events (0.0%)
Rejected 0 out of 66 events (0.0%)
Rejected 0 out of 282 events (0.0%)
Rejected 0 out of 74 events (0.0%)
Rejected 0 out of 312 events (0.0%)
Rejected 0 out of 70 events (0.0%)
Rejected 0 out of 246 events (0.0%)
Rejected 0 out of 113 events (0.0%)
Rejected 0 out of 278 events (0.0%)
Rejected 0 out of 71 events (0.0%)
Rejected 0 out of 294 events (0.0%)
Rejected 0 out of 95 events (0.0%)
Rejected 0 out of 308 events (0.0%)
Rejected 0 out of 62 events (0.0%)
Rejected 0 out of 327 events (0.0%)
Rejected 0 out of 66 events (0.0%)
Rejected 0 out of 330 events (0.0%)
Rejected 0 out of 58 events (0.0%)
Rejected 0 out of 257 events (0.0%)
Rejected 0 out of 100 events (0.0%)

[5]:

n_permutations = 1000

group_observed_stats = []
group_null_stats = []
group_pvals = []


for subj_id in circs:

    correct_targets_tmp = circs[subj_id]["correct"]
    incorrect_targets_tmp = circs[subj_id]["incorrect"]

    # Run subject-level permutation test
    obs, pval, null_distribution = permutation_test_between_samples(
        incorrect_targets_tmp.data,
        correct_targets_tmp.data,
        n_permutations=n_permutations,
        verbose=False,
        return_null_distribution=True,
    )

    group_observed_stats.append(obs)
    group_null_stats.append(null_distribution)
    group_pvals.append(pval)


# Group-level observed statistic: average of subject-level observed stats
# On average, do subjects show evidence of difference in phase for irregular and regular targets?
group_obs = np.mean(group_observed_stats)

# Group-level null: average of null stats across subjects, per null_stat
# Reference distribution of what kind of group-level effect you’d expect if there was no difference between irregular and regular targets
group_null = [
    np.mean([nulls[i] for nulls in group_null_stats])
    for i in range(n_permutations)
]

# Group-level p-value
p_val = (np.sum(np.array(group_null) >= group_obs) + 1) / (len(group_null) + 1)

print(f"Group-level statistic = {group_obs:.3f}, p = {p_val:.4f}")

Group-level statistic = 0.078, p = 0.2827