dspm_localization/causation.py at main · InatiLab/dspm_localization · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Sep 20 15:18:23 2024
@author: Leela Srinivasan

Functions to assess Granger causality of time series data
"""

#Imports
import random
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import grangercausalitytests


def granger_causation_matrix(data, variables, test='ssr_chi2test', verbose=False):
    """


    Parameters
    ----------
    data : mne object
        time series data contained in stc.data
    variables : array
        array of strings for causation matrix labels.
    test : string, optional
        test type. The default is 'ssr_chi2test'.
    verbose : Boolean, optional
        The default is False.

    Returns
    -------
    df : df
        causation df/matrix

    """
    df = pd.DataFrame(np.zeros((len(variables), len(variables))), columns=variables, index=variables)
    MAXLAG=3


    for c in df.columns:
        for r in df.index:
            test_result = grangercausalitytests(data[[r, c]], maxlag=MAXLAG, verbose=False)
            p_values = [round(test_result[i+1][0][test][1],4) for i in range(MAXLAG)]
            if verbose: print(f'Y = {r}, X = {c}, P Values = {p_values}')
            min_p_value = np.min(p_values)
            df.loc[r, c] = min_p_value


    #Label causation matrix
    df.columns = [var + '_x' for var in variables]
    df.index = [var + '_y' for var in variables]
    return df


def granger_causation_proportion(wm_df, stc, start_window, end_window):
    """


    Parameters
    ----------
    wm_df : df
        df of sensors involved in proposed white matter info.
    stc : mne object
        source timecourse.
    start_window : int
        start window for AUC integration.
    end_window : int
        end window for AUC integration.

    Returns
    -------
    granger_causal_proportion : int
        proportion of white matter connections that are Granger causal

    """
    timeseries_df = pd.DataFrame()
    NON_CASUAL = 0
    for row_num in range(wm_df.shape[0]):
        timeseries_df['source_stc'] = stc.data[wm_df['Source Sensor'][row_num], start_window:end_window]
        timeseries_df['sink_stc'] = stc.data[wm_df['Destination Sensor'][row_num], start_window:end_window]


        mat = granger_causation_matrix(timeseries_df, variables = timeseries_df.columns)
        causality = mat.loc['sink_stc_y', 'source_stc_x']
        if causality > 0.00001:
            NON_CASUAL += 1


    granger_causal_proportion = 1-(NON_CASUAL/wm_df.shape[0])
    return granger_causal_proportion


def chance_causation_proportion(stc, start_window, end_window):
    """


    Parameters
    ----------
    stc : mne object
        source timecourse.
    start_window : int
        start window for AUC integration.
    end_window : int
        end window for AUC integration.

    Returns
    -------
    prop : int
        proportion of chance pairings that are Granger Causal

    """
    #Initialize empty variables
    chance_prop=[]
    chance_timeseries_df=pd.DataFrame()


    #Set params
    SOURCE_VS=5124
    SIG_CUTOFF=0.00001
    CHANCE_ITER=1000

    #Iterate through desired number to test chance
    for iter in np.arange(CHANCE_ITER):


        #Randomly select virtual sensor timecourse
        sensor1 =random.randint(0,SOURCE_VS-1)
        sensor2 =random.randint(0,SOURCE_VS-1)
        if sensor1 != sensor2:
            chance_timeseries_df['source_stc'] = stc.data[sensor1, start_window:end_window]
            chance_timeseries_df['sink_stc'] = stc.data[sensor2, start_window:end_window]


            #Compute chance and causality in both directions
            chance_mat = granger_causation_matrix(chance_timeseries_df, variables = chance_timeseries_df.columns)
            chance_xy_causality = chance_mat.loc['sink_stc_y', 'source_stc_x']
            chance_yx_causality = chance_mat.loc['source_stc_y', 'sink_stc_x']


            #Append
            if chance_xy_causality > SIG_CUTOFF or chance_yx_causality > SIG_CUTOFF:
                chance_prop.append(0)
            else:
                chance_prop.append(1)


    #Calculate proportion of 1000 iterations
    prop=chance_prop.count(1)/(chance_prop.count(1)+chance_prop.count(0))
    return prop