Dyadic Representational Similarity Analysis (DyRSA)¶
The following code implements the method described in:
Rhoads, S. A. & Marsh, A. A. (2023). Modeling the minds of friends: Dorsomedial prefrontal cortex encodes accurate dyadic meta-evaluations about how others perceive us. PsyArXiv. https://doi.org/10.31234/osf.io/g9e3z
In [ ]:
%matplotlib inline
%reload_ext autoreload
%autoreload 2
In [ ]:
!pip install nltools nilearn rsatoolbox
import numpy as np, pickle
import matplotlib.pyplot as plt, seaborn as sns
from scipy.stats import spearmanr
from urllib.request import urlopen
import nibabel as nib
from nilearn.image import new_img_like
from nilearn.plotting import view_img_on_surf
from nltools.stats import fdr
from joblib import Parallel, delayed
from tqdm import tqdm
from copy import deepcopy
In [ ]:
!git clone https://github.com/shawnrhoads/dyrsa.git
from dyrsa.code import dyrsa
In [ ]:
# Load participant information
partner_dict = {'sub-1001':'sub-2001','sub-2001':'sub-1001',
'sub-1002':'sub-2002','sub-2002':'sub-1002',
'sub-1003':'sub-2003','sub-2003':'sub-1003',
'sub-1004':'sub-2004','sub-2004':'sub-1004',
'sub-1005':'sub-2005','sub-2005':'sub-1005',
'sub-1006':'sub-2006','sub-2006':'sub-1006',
'sub-1007':'sub-2007','sub-2007':'sub-1007',
'sub-1008':'sub-2008','sub-2008':'sub-1008',
'sub-1009':'sub-2009','sub-2009':'sub-1009',
'sub-1010':'sub-2010','sub-2010':'sub-1010',
'sub-1011':'sub-2011','sub-2011':'sub-1011',
'sub-1012':'sub-2012','sub-2012':'sub-1012',
'sub-1013':'sub-2013','sub-2013':'sub-1013',
'sub-1015':'sub-2015','sub-2015':'sub-1015'}
dyad_ids = [1, 1, 2, 2, 3, 3, 4, 4,
5, 5, 6, 6, 7, 7, 8, 8,
9, 9, 10, 10, 11, 11,
12, 12, 13, 13, 14, 14]
subject_ids = ['sub-1001','sub-2001',
'sub-1002','sub-2002',
'sub-1003','sub-2003',
'sub-1004','sub-2004',
'sub-1005','sub-2005',
'sub-1006','sub-2006',
'sub-1007','sub-2007',
'sub-1008','sub-2008',
'sub-1009','sub-2009',
'sub-1010','sub-2010',
'sub-1011','sub-2011',
'sub-1012','sub-2012',
'sub-1013','sub-2013',
'sub-1015','sub-2015']
nsubj = len(subject_ids)
In [ ]:
# create friend matrix and nonfriend_matrix
friend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
friend_matrix[:] = np.nan
nonfriend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
nonfriend_matrix[:] = np.nan
for i, sub_i in enumerate(subject_ids):
for j, sub_j in enumerate(subject_ids):
if (sub_i == partner_dict[sub_j]) or (partner_dict[sub_i] == sub_j):
friend_matrix[i, j] = 1
friend_matrix[j, i] = 1
else:
nonfriend_matrix[i, j] = 1
nonfriend_matrix[j, i] = 1
# also exclude self
if i == j:
nonfriend_matrix[i, j] = np.nan
friend_matrix[i, j] = np.nan
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))
sns.heatmap(friend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[0], cbar=False, cmap='RdBu')
axs[0].set_title('Friend Pairs')
sns.heatmap(nonfriend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[1], cbar=False, cmap='RdBu')
axs[1].set_title('Non-friend Pairs')
plt.show()
In [ ]:
# Load dictionary of dictionary of searchlight objects contains RDMs based on the cross-validated Euclidean distance between the neural patterns across trials and runs
# Use `rsatoolbox.util.searchlight.get_searchlight_RDMs` to generate each searchlight object (https://rsatoolbox.readthedocs.io/en/latest/rsatoolbox.util.searchlight.html#rsatoolbox.util.searchlight.get_searchlight_RDMs)
# `dyad_rdms` is formatted like: dyad_rdms[sub_id][cond_id]
# e.g., {'sub-1001' : {'Meta': searchlight_object, 'Other': searchlight_object},
# 'sub-1002' : {'Meta': searchlight_object, 'Other': searchlight_object},
# ...}
# pickled file is available at https://osf.io/t5cmd/
try:
with open('../data/dyad_rdms.pickle', 'rb') as handle:
dyad_rdms = pickle.load(handle)
except:
# load from OSF
dyad_rdms = pickle.load(urlopen('https://osf.io/download/yctjb/'))
In [ ]:
# Load mask
mask_fn = 'dyrsa/code/MNI152NLin2009c-Schaefer2018-DMN21-k100_3mm.nii.gz'
mask_file = nib.load(mask_fn)
mask_img = mask_file.get_fdata()
mask = mask_img==1
x, y, z = mask_img.shape
In [ ]:
# Run DyRSA for every possible pair, then can filter friend vs non-friend pairs later
# first grab nvoxels from first subject
voxel_indices = list(dyad_rdms['sub-1001']['Meta'].rdm_descriptors['voxel_index'])
nvoxels = len(voxel_indices)
# initialize result matrices
pairwise_zmat_brain = np.empty((nvoxels, nsubj, nsubj))
pairwise_zmat_brain[:] = np.nan
# loop through all pairs and compute DyRSA
for subject_index, subject_id in tqdm(enumerate(subject_ids)):
for subject_index2, subject_id2 in enumerate(subject_ids):
if subject_index == subject_index2:
pairwise_zmat_brain[:, subject_index, subject_index2] = np.nan
else:
# Neural Meta(A(B(A))) ~ Neural Other(B(A))
subj_sl_RDM = dyad_rdms[subject_id]['Meta']
partner_sl_RDM = dyad_rdms[subject_id2]['Other'] #can also replace with other model RDMs
slresults = Parallel(n_jobs=-1)(delayed(spearmanr)(partner_sl_RDM.dissimilarities[x_index],
x.dissimilarities[0,]) for x_index, x in enumerate(subj_sl_RDM))
spearmanrhos_brain = [float(e[0]) for e in slresults]
pairwise_zmat_brain[:, subject_index, subject_index2] = np.array([np.arctanh(r) for r in deepcopy(spearmanrhos_brain)])
In [ ]:
# remove friend values and flatten across voxels to shape nvox, sum(nonfriend_matrix)=728
all_nonfriend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0],
int(np.sum(nonfriend_matrix==1))))
# remove non-friend values and flatten across voxels to shape nvox, nsubj
all_friend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0], nsubj))
# loop through voxels
for v in range(pairwise_zmat_brain.shape[0]):
all_nonfriend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][nonfriend_matrix==1]
all_friend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][friend_matrix==1]
In [ ]:
# compute the mean of all_friend_pairs_flat
actual_group_mean = np.nanmean(all_friend_pairs_flat, axis=1) # mean at every voxel
# generate null distribution of bootstrap samples
nboot = 5000
zstat, perm_p = dyrsa.run_permutations(nsubj, nboot,
actual_group_mean,
all_nonfriend_pairs_flat,
replacement=True,
write_zstat=True)
In [ ]:
# Write unthresholded result
results_brain = np.zeros((x*y*z))
results_brain[voxel_indices] = zstat
results_brain = results_brain.reshape((x, y, z))
unthresh_img = new_img_like(mask_file, results_brain)
In [ ]:
# interactive plot
view_img_on_surf(unthresh_img,threshold=1)
Out[ ]:
In [ ]:
# Apply FDR correction and write thresholded result
mask_pval = [pv for mask_ind, pv in zip(mask.ravel(), perm_p) if mask_ind]
fdr_thresh = fdr(np.array(mask_pval), q=.05)
results_brain2 = np.zeros((x*y*z))
results_brain2[voxel_indices] = np.array([m if p<fdr_thresh else np.nan for m, p in zip(zstat, perm_p)])
results_brain2 = results_brain2.reshape((x, y, z))
thresh_img = new_img_like(mask_file, results_brain2)
In [ ]:
# interactive plot
view_img_on_surf(thresh_img,threshold=.001)
Out[ ]:
In [ ]: