Dominance analyses

[hellenmir]

Hello,

Really amazing work and code! Thanks a lot for making it accessible I just had a question on the dominance analyses in the disease.py script. How did you calculate the significant dominance (pspin < 0.05) values? And do you have any code for that? Apologies if I missed it somewhere

Thank you!

[justinehansen]

Hello! Sorry I did not see this until now!

I calculate the null distribution by rerunning the analysis but rotating the brain regions of the variable being predicted (Y) and saving the null dominances (which make up the null distribution for each dominance element). Here's the code I use (I parallelize the process because otherwise it can take a long time):

from joblib import Parallel, delayed
from netneurotools import stats, datasets, utils
from scipy.stats import zscore
import numpy as np
import pandas as pd

def run_dominance_null(X, y, spins, nspins):
    dominance_null = np.zeros((X.shape[1], nspins))
    for k in range(nspins):
        # autorad method
        # notnans = np.invert(spins[:, k] == 33)
        # s = spins[:, k]
        # s = s[:-1]
        # notnans = np.invert(s == 33)
        # s = s[notnans]
        mm, _ = stats.get_dominance_stats(X, y[spins[:, k]])
        dominance_null[:, k] = mm['total_dominance']
    return dominance_null

scale = 'scale033'

cammoun = datasets.fetch_cammoun2012()
info = pd.read_csv(cammoun['info'])
cortex = info.query('scale == @scale & structure == "cortex"')['id']
cortex = np.array(cortex) - 1  # python indexing
nnodes = len(cortex)
hemiid = np.array(info.query('scale == @scale')['hemisphere'])
hemiid = hemiid == 'R'
coords = utils.get_centroids(cammoun[scale], image_space=True)
coords = coords[cortex, :]
spins = stats.gen_spinsamples(coords, hemiid[cortex], seed=1234)
nspins = 1000

power = np.genfromtxt('power_'+scale+'.csv', delimiter=',')
receptor_data = np.genfromtxt('receptor_data_'+scale+'.csv', delimiter=',')
# autorad_data = np.load('/home/jhansen/run_dominance_null/autorad_data.npy')

X = zscore(receptor_data)
y = zscore(power)  # or the disease maps 

null = Parallel(n_jobs=36)(delayed(run_dominance_null)(X, y[:, i], spins, nspins) for i in range(y.shape[1]))
np.save('results/dominance_null_power.npy', np.stack(null, axis=2))

I hope that helps! Let me know if you have any questions.

Justine

[hellenmir]

Great - Thank you very much for this!