OSCR

A multidimensional framework for dissociating the neuroplasticity of auditory and early language deprivation.

Code ↔ Paper

10 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 10 matches
  1. [1] § Materials and methods › Cluster analysis of intrinsic functional dimensions across cortical networks ↔ Code/S03_stats_the_rsFC_gradient_2group_model_AD_effect.ipynb, lines 508 to 542 · score 0.79 · co occurrence matrix, n_components, silhouette scoring, Ward, linkage, PCA · Report a mismatch
  2. [2] § Materials and methods › Cluster analysis of intrinsic functional dimensions across cortical networks ↔ Code/S04_stats_the_rsFC_gradient_2group_model_ELD_effect.ipynb, lines 509 to 543 · score 0.79 · co occurrence matrix, n_components, silhouette scoring, Ward, linkage, PCA · Report a mismatch
  3. [3] § Materials and methods › Reconstruction of brain activity with intrinsic functional dimensions ↔ Code/S05_the_rsFC_gradient_rest_predict_task_activation.ipynb, lines 348 to 461 · score 0.79 · predicted task activation, Pearson correlation, FC gradient features, 1–10, incrementally, regression · Report a mismatch
  4. [4] § Materials and methods › Reconstruction of brain activity with intrinsic functional dimensions ↔ Code/S05_the_rsFC_gradient_rest_predict_task_activation.ipynb, lines 348 to 461 · score 0.70 · linear regression model, task activation, 1–10, FC gradients, predict, scored · Report a mismatch
  5. [5] § Materials and methods › Building ridge classifier models based on intrinsic functional dimensions ↔ Code/S02_Calculate_the_rsFC_gradient_2group_classify_model.ipynb, lines 198 to 316 · score 0.66 · Ridge classifier, linear explainer, trained, SHAP, score, ROI · Report a mismatch
  6. [6] § Materials and methods › Computation of the gradients of the resting-state FC ↔ Code/S01_Calculate_the_resting_state_gradient_2group.ipynb, lines 449 to 470 · score 0.61 · Procrustes alignment, FC gradients, template, deaf · Report a mismatch
  7. [7] § Results › The effects of auditory and early language deprivation exhibit dissociation between the dominant and subdominant dimensions › Global feature analysis of the dominant and subdominant dimensio ↔ Code/S04_stats_the_rsFC_gradient_2group_model_ELD_effect.ipynb, lines 1746 to 1883 · score 0.56 · G1 G3, G4 G10, Dom, Sub, LAN, ELD · Report a mismatch
  8. [8] § Results › The effects of auditory and early language deprivation exhibit dissociation between the dominant and subdominant dimensions › Global feature analysis of the dominant and subdominant dimensio ↔ Code/S03_stats_the_rsFC_gradient_2group_model_AD_effect.ipynb, lines 1430 to 1568 · score 0.55 · G1 G3, G4 G10, Dom, Sub, AD, clusters · Report a mismatch
  9. [9] § Materials and methods › Volumetric-to-surface mapping ↔ Code/S05_the_rsFC_gradient_rest_predict_task_activation.ipynb, lines 16 to 47 · score 0.51 · fs lr, Geometry, lh, rh, surf, templates · Report a mismatch
  10. [10] § Materials and methods › Volumetric-to-surface mapping ↔ Code/S01_Calculate_the_resting_state_gradient_2group.ipynb, lines 16 to 47 · score 0.51 · fs lr, Geometry, lh, rh, surf, templates · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,605 lines · 56 KB · no license · 3 matches

  1. # %%
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. import pandas as pd
  5. from scipy.stats import zscore
  6. from sklearn.metrics.pairwise import cosine_similarity
  7. import seaborn as sns
  8. from collections import Counter
  9. import ActflowToolbox as actflow
  11. # %%
  12. ##加载cole模板网络标签
  13. roi_path = 'F:/early_language__experienced_data/经验期待和经验依赖分析/ROI/ColeAnticevicNetPartition/'
  14. labels = pd.read_excel(roi_path + 'CAB-NP_v1.1_Labels-ReorderedbyNetworks.xlsx')['NETWORKKEY'][0:360]
  16. # %%
  17. import h5py, sys, io, numpy as np
  18. import nibabel as nib
  20. path = 'F:/early_language__experienced_data/经验期待和经验依赖分析/功能梯度/code/gradient_analysis-master/gradient_data/templates/'
  22. surfmL = nib.freesurfer.read_geometry(path + 'S900.L.midthickness_MSMAll.32k_fs_LR.surf')
  23. surfiL = nib.freesurfer.read_geometry(path + 'S900.L.very_inflated_MSMAll.32k_fs_LR.surf')
  24. surfL = []
  25. surfL.append(np.array(surfmL[0]*0.3 + surfiL[0]*0.7))
  26. surfL.append(surfmL[1])
  28. surfmR = nib.freesurfer.read_geometry(path + 'S900.R.midthickness_MSMAll.32k_fs_LR.surf')
  29. surfiR = nib.freesurfer.read_geometry(path + 'S900.R.very_inflated_MSMAll.32k_fs_LR.surf')
  30. surfR = []
  31. surfR.append(np.array(surfmR[0]*0.3 + surfiR[0]*0.7))
  32. surfR.append(surfmR[1])
  34. res = nib.load(path +'hcp.tmp.lh.dscalar.nii').get_fdata()
  35. cortL = np.squeeze(np.array(np.where(res != 0)[1], dtype=np.int32))
  36. res = nib.load(path +'hcp.tmp.rh.dscalar.nii').get_fdata()
  37. cortR = np.squeeze(np.array(np.where(res != 0)[1], dtype=np.int32))
  38. cortLen = len(cortL) + len(cortR)
  39. del res
  41. sulcL = np.zeros(len(surfL[0]))
  42. sulcR = np.zeros(len(surfR[0]))
  43. #MSMAll是一种表皮功能对齐的工具
  44. sulcL[cortL] = -1 * nib.load(path + 'S1200.sulc_MSMAll.32k_fs_LR.dscalar.nii').get_fdata()[0,:len(cortL)]
  45. sulcR[cortR] = -1 * nib.load(path + 'S1200.sulc_MSMAll.32k_fs_LR.dscalar.nii').get_fdata()[0,len(cortL)::]
  46. sulcL[np.setdiff1d(range(32492),cortL)] = -1 * nib.load(path + 'Q1-Q6_R440.sulc.32k_fs_LR.dscalar.nii').get_fdata()[0,np.setdiff1d(range(32492),cortL)]
  47. sulcR[np.setdiff1d(range(32492),cortR)] = -1 * nib.load(path + 'Q1-Q6_R440.sulc.32k_fs_LR.dscalar.nii').get_fdata()[0,32492+np.setdiff1d(range(32492),cortR)]
  49. # %% [markdown]
  50. # # 绘制语言激活图
  52. # %%
  53. # 静息态平均梯度画图
  54. from brainspace.gradient import ProcrustesAlignment
  55. import numpy as np
  56. import pickle
  57. import matplotlib.pyplot as plt
  58. from neuromaps.datasets import fetch_fslr
  59. from surfplot import Plot
  60. from brainspace.datasets import load_parcellation
  61. from nilearn.plotting.cm import _cmap_d as nilearn_cmaps
  62. from matplotlib.colors import LinearSegmentedColormap,ListedColormap
  64. # Load and prepare task
  65. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  66. path = 'J:/HCP_analysis/DEAF/data_DF_cole/BLOCK/'
  67. task = np.load(path + 'combined_all_tasks.npy') #91,360,8non language
  68. task_mean = np.mean(task[61:, :, :8], axis = (0,2))
  69. task_mean = zscore(task_mean)
  70. emb = task_mean[:,np.newaxis]
  71. ##加载cole模板
  72. surfaces = fetch_fslr()
  73. lh, rh = surfaces['veryinflated']
  74. p = Plot(lh, rh,brightness=.6)
  76. Cole_path = 'F:/early_language__experienced_data/经验期待和经验依赖分析/ROI/ColeAnticevicNetPartition/'
  77. coleLR = nib.load(Cole_path+'Glasser360Indices_LR.dscalar.nii').get_fdata()[:, 0:64984]
  78. #把模板填入到完整的皮层中
  79. for i in range(360):
  80. coleLR[coleLR==i+1] = emb[i,0]
  81. cifti = coleLR.T
  82. dataL = np.zeros((len(surfL[0]), cifti.shape[1]))
  83. dataL[cortL,:] = cifti[0:len(cortL),:]
  84. dataR = np.zeros((len(surfR[0]), cifti.shape[1]))
  85. dataR[cortR,:] = cifti[len(cortL):cortLen,:]
  87. cmap = nilearn_cmaps['brown_cyan']
  88. # 翻转色阶,将深色设置在两端
  89. colors = cmap(np.linspace(0, 1, cmap.N))
  90. colors = np.concatenate([colors[-int(cmap.N/2):], colors[:int(cmap.N/2)]], axis=0)
  91. new_cmap = LinearSegmentedColormap.from_list('custom', colors)
  92. p.add_layer({'left': dataL, 'right': dataR}, cmap=new_cmap, cbar=True)
  93. fig = p.build()
  94. plt.savefig(out_path + 'task_language_activation_NON-ELA.png', dpi=300, bbox_inches='tight')
  95. fig.show()
  97. # %%
  98. # 静息态平均梯度画图
  99. from brainspace.gradient import ProcrustesAlignment
  100. import numpy as np
  101. import pickle
  102. import matplotlib.pyplot as plt
  103. from neuromaps.datasets import fetch_fslr
  104. from surfplot import Plot
  105. from brainspace.datasets import load_parcellation
  106. from nilearn.plotting.cm import _cmap_d as nilearn_cmaps
  107. from matplotlib.colors import LinearSegmentedColormap,ListedColormap
  109. # Load and prepare task
  110. path = 'J:/HCP_analysis/DEAF/1_gradient/'
  111. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  112. emb_all = np.load(path + 'restFC_align_Gradients_10.npy') # emb_all.shape(91, 360, 10)
  114. path = 'J:/HCP_analysis/DEAF/data_DF_cole/BLOCK/'
  115. task1 = np.load(path + 'combined_all_tasks.npy')[:, :, :8] #91,360,8 language
  116. task1 = np.mean(task1, axis = 2)
  117. task2 = np.load(path + 'combined_all_tasks.npy')[:, :, 8] #91,360,1 non language
  119. task = task1 - task2
  120. task_mean = np.mean(task[:34], axis = 0)
  121. emb = task_mean[:,np.newaxis]
  122. ##加载cole模板
  123. surfaces = fetch_fslr()
  124. lh, rh = surfaces['veryinflated']
  125. p = Plot(lh, rh,brightness=.6)
  127. Cole_path = 'F:/early_language__experienced_data/经验期待和经验依赖分析/ROI/ColeAnticevicNetPartition/'
  128. coleLR = nib.load(Cole_path+'Glasser360Indices_LR.dscalar.nii').get_fdata()[:, 0:64984]
  129. #把模板填入到完整的皮层中
  130. for i in range(360):
  131. coleLR[coleLR==i+1] = emb[i,0]
  132. cifti = coleLR.T
  133. dataL = np.zeros((len(surfL[0]), cifti.shape[1]))
  134. dataL[cortL,:] = cifti[0:len(cortL),:]
  135. dataR = np.zeros((len(surfR[0]), cifti.shape[1]))
  136. dataR[cortR,:] = cifti[len(cortL):cortLen,:]
  138. cmap = nilearn_cmaps['brown_cyan']
  139. # 翻转色阶,将深色设置在两端
  140. colors = cmap(np.linspace(0, 1, cmap.N))
  141. colors = np.concatenate([colors[-int(cmap.N/2):], colors[:int(cmap.N/2)]], axis=0)
  142. new_cmap = LinearSegmentedColormap.from_list('custom', colors)
  143. p.add_layer({'left': dataL, 'right': dataR}, cmap=new_cmap, cbar=True)
  144. fig = p.build()
  145. plt.savefig(out_path + 'task_language_activation_hearing_reverse.png', dpi=300, bbox_inches='tight')
  146. fig.show()
  148. # %%
  149. # 静息态平均梯度画图
  150. from brainspace.gradient import ProcrustesAlignment
  151. import numpy as np
  152. import pickle
  153. import matplotlib.pyplot as plt
  154. from neuromaps.datasets import fetch_fslr
  155. from surfplot import Plot
  156. from brainspace.datasets import load_parcellation
  157. from nilearn.plotting.cm import _cmap_d as nilearn_cmaps
  158. from matplotlib.colors import LinearSegmentedColormap,ListedColormap
  160. # Load and prepare task
  161. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  162. predict_y = np.load(out_path + 'Y_pred_gradient_predict_task.npy') # emb_all.shape(91, 360, 10)
  163. predict_y_mean = np.mean(predict_y[:34], axis = (0,1))
  164. emb = predict_y_mean[:,np.newaxis]
  165. ##加载cole模板
  166. surfaces = fetch_fslr()
  167. lh, rh = surfaces['veryinflated']
  168. p = Plot(lh, rh,brightness=.6)
  170. Cole_path = 'F:/early_language__experienced_data/经验期待和经验依赖分析/ROI/ColeAnticevicNetPartition/'
  171. coleLR = nib.load(Cole_path+'Glasser360Indices_LR.dscalar.nii').get_fdata()[:, 0:64984]
  172. #把模板填入到完整的皮层中
  173. for i in range(360):
  174. coleLR[coleLR==i+1] = emb[i,0]
  175. cifti = coleLR.T
  176. dataL = np.zeros((len(surfL[0]), cifti.shape[1]))
  177. dataL[cortL,:] = cifti[0:len(cortL),:]
  178. dataR = np.zeros((len(surfR[0]), cifti.shape[1]))
  179. dataR[cortR,:] = cifti[len(cortL):cortLen,:]
  181. cmap = nilearn_cmaps['brown_cyan']
  182. # 翻转色阶,将深色设置在两端
  183. colors = cmap(np.linspace(0, 1, cmap.N))
  184. colors = np.concatenate([colors[-int(cmap.N/2):], colors[:int(cmap.N/2)]], axis=0)
  185. new_cmap = LinearSegmentedColormap.from_list('custom', colors)
  186. p.add_layer({'left': dataL, 'right': dataR}, cmap=new_cmap, cbar=True)
  187. fig = p.build()
  188. plt.savefig(out_path + 'task_language_predict_activation_hearing.png', dpi=300, bbox_inches='tight')
  189. fig.show()
  191. # %% [markdown]
  192. # # 使用静息态梯度的神经可塑性重构实际任务的激活差值
  194. # %%
  195. import numpy as np
  196. import shap
  197. from sklearn.linear_model import Ridge
  198. from scipy.stats import pearsonr
  199. import matplotlib.pyplot as plt
  200. from sklearn.linear_model import LinearRegression
  201. from scipy.stats import zscore
  203. # Load and prepare task
  204. path = 'J:/HCP_analysis/DEAF/1_gradient/'
  205. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  206. emb_all = np.load(path + 'restFC_align_Gradients_10.npy') # emb_all.shape(91, 360, 10)
  208. path = 'J:/HCP_analysis/DEAF/data_DF_cole/BLOCK/'
  209. task_combined = np.load(path + 'combined_all_tasks.npy')[:, :, :8] #91,360,4 non language
  211. subs = 91
  212. features = 10
  213. blocks = 8
  214. tasks = 2
  215. rois = 360
  216. # Initialize result arrays
  217. corr_predict_ind = np.zeros((subs, tasks, blocks, 3))
  218. Y_pred_all = np.zeros((subs, tasks, blocks, rois))
  219. shap_values_all = np.zeros((subs, tasks, blocks, features)) # For storing SHAP values
  221. for d, task in enumerate(task_combined):
  222. for sub in range(subs):
  223. for block in range(blocks):
  224. # Prepare data
  225. X = zscore(emb_all[sub], axis=0) # Resting-state gradients (360, 10)
  226. Y = zscore(task[sub, :, block]) # Task activation (360,)
  228. # Fit linear regression model
  229. model = LinearRegression()
  230. model.fit(X, Y)
  232. # Compute predictions
  233. Y_pred = model.predict(X)
  234. Y_pred_all[sub, :, d] = Y_pred
  236. # Calculate statistics
  237. r_squared = model.score(X, Y)
  238. corr, p_value = pearsonr(Y, Y_pred)
  240. # Store results
  241. corr_predict_ind[sub, d, block, 0] = corr
  242. corr_predict_ind[sub, d, block, 1] = p_value
  243. corr_predict_ind[sub, d, block, 2] = r_squared
  245. # SHAP analysis
  246. explainer = shap.Explainer(model.predict, X)
  247. shap_values = explainer(X)
  249. # Average absolute SHAP values across ROIs
  250. shap_values_all[sub, d, block] = np.mean(np.abs(shap_values.values), axis=0)
  252. # Save results
  253. np.save(out_path + 'gradient_predict_task.npy', corr_predict_ind)
  254. np.save(out_path + 'Y_pred_gradient_predict_task.npy', Y_pred_all)
  255. np.save(out_path + 'shap_feature_importance.npy', shap_values_all)
  257. # Calculate and save mean SHAP values across subjects and tasks
  258. mean_shap = np.mean(shap_values_all, axis=(0, 1, 2))
  259. np.save(out_path + 'mean_shap_feature_importance.npy', mean_shap)
  261. # Plot mean SHAP importance
  262. plt.figure(figsize=(10, 6))
  263. plt.bar(range(features), mean_shap)
  264. plt.title('Average SHAP Feature Importance')
  265. plt.xlabel('Gradient Feature')
  266. plt.ylabel('Mean |SHAP value|')
  267. plt.savefig(out_path + 'mean_shap_feature_importance.png')
  268. plt.close()
  270. print("Analysis completed and results saved.")
  272. # %%
  273. import numpy as np
  274. import shap
  275. from sklearn.linear_model import Ridge
  276. from scipy.stats import pearsonr
  277. import matplotlib.pyplot as plt
  278. from sklearn.linear_model import LinearRegression
  279. from scipy.stats import zscore
  281. # Load and prepare task
  282. path = 'J:/HCP_analysis/DEAF/1_gradient/'
  283. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  284. emb_all = np.load(path + 'restFC_align_Gradients_10.npy') # emb_all.shape(91, 360, 10)
  286. path = 'J:/HCP_analysis/DEAF/data_DF_cole/BLOCK/'
  287. task = np.load(path + 'combined_all_tasks.npy')[:, :, :8] #91,360,8 non language
  289. subs = 91
  290. features = 10
  291. rois = 360
  292. blocks = 8
  293. # Initialize result arrays
  294. corr_predict_ind = np.zeros((subs, blocks, 3))
  295. Y_pred_all = np.zeros((subs, blocks, rois))
  296. shap_values_all = np.zeros((subs, blocks, features)) # For storing SHAP values
  298. for sub in range(subs):
  299. for block in range(blocks):
  300. # Prepare data
  301. X = zscore(emb_all[sub], axis=0) # Resting-state gradients (360, 10)
  302. Y = zscore(task[sub, :, block]) # Task activation (360,)
  304. # Fit linear regression model
  305. model = LinearRegression()
  306. model.fit(X, Y)
  308. # Compute predictions
  309. Y_pred = model.predict(X)
  310. Y_pred_all[sub, block, :] = Y_pred
  312. # Calculate statistics
  313. r_squared = model.score(X, Y)
  314. corr, p_value = pearsonr(Y, Y_pred)
  316. # Store results
  317. corr_predict_ind[sub, block, 0] = corr
  318. corr_predict_ind[sub, block, 1] = p_value
  319. corr_predict_ind[sub, block, 2] = r_squared
  321. # SHAP analysis
  322. explainer = shap.Explainer(model.predict, X)
  323. shap_values = explainer(X)
  325. # Average absolute SHAP values across ROIs
  326. shap_values_all[sub, block] = np.mean(np.abs(shap_values.values), axis=0)
  328. # Save results
  329. np.save(out_path + 'gradient_predict_task.npy', corr_predict_ind)
  330. np.save(out_path + 'Y_pred_gradient_predict_task.npy', Y_pred_all)
  331. np.save(out_path + 'shap_feature_importance.npy', shap_values_all)
  333. # Calculate and save mean SHAP values across subjects and tasks
  334. mean_shap = np.mean(shap_values_all, axis=(0, 1, 2))
  335. np.save(out_path + 'mean_shap_feature_importance.npy', mean_shap)
  337. # Plot mean SHAP importance
  338. plt.figure(figsize=(10, 6))
  339. plt.bar(range(features), mean_shap)
  340. plt.title('Average SHAP Feature Importance')
  341. plt.xlabel('Gradient Feature')
  342. plt.ylabel('Mean |SHAP value|')
  343. plt.savefig(out_path + 'mean_shap_feature_importance.png')
  344. plt.close()
  346. print("Analysis completed and results saved.")
  348. # %% [markdown]
  349. # # 循环增加特征
  351. # %%
  352. import numpy as np
  353. import shap
  354. from sklearn.linear_model import Ridge
  355. from scipy.stats import pearsonr
  356. import matplotlib.pyplot as plt
  357. from sklearn.linear_model import LinearRegression
  358. from scipy.stats import zscore
  360. # Load and prepare task
  361. path = 'J:/HCP_analysis/DEAF/1_gradient/'
  362. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  363. emb_all = np.load(path + 'restFC_align_Gradients_10.npy') # emb_all.shape(91, 360, 10)
  365. path = 'J:/HCP_analysis/DEAF/data_DF_cole/BLOCK/'
  366. task = np.load(path + 'combined_all_tasks.npy')[:, :, :8] #91,360,4 non language
  368. subs = 91
  369. features = 10
  370. blocks = 8
  371. rois = 360 # Added this since it was used but not defined
  373. # Initialize result arrays for all feature combinations
  374. # We'll have results for each feature count (1-10), each subject, each block
  375. corr_predict_ind = np.zeros((features, subs, blocks, 3)) # Now shape (10, 91, 4, 3)
  376. Y_pred_all = np.zeros((features, subs, blocks, rois)) # Now shape (10, 91, 4, 360)
  377. shap_values_all = np.zeros((features, subs, blocks, features)) # Now shape (10, 91, 4, 10)
  379. for num_features in range(1, features+1): # Loop from 1 to 10 features
  380. print(f"Processing models with {num_features} features...")
  382. for sub in range(subs):
  383. for block in range(blocks):
  384. # Prepare data - use only the first num_features features
  385. X = zscore(emb_all[sub, :, :num_features], axis=0) # (360, num_features)
  386. Y = zscore(task[sub, :, block]) # Task activation (360,)
  388. # Fit linear regression model
  389. model = LinearRegression()
  390. model.fit(X, Y)
  392. # Compute predictions
  393. Y_pred = model.predict(X)
  394. Y_pred_all[num_features-1, sub, block, :] = Y_pred # num_features-1 for 0-based index
  396. # Calculate statistics
  397. r_squared = model.score(X, Y)
  398. corr, p_value = pearsonr(Y, Y_pred)
  400. # Store results
  401. corr_predict_ind[num_features-1, sub, block, 0] = corr
  402. corr_predict_ind[num_features-1, sub, block, 1] = p_value
  403. corr_predict_ind[num_features-1, sub, block, 2] = r_squared
  405. # SHAP analysis
  406. explainer = shap.Explainer(model.predict, X)
  407. shap_values = explainer(X)
  409. # Store SHAP values (only the first num_features will have values)
  410. shap_values_all[num_features-1, sub, block, :num_features] = np.mean(np.abs(shap_values.values), axis=0)
  412. # Save results
  413. np.save(out_path + 'gradient_predict_task_by_feature_count.npy', corr_predict_ind)
  414. np.save(out_path + 'Y_pred_gradient_predict_task_by_feature_count.npy', Y_pred_all)
  415. np.save(out_path + 'shap_feature_importance_by_feature_count.npy', shap_values_all)
  417. # Calculate and save mean results across subjects and tasks for each feature count
  418. mean_corr = np.mean(corr_predict_ind[:,:,:,0], axis=(1,2)) # Mean correlation for each feature count
  419. mean_r2 = np.mean(corr_predict_ind[:,:,:,2], axis=(1,2)) # Mean R-squared for each feature count
  420. mean_shap = np.zeros((features, features)) # Will store mean SHAP for each feature at each feature count
  422. for num_features in range(features):
  423. # For each feature count, calculate mean SHAP values for the features used
  424. mean_shap[num_features, :num_features+1] = np.mean(shap_values_all[num_features, :, :, :num_features+1], axis=(0,1))
  426. np.save(out_path + 'mean_correlation_by_feature_count.npy', mean_corr)
  427. np.save(out_path + 'mean_r2_by_feature_count.npy', mean_r2)
  428. np.save(out_path + 'mean_shap_by_feature_count.npy', mean_shap)
  430. # Plot results
  431. plt.figure(figsize=(12, 8))
  433. # Plot mean correlation by feature count
  434. plt.subplot(2, 2, 1)
  435. plt.plot(range(1, features+1), mean_corr, marker='o')
  436. plt.title('Mean Prediction Correlation by Number of Features')
  437. plt.xlabel('Number of Features')
  438. plt.ylabel('Mean Pearson Correlation')
  440. # Plot mean R-squared by feature count
  441. plt.subplot(2, 2, 2)
  442. plt.plot(range(1, features+1), mean_r2, marker='o')
  443. plt.title('Mean R-squared by Number of Features')
  444. plt.xlabel('Number of Features')
  445. plt.ylabel('Mean R-squared')
  447. # Plot SHAP importance heatmap
  448. plt.subplot(2, 1, 2)
  449. plt.imshow(mean_shap.T, aspect='auto', cmap='viridis')
  450. plt.colorbar(label='Mean |SHAP value|')
  451. plt.title('SHAP Feature Importance by Number of Features Used')
  452. plt.xlabel('Number of Features Used in Model')
  453. plt.ylabel('Gradient Feature')
  454. plt.yticks(range(features), range(1, features+1))
  455. plt.xticks(range(features), range(1, features+1))
  457. plt.tight_layout()
  458. plt.savefig(out_path + 'feature_incremental_results.png')
  459. plt.close()
  461. print("Analysis completed and results saved.")
  463. # %% [markdown]
  464. # # 断点回归分析
  466. # %%
  467. import numpy as np
  468. import matplotlib.pyplot as plt
  469. from sklearn.linear_model import LinearRegression
  470. from scipy.stats import ttest_1samp, ttest_rel, zscore
  471. import seaborn as sns
  473. # Load data
  474. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  475. corr = np.load(out_path + 'gradient_predict_task_by_feature_count.npy')[:, :, :, 0] # (10, 91, 9)
  477. # Preprocess: Fisher z-transform but don't average across blocks
  478. corr_z = np.arctanh(corr) # (10, 91, 8)
  480. # Define breakpoint (cutoff=3)
  481. cutoff = 2
  482. dimensions = np.arange(10) # 0-9
  484. # Initialize storage for results
  485. n_subjects = corr_z.shape[1]
  486. n_blocks = corr_z.shape[2]
  487. slopes_low = np.zeros((n_subjects, n_blocks)) # Low-dim slopes for each subject and block
  488. slopes_high = np.zeros((n_subjects, n_blocks)) # High-dim slopes
  489. jump_effects = np.zeros((n_subjects, n_blocks)) # Jump effects
  491. for subj in range(n_subjects):
  492. for block in range(n_blocks):
  493. y = corr_z[:, subj, block]
  495. # Fit low-dimension (0-3) regression
  496. X_low = dimensions[dimensions <= cutoff].reshape(-1, 1)
  497. reg_low = LinearRegression().fit(X_low, y[dimensions <= cutoff])
  498. y_pred_low = reg_low.predict([[cutoff]])
  499. slopes_low[subj, block] = reg_low.coef_[0]
  501. # Fit high-dimension (3-9) regression
  502. X_high = dimensions[dimensions >= cutoff].reshape(-1, 1)
  503. reg_high = LinearRegression().fit(X_high, y[dimensions >= cutoff])
  504. y_pred_high = reg_high.predict([[cutoff]])
  505. slopes_high[subj, block] = reg_high.coef_[0]
  507. # Calculate jump effect
  508. jump_effects[subj, block] = y_pred_high - y_pred_low
  510. # =================================================================
  511. # 1. Analyze results across subjects and blocks
  512. # =================================================================
  513. # Flatten the results across blocks for subject-level analysis
  514. slopes_low_flat = slopes_low.ravel()
  515. slopes_high_flat = slopes_high.ravel()
  516. jump_effects_flat = jump_effects.ravel()
  518. # Low-dimension slopes
  519. t_low, p_low = ttest_1samp(slopes_low_flat, 0)
  520. print("=== Low-dimension regression (0-3) ===")
  521. print(f"Mean Slope: {np.mean(slopes_low_flat):.4f} ± {np.std(slopes_low_flat, ddof=1):.4f}")
  522. print(f"t = {t_low:.4f}, p = {p_low:.4f}\n")
  524. # High-dimension slopes
  525. t_high, p_high = ttest_1samp(slopes_high_flat, 0)
  526. print("=== High-dimension regression (3-9) ===")
  527. print(f"Mean Slope: {np.mean(slopes_high_flat):.4f} ± {np.std(slopes_high_flat, ddof=1):.4f}")
  528. print(f"t = {t_high:.4f}, p = {p_high:.4f}\n")
  530. # Slope comparison (paired t-test)
  531. t_slope_diff, p_slope_diff = ttest_rel(slopes_low_flat, slopes_high_flat)
  532. print("=== Low vs High dimension slope comparison ===")
  533. print(f"Mean Slope Difference: {np.mean(slopes_low_flat - slopes_high_flat):.4f}")
  534. print(f"t = {t_slope_diff:.4f}, p = {p_slope_diff:.4f}\n")
  536. # Jump effects
  537. mean_jump = np.mean(jump_effects_flat)
  538. std_jump = np.std(jump_effects_flat, ddof=1)
  539. t_jump, p_jump = ttest_1samp(jump_effects_flat, 0)
  540. print("=== Breakpoint jump effect ===")
  541. print(f"Mean Jump: {mean_jump:.4f} ± {std_jump:.4f}")
  542. print(f"t = {t_jump:.4f}, p = {p_jump:.4f}")
  544. # =================================================================
  545. # 2. Mixed-effects analysis (optional)
  546. # =================================================================
  547. # For a more sophisticated analysis, you could use mixed-effects models
  548. # that account for both within-subject and between-subject variability
  549. # (would require packages like statsmodels or pingouin)
  551. # =================================================================
  552. # 3. Visualization
  553. # =================================================================
  554. plt.figure(figsize=(15, 5))
  556. # Slope distributions
  557. plt.subplot(1, 3, 1)
  558. sns.histplot(slopes_low_flat, color='blue', kde=True, label='Low Dim (0-3)')
  559. sns.histplot(slopes_high_flat, color='orange', kde=True, label='High Dim (3-9)')
  560. plt.axvline(0, color='black', linestyle='--')
  561. plt.xlabel('Regression Slope')
  562. plt.ylabel('Frequency')
  563. plt.title('Slope Distribution by Dimension')
  564. plt.legend()
  566. # Slope differences
  567. plt.subplot(1, 3, 2)
  568. sns.histplot(slopes_low_flat - slopes_high_flat, color='green', kde=True)
  569. plt.axvline(0, color='black', linestyle='--')
  570. plt.xlabel('Slope Difference (Low - High)')
  571. plt.title('Slope Difference Distribution')
  573. # Jump effects
  574. plt.subplot(1, 3, 3)
  575. sns.histplot(jump_effects_flat, color='red', kde=True)
  576. plt.axvline(mean_jump, color='black', linestyle='--', label=f'Mean Jump = {mean_jump:.3f}')
  577. plt.xlabel('Jump Effect at Cutoff')
  578. plt.title('Jump Effect Distribution')
  579. plt.legend()
  581. plt.tight_layout()
  582. plt.show()
  584. # %%
  585. import numpy as np
  586. import shap
  587. from sklearn.linear_model import Ridge
  588. from scipy.stats import pearsonr
  589. import matplotlib.pyplot as plt
  590. from sklearn.linear_model import LinearRegression
  591. from scipy.stats import zscore
  593. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  594. corr = np.load(out_path + 'gradient_predict_task_by_feature_count.npy')[:, :, :, 0] # (10, 91, 4)
  595. # 计算均值和误差（标准差或标准误）
  596. mean_corr = np.mean(corr, axis=(1, 2)) # 跨被试和block的平均
  597. std_corr = np.std(corr, axis=(1, 2)) # 标准差
  599. mean_shap = np.load(out_path + 'mean_shap_by_feature_count.npy')
  600. subs = 91
  601. features = 10
  602. blocks = 8
  603. rois = 360 # Added this since it was used but not defined
  605. # 绘图设置
  606. plt.figure(figsize=(8, 5))
  607. title_font = {'fontsize': 18, 'fontweight': 'bold'}
  608. label_font = {'fontsize': 18, 'fontweight': 'bold'}
  609. tick_font = {'fontsize': 18, 'fontweight': 'bold'}
  611. # 绘制均值线及误差线（这里使用标准误，如需标准差请替换为std_corr）
  612. plt.errorbar(
  613. range(1, 11),
  614. mean_corr,
  615. yerr=std_corr, # 误差线长度（标准误）
  616. fmt='-o', # 线型+标记
  617. linewidth=3,
  618. markersize=8,
  619. color='gray',
  620. capsize=5, # 误差线端帽长度
  621. capthick=2, # 端帽线宽
  622. ecolor='gray' # 误差线颜色
  623. )
  625. # 美化图形
  626. plt.xlabel('Number of Features', **label_font)
  627. plt.ylabel('Correlation (±SD)', **label_font)
  628. plt.xticks(range(1, 11), **tick_font)
  629. plt.yticks(**tick_font)
  630. plt.grid(False)
  632. ax = plt.gca()
  633. ax.spines['top'].set_visible(False)
  634. ax.spines['right'].set_visible(False)
  635. ax.spines['bottom'].set_linewidth(3)
  636. ax.spines['left'].set_linewidth(3)
  637. ax.spines['bottom'].set_color('black')
  638. ax.spines['left'].set_color('black')
  640. plt.tight_layout()
  641. plt.savefig(out_path + 'mean_correlation_plot_with_errorbars.png', dpi=300, bbox_inches='tight')
  642. plt.close()
  644. # Second plot: SHAP Heatmap
  645. plt.figure(figsize=(7, 4))
  646. im = plt.imshow(mean_shap.T, aspect='auto', cmap='viridis')
  647. #plt.title('Feature Importance by Number of Features', **title_font)
  648. plt.xlabel('Number of Features Used in Model', **label_font)
  649. #plt.ylabel('Gradient Feature', **label_font)
  650. plt.yticks(range(features), [f'G{i+1}' for i in range(10)], **tick_font)
  651. plt.xticks(range(features), range(1, features+1), **tick_font)
  653. cbar = plt.colorbar(im, pad=0.05)
  654. cbar.set_label('Mean SHAP Value', **label_font)
  655. cbar.ax.tick_params(
  656. labelsize=15,
  657. which='both',
  658. width=2,
  659. length=6
  660. )
  661. for label in cbar.ax.get_yticklabels():
  662. label.set_fontweight('bold')
  664. plt.tight_layout()
  665. plt.savefig(out_path + 'shap_heatmap.png', dpi=300, bbox_inches='tight')
  666. plt.close()
  668. # %%
  669. import numpy as np
  670. import matplotlib.pyplot as plt
  671. from matplotlib.colors import LinearSegmentedColormap
  672. from scipy.stats import linregress
  674. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  675. corr = np.load(out_path + 'gradient_predict_task_by_feature_count.npy')[:, :, :, 0] # (10, 91, 8)
  677. # Define 10-step red gradient colormap
  678. red_gradient = [
  679. '#400000', # 深黑红
  680. '#800000', # 暗红
  681. '#b30000', # 标准红
  682. '#e60000', # 亮红
  683. '#ff3333', # 浅红
  684. '#ff6666', # 淡红
  685. '#ff9999', # 粉红
  686. '#ffb3b3', # 浅粉
  687. '#ffd6d6', # 极浅粉
  688. '#ffcccc' # 淡灰红
  689. ]
  690. n_dims = 10
  692. # Calculate statistics
  693. mean_corr = np.mean(corr, axis=(1, 2)) # Average across subjects and blocks
  694. std_corr = np.std(corr, axis=(1, 2)) # Standard deviation
  696. # Plot settings
  697. plt.figure(figsize=(8, 5))
  698. label_font = {'fontsize': 24, 'fontweight': 'bold'}
  699. tick_font = {'fontsize': 24, 'fontweight': 'bold'}
  701. # Plot jittered individual points with gradient colors
  702. for i in range(10):
  703. # Create jittered x positions
  704. x_jitter = np.random.normal(i+1, 0.08, size=corr[i].size)
  705. plt.scatter(
  706. x_jitter,
  707. corr[i].flatten(),
  708. color=red_gradient[i],
  709. alpha=0.7,
  710. s=25,
  711. edgecolor='white',
  712. linewidth=0.5,
  713. zorder=1
  714. )
  716. # Add boxplot with thicker elements and black borders
  717. box_data = [corr[i].flatten() for i in range(10)]
  718. box = plt.boxplot(
  719. box_data,
  720. positions=range(1, 11),
  721. widths=0.5,
  722. patch_artist=True,
  723. showfliers=False,
  724. medianprops={'color': 'black', 'linewidth': 3}, # Thicker median line
  725. boxprops={'color': 'black', 'linewidth': 3}, # Black border for box
  726. whiskerprops={'color': 'black', 'linewidth': 3}, # Black whiskers
  727. capprops={'color': 'black', 'linewidth': 3}, # Black caps
  728. zorder=2
  729. )
  731. # Color boxes with red gradient
  732. for patch, color in zip(box['boxes'], red_gradient):
  733. patch.set_facecolor(color)
  734. patch.set_alpha(0.4)
  736. # Calculate and plot regression lines
  737. x_vals = np.array(range(1, 11))
  739. # Regression for features 1-3 (dark red)
  740. x1 = x_vals[:3]
  741. y1 = mean_corr[:3]
  742. slope1, intercept1, _, _, _ = linregress(x1, y1)
  743. plt.plot(x1, intercept1 + slope1*x1,
  744. color='#800000', linewidth=5, linestyle='-',
  745. label='Features 1-3 Trend')
  747. # Regression for features 3-10 (pink)
  748. x2 = x_vals[2:] # From feature 3 to 10
  749. y2 = mean_corr[2:]
  750. slope2, intercept2, _, _, _ = linregress(x2, y2)
  751. plt.plot(x2, intercept2 + slope2*x2,
  752. color='#ff9999', linewidth=5, linestyle='-',
  753. label='Features 3-10 Trend')
  755. # Plot mean points without connecting line
  756. plt.scatter(
  757. range(1, 11),
  758. mean_corr,
  759. s=100, # Larger point size
  760. color='black',
  761. edgecolor='white',
  762. linewidth=2,
  763. zorder=4
  764. )
  766. # Add error bars separately (without connecting line)
  767. plt.errorbar(
  768. range(1, 11),
  769. mean_corr,
  770. fmt='none', # No markers or connecting line
  771. color='black',
  772. capsize=5,
  773. capthick=2,
  774. linewidth=1.5,
  775. zorder=3
  776. )
  778. # Formatting
  779. plt.xlabel('Number of Features', **label_font)
  780. plt.ylabel('Correlation (±SD)', **label_font)
  781. plt.xticks(range(1, 11), **tick_font)
  782. plt.yticks(**tick_font)
  783. plt.grid(False)
  785. ax = plt.gca()
  786. ax.spines['top'].set_visible(False)
  787. ax.spines['right'].set_visible(False)
  788. ax.spines['bottom'].set_linewidth(3)
  789. ax.spines['left'].set_linewidth(3)
  790. ax.spines['bottom'].set_color('black')
  791. ax.spines['left'].set_color('black')
  793. plt.tight_layout()
  794. plt.savefig(out_path + 'feature_correlation_plot.png', dpi=300, bbox_inches='tight', transparent=True)
  795. plt.close()
  797. # %%
  798. import numpy as np
  799. import matplotlib.pyplot as plt
  800. from matplotlib.colors import LinearSegmentedColormap
  801. from scipy.stats import linregress
  803. out_path = 'J:/HCP_analysis/DEAF/10_rest_predict_task/'
  804. corr = np.load(out_path + 'gradient_predict_task_by_feature_count.npy')[:, :, :, 0] # (10, 91, 4)
  806. # Define 10-step red gradient colormap
  807. red_gradient = [
  808. '#400000', # 深黑红
  809. '#800000', # 暗红
  810. '#b30000', # 标准红
  811. '#e60000', # 亮红
  812. '#ff3333', # 浅红
  813. '#ff6666', # 淡红
  814. '#ff9999', # 粉红
  815. '#ffb3b3', # 浅粉
  816. '#ffd6d6', # 极浅粉
  817. '#ffcccc' # 淡灰红
  818. ]
  819. n_dims = 10
  821. # Calculate statistics
  822. mean_corr = np.mean(corr, axis=(1, 2)) # Average across subjects and blocks
  823. std_corr = np.std(corr, axis=(1, 2)) # Standard deviation
  825. # Plot settings
  826. plt.figure(figsize=(8, 5))
  827. label_font = {'fontsize': 24, 'fontweight': 'bold'}
  828. tick_font = {'fontsize': 24, 'fontweight': 'bold'}
  830. # Plot jittered individual points with gradient colors
  831. for i in range(10):
  832. # Create jittered x positions
  833. x_jitter = np.random.normal(i+1, 0.08, size=corr[i].size)
  834. plt.scatter(
  835. x_jitter,
  836. corr[i].flatten(),
  837. color=red_gradient[i],
  838. alpha=0.7,
  839. s=25,
  840. edgecolor='white',
  841. linewidth=0.5,
  842. zorder=1
  843. )
  845. # Plot mean points with error bars showing standard deviation
  846. plt.errorbar(
  847. range(1, 11),
  848. mean_corr,
  849. yerr=std_corr,
  850. fmt='o',
  851. markersize=10,
  852. capsize=5,
  853. capthick=3,
  854. elinewidth=3,
  855. color='black',
  856. ecolor='black',
  857. markerfacecolor='black',
  858. markeredgecolor='white',
  859. markeredgewidth=2,
  860. zorder=4
  861. )
  863. # Calculate and plot regression lines
  864. x_vals = np.array(range(1, 11))
  866. # Regression for features 1-3 (dark red)
  867. x1 = x_vals[:3]
  868. y1 = mean_corr[:3]
  869. slope1, intercept1, _, _, _ = linregress(x1, y1)
  870. plt.plot(x1, intercept1 + slope1*x1,
  871. color='#800000', linewidth=5, linestyle='-',
  872. label='Features 1-3 Trend')
  874. # Regression for features 3-10 (pink)
  875. x2 = x_vals[2:] # From feature 3 to 10
  876. y2 = mean_corr[2:]
  877. slope2, intercept2, _, _, _ = linregress(x2, y2)
  878. plt.plot(x2, intercept2 + slope2*x2,
  879. color='#ff9999', linewidth=5, linestyle='-',
  880. label='Features 3-10 Trend')
  882. # Formatting
  883. plt.xlabel('Number of Features', **label_font)
  884. plt.ylabel('Correlation (±SD)', **label_font)
  885. plt.xticks(range(1, 11), **tick_font)
  886. plt.yticks(**tick_font)
  887. plt.grid(False)
  889. ax = plt.gca()
  890. ax.spines['top'].set_visible(False)
  891. ax.spines['right'].set_visible(False)
  892. ax.spines['bottom'].set_linewidth(3)
  893. ax.spines['left'].set_linewidth(3)
  894. ax.spines['bottom'].set_color('black')
  895. ax.spines['left'].set_color('black')
  897. plt.tight_layout()
  898. plt.savefig(out_path + 'feature_correlation_plot_mean_SD.png', dpi=300, bbox_inches='tight', transparent=True)
  899. plt.close()
  901. # %% [markdown]
  902. # # 统计检验
  904. # %%
  905. import numpy as np
  906. import matplotlib.pyplot as plt
  907. import os
  908. import pandas as pd
  909. import statsmodels.api as sm
  910. import statsmodels.formula.api as smf
  911. from scipy import stats
  913. # Set paths
  914. input_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/'
  915. output_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/results/'
  916. os.makedirs(output_path, exist_ok=True)
  918. # Load and prepare data
  919. corr_values_all = np.load(input_path + 'gradient_predict_task.npy') # shape: (91, 4, 3)
  920. corr_values_all = np.arctanh(corr_values_all)
  922. # Prepare DataFrame for mixed model
  923. data_list = []
  924. subjects = np.arange(corr_values_all.shape[0])
  925. blocks = np.arange(corr_values_all.shape[1])
  927. for subj in subjects:
  928. for block in blocks:
  929. group = 'Hearing' if subj < 34 else 'DNS'
  930. data_list.append({
  931. 'Subject': f'sub-{subj:03d}',
  932. 'Block': block,
  933. 'Group': group,
  934. 'Correlation': corr_values_all[subj, block, 0]
  935. })
  937. df = pd.DataFrame(data_list)
  939. # Fit mixed effects model
  940. model = smf.mixedlm("Correlation ~ Group",
  941. data=df,
  942. groups=df["Subject"],
  943. re_formula="~1").fit()
  945. # Get p-value for group effect
  946. p_value = model.pvalues['Group[T.Hearing]']
  948. # Prepare plotting data
  949. hearing = corr_values_all[:34,:,0].flatten()
  950. deaf = corr_values_all[34:,:,0].flatten()
  952. # Create figure
  953. plt.figure(figsize=(8, 6))
  954. plt.rcParams.update({
  955. 'font.size': 12,
  956. 'axes.linewidth': 1.5,
  957. 'xtick.major.width': 1.5,
  958. 'ytick.major.width': 1.5
  959. })
  961. # Plot boxplots with thicker borders (linewidth=3)
  962. boxprops = dict(linewidth=3, color='black') # Thicker borders
  963. whiskerprops = dict(linewidth=3, color='black') # Thicker whiskers
  964. capprops = dict(linewidth=3, color='black') # Thicker caps
  965. medianprops = dict(linewidth=3, color='black') # Thicker median line
  967. bp = plt.boxplot([hearing, deaf],
  968. labels=['Hearing', 'DNS'],
  969. patch_artist=True,
  970. widths=0.6,
  971. boxprops=boxprops,
  972. whiskerprops=whiskerprops,
  973. capprops=capprops,
  974. medianprops=medianprops)
  976. # Set box fill colors to #F08E64 (orange) and #5E93CB (blue)
  977. colors = ['#F08E64', '#5E93CB']
  978. for patch, color in zip(bp['boxes'], colors):
  979. patch.set_facecolor(color)
  980. patch.set_alpha(0.7)
  982. # Plot individual data points with jitter (larger size)
  983. def plot_jitter(data, x_pos, color, alpha=0.6):
  984. jitter = np.random.normal(0, 0.05, size=len(data))
  985. plt.scatter(x_pos + jitter, data, color=color, alpha=alpha, s=250, edgecolor='none')
  987. plot_jitter(hearing, 1, '#F08E64', alpha=0.6)
  988. plot_jitter(deaf, 2, '#5E93CB', alpha=0.6)
  990. # Calculate standard errors for error bars
  991. se_hearing = stats.sem(hearing)
  992. se_deaf = stats.sem(deaf)
  994. # Mark means with white diamonds and thicker error bars
  995. mean_hearing = np.mean(hearing)
  996. mean_deaf = np.mean(deaf)
  998. # Plot mean markers with thicker error bars (linewidth=3)
  999. plt.errorbar([1], [mean_hearing], #yerr=se_hearing,
  1000. fmt='none', ecolor='black', elinewidth=3, capsize=10, capthick=3)
  1001. plt.errorbar([2], [mean_deaf], #yerr=se_deaf,
  1002. fmt='none', ecolor='black', elinewidth=3, capsize=10, capthick=3)
  1004. plt.scatter(1, mean_hearing, color='white', s=250, edgecolor='black',
  1005. marker='D', zorder=10, linewidth=2, label='Mean')
  1006. plt.scatter(2, mean_deaf, color='white', s=250, edgecolor='black',
  1007. marker='D', zorder=10, linewidth=2)
  1009. # Add significance annotation
  1010. y_max = max(np.max(hearing), np.max(deaf)) + 0.1
  1011. plt.plot([1, 1, 2, 2], [y_max-0.05, y_max, y_max, y_max-0.05],
  1012. color='black', lw=3) # Thicker significance line
  1014. if p_value < 0.001:
  1015. sig_text = '***'
  1016. elif p_value < 0.01:
  1017. sig_text = '**'
  1018. elif p_value < 0.05:
  1019. sig_text = '*'
  1020. else:
  1021. sig_text = 'n.s.'
  1023. plt.text(1.5, y_max, sig_text,
  1024. ha='center', va='bottom',
  1025. fontsize=48, fontweight='bold')
  1027. # Customize plot
  1028. plt.xticks(fontsize=46, fontweight='bold')
  1029. plt.yticks(fontsize=46, fontweight='bold')
  1031. # Adjust spines (thicker)
  1032. ax = plt.gca()
  1033. ax.spines['top'].set_visible(False)
  1034. ax.spines['right'].set_visible(False)
  1035. ax.spines['bottom'].set_linewidth(4)
  1036. ax.spines['left'].set_linewidth(4)
  1037. ax.spines['bottom'].set_color('black')
  1038. ax.spines['left'].set_color('black')
  1040. plt.tight_layout()
  1041. plt.savefig(output_path + 'hearing_vs_deaf_boxplot_thick_borders.png',
  1042. dpi=300, bbox_inches='tight', transparent=True)
  1043. plt.show()
  1045. # Print model summary
  1046. print(model.summary())
  1048. # %%
  1049. import numpy as np
  1050. import os
  1051. import matplotlib.pyplot as plt
  1052. import seaborn as sns
  1053. import pandas as pd
  1054. import statsmodels.api as sm
  1055. import statsmodels.formula.api as smf
  1056. from scipy.stats import zscore
  1057. from matplotlib.colors import LinearSegmentedColormap
  1058. from statsmodels.stats.multitest import multipletests # For p-value correction
  1060. # Set paths
  1061. input_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/'
  1062. output_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/results/'
  1063. os.makedirs(output_path, exist_ok=True)
  1065. # Define gradient color schemes
  1066. hearing_colors = LinearSegmentedColormap.from_list('hearing', ['#F08E64', '#F8C8B0'])(np.linspace(0, 1, 10))
  1067. deaf_colors = LinearSegmentedColormap.from_list('DNS', ['#5E93CB', '#B8D1E8'])(np.linspace(0, 1, 10))
  1069. # Load and standardize data
  1070. shap_values_all = np.load(input_path + 'shap_feature_importance.npy') # shape: (91, 4, 10) or (91, 8, 10)
  1071. shap_values_all = zscore(shap_values_all, axis=2) # Standardize within each dimension
  1073. # Define group sizes
  1074. n_Hearing = 34
  1075. n_Deaf = 27
  1076. n_subjects = n_Hearing + n_Deaf
  1077. n_blocks = 8
  1078. n_dims = 10
  1080. # Prepare DataFrame for mixed model
  1081. data_list = []
  1082. for subj in range(n_subjects):
  1083. group = 'Hearing' if subj < n_Hearing else 'DNS'
  1084. for block in range(n_blocks):
  1085. for dim in range(n_dims):
  1086. data_list.append({
  1087. 'subject': f'subj_{subj}',
  1088. 'group': group,
  1089. 'block': f'block_{block}',
  1090. 'dimension': dim,
  1091. 'shap_value': shap_values_all[subj, block, dim]
  1092. })
  1094. df = pd.DataFrame(data_list)
  1096. # Initialize storage for results
  1097. lmm_results = []
  1099. # Fit mixed models for each SHAP dimension
  1100. for dim in range(n_dims):
  1101. dim_df = df[df['dimension'] == dim].copy()
  1103. # Linear Mixed Effects Model with subject random intercept
  1104. model = smf.mixedlm(
  1105. "shap_value ~ group",
  1106. data=dim_df,
  1107. groups=dim_df["subject"],
  1108. re_formula="~1" # Random intercept per subject
  1109. )
  1111. try:
  1112. result = model.fit()
  1113. coef = result.fe_params['group[T.Hearing]'] # Effect of Hearing vs Deaf
  1114. pval = result.pvalues['group[T.Hearing]']
  1115. t_stat = result.tvalues['group[T.Hearing]']
  1117. lmm_results.append({
  1118. 'dimension': dim,
  1119. 'coefficient': coef,
  1120. 't_statistic': t_stat,
  1121. 'p_value': pval,
  1122. 'significant': pval < 0.05 # Direct significance without correction
  1123. })
  1124. except Exception as e:
  1125. print(f"Failed to fit model for dimension {dim}: {str(e)}")
  1126. lmm_results.append({
  1127. 'dimension': dim,
  1128. 'coefficient': np.nan,
  1129. 't_statistic': np.nan,
  1130. 'p_value': np.nan,
  1131. 'significant': False
  1132. })
  1134. # Convert to DataFrame
  1135. results_df = pd.DataFrame(lmm_results)
  1137. # Apply p-value correction (FDR)
  1138. pvals = results_df['p_value'].values
  1139. rejected, pvals_corrected, _, _ = multipletests(pvals, alpha=0.05, method='fdr_bh')
  1140. results_df['p_value_corrected'] = pvals_corrected
  1141. results_df['significant_corrected'] = rejected # Update significance based on corrected p-values
  1143. # Save corrected results
  1144. results_df.to_csv(output_path + 'lmm_results_corrected.csv', index=False)
  1146. # Visualization with more pronounced gradients
  1147. plt.figure(figsize=(12, 6.5)) # Larger figure size
  1149. # Calculate group means (averaged across blocks)
  1150. group_means = df.groupby(['group', 'dimension'])['shap_value'].mean().unstack()
  1151. group_sems = df.groupby(['group', 'dimension'])['shap_value'].sem().unstack()
  1153. # ===================== 新增：计算“每个被试在每个维度”的单个点（跨block平均） =====================
  1154. subj_means = df.groupby(['group', 'subject', 'dimension'])['shap_value'].mean()
  1155. hearing_subjects = subj_means.loc['Hearing'].index.get_level_values('subject').unique()
  1156. dns_subjects = subj_means.loc['DNS'].index.get_level_values('subject').unique()
  1158. np.random.seed(0) # 抖动可复现
  1159. # =============================================================================================
  1161. # Plot with more pronounced gradient colors
  1162. x = range(n_dims)
  1163. for dim in range(n_dims):
  1165. # ===================== 新增：叠加被试散点（颜色与对应维度一致，透明，不遮挡误差线） =====================
  1166. # Hearing 所有被试点（向左轻微抖动）
  1167. y_hearing_pts = np.array([subj_means.loc[('Hearing', s, dim)] for s in hearing_subjects])
  1168. x_hearing_pts = dim + np.random.uniform(-0.18, -0.06, size=len(y_hearing_pts))
  1169. plt.scatter(x_hearing_pts, y_hearing_pts,
  1170. s=75, color=hearing_colors[dim], alpha=0.88,
  1171. edgecolors='none', zorder=1)
  1173. # DNS 所有被试点（向右轻微抖动）
  1174. y_dns_pts = np.array([subj_means.loc[('DNS', s, dim)] for s in dns_subjects])
  1175. x_dns_pts = dim + np.random.uniform(0.06, 0.18, size=len(y_dns_pts))
  1176. plt.scatter(x_dns_pts, y_dns_pts,
  1177. s=75, color=deaf_colors[dim], alpha=0.88,
  1178. edgecolors='none', zorder=1)
  1179. # ==========================================================================================
  1181. # Hearing markers - darker to lighter
  1182. plt.errorbar(dim, group_means.loc['Hearing', dim],
  1183. yerr=group_sems.loc['Hearing', dim],
  1184. fmt='o', color=hearing_colors[dim],
  1185. markersize=20, capsize=6, capthick=3,
  1186. elinewidth=3, markeredgewidth=2,
  1187. markeredgecolor='black', zorder=7)
  1189. # Deaf markers - darker to lighter
  1190. plt.errorbar(dim, group_means.loc['DNS', dim],
  1191. yerr=group_sems.loc['DNS', dim],
  1192. fmt='s', color=deaf_colors[dim],
  1193. markersize=20, capsize=6, capthick=3,
  1194. elinewidth=3, markeredgewidth=2,
  1195. markeredgecolor='black', zorder=7)
  1197. # Connect the dots with gradient-colored lines
  1198. for dim in range(n_dims-1):
  1199. plt.plot([dim, dim+1],
  1200. [group_means.loc['Hearing', dim], group_means.loc['Hearing', dim+1]],
  1201. '-', color=hearing_colors[dim], linewidth=4, alpha=0.8)
  1202. plt.plot([dim, dim+1],
  1203. [group_means.loc['DNS', dim], group_means.loc['DNS', dim+1]],
  1204. '-', color=deaf_colors[dim], linewidth=4, alpha=0.8)
  1206. # Mark significant dimensions with larger stars (based on corrected p-values)
  1207. sig_dims_corrected = results_df[results_df['significant_corrected']]['dimension']
  1208. for dim in sig_dims_corrected:
  1209. y_pos = max(group_means.loc['Hearing', dim], group_means.loc['DNS', dim]) + 0.15
  1210. plt.scatter(dim, y_pos, marker='*', color='black', s=600, zorder=10)
  1212. # Customize plot appearance
  1213. plt.ylabel('Standardized SHAP Value', fontsize=29, fontweight='bold', labelpad=15)
  1214. plt.xticks(range(10), [f'G{i+1}' for i in range(10)], fontsize=30, fontweight='bold')
  1215. plt.yticks(fontsize=30, fontweight='bold')
  1217. # Create custom legend
  1218. from matplotlib.lines import Line2D
  1219. legend_elements = [
  1220. Line2D([0], [0], marker='o', color='w', label='Hearing',
  1221. markerfacecolor=hearing_colors[0], markersize=25, markeredgewidth=2),
  1222. Line2D([0], [0], marker='s', color='w', label='DNS',
  1223. markerfacecolor=deaf_colors[0], markersize=25, markeredgewidth=2)
  1224. ]
  1225. plt.legend(handles=legend_elements, fontsize=48, prop={'weight':'bold', 'size':26},
  1226. framealpha=0.9, loc='upper right')
  1228. # Customize spines
  1229. ax = plt.gca()
  1230. ax.spines['top'].set_visible(False)
  1231. ax.spines['right'].set_visible(False)
  1232. ax.spines['bottom'].set_linewidth(4)
  1233. ax.spines['left'].set_linewidth(4)
  1234. ax.spines['bottom'].set_color('black')
  1235. ax.spines['left'].set_color('black')
  1236. plt.grid(False)
  1237. plt.tight_layout()
  1239. plt.savefig(output_path + 'shap_comparison_gradient_colors_corrected.png',
  1240. dpi=300, bbox_inches='tight', transparent=True)
  1241. plt.close()
  1243. print("Analysis complete. Results saved to:", output_path)
  1245. # %% [markdown]
  1246. # # 对比有无语言经验的差异
  1248. # %%
  1249. import numpy as np
  1250. import matplotlib.pyplot as plt
  1251. import os
  1252. import pandas as pd
  1253. import statsmodels.api as sm
  1254. import statsmodels.formula.api as smf
  1255. from scipy import stats
  1257. # Set paths
  1258. input_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/'
  1259. output_path = 'J:/HCP_analysis/Deaf/10_rest_predict_task/results_lan_ELA/'
  1260. os.makedirs(output_path, exist_ok=True)
  1262. # Load and prepare data
  1263. corr_values_all = np.load(input_path + 'gradient_predict_task.npy') # shape: (91, 4, 3)
  1264. corr_values_all = np.arctanh(corr_values_all)
  1266. # Prepare DataFrame for mixed model
  1267. data_list = []
  1268. subjects = np.arange(34, corr_values_all.shape[0]) # Only deaf subjects (34-90)
  1269. blocks = np.arange(corr_values_all.shape[1])
  1271. for i, subj in enumerate(subjects):
  1272. for block in blocks:
  1273. group = 'ELA' if i < 27 else 'Non_ELA' # ELA:34-60 (27), Non_ELA:61-90 (30)
  1274. data_list.append({
  1275. 'Subject': f'sub-{subj:03d}',
  1276. 'Block': block,
  1277. 'Group': group,
  1278. 'Correlation': corr_values_all[subj, block, 0]
  1279. })
  1281. df = pd.DataFrame(data_list)
  1283. # Fit mixed effects model
  1284. model = smf.mixedlm("Correlation ~ Group",
  1285. data=df,
  1286. groups=df["Subject"],
  1287. re_formula="~1").fit()
  1289. # Get p-value for group effect
  1290. p_value = model.pvalues['Group[T.Non_ELA]'] # Reference is ELA
  1292. # Prepare plotting data
  1293. ELA = corr_values_all[34:61,:,0].flatten()
  1294. Non_ELA = corr_values_all[61:,:,0].flatten()
  1296. # Create figure with larger size
  1297. plt.figure(figsize=(8, 6))
  1298. plt.rcParams.update({
  1299. 'font.size': 12,
  1300. 'axes.linewidth': 2,
  1301. 'xtick.major.width': 2,
  1302. 'ytick.major.width': 2
  1303. })
  1305. # Plot boxplots with thicker elements (linewidth=3)
  1306. boxprops = dict(linewidth=3, color='black') # Thicker box borders
  1307. whiskerprops = dict(linewidth=3, color='black') # Thicker whiskers
  1308. capprops = dict(linewidth=3, color='black') # Thicker caps
  1309. medianprops = dict(linewidth=3, color='black') # Thicker median line
  1311. bp = plt.boxplot([ELA, Non_ELA],
  1312. labels=['DNS', 'DLS'],
  1313. patch_artist=True,
  1314. widths=0.6,
  1315. boxprops=boxprops,
  1316. whiskerprops=whiskerprops,
  1317. capprops=capprops,
  1318. medianprops=medianprops)
  1320. # Set box fill colors with new color scheme
  1321. colors = ['#5E93CB', '#9A4942'] # DNS and DLS colors '#5E93CB' '#79902D'
  1322. for patch, color in zip(bp['boxes'], colors):
  1323. patch.set_facecolor(color)
  1324. patch.set_alpha(0.7)
  1326. # Plot individual data points with jitter (larger size)
  1327. def plot_jitter(data, x_pos, color, alpha=0.6):
  1328. jitter = np.random.normal(0, 0.05, size=len(data))
  1329. plt.scatter(x_pos + jitter, data, color=color, alpha=alpha, s=250, edgecolor='none')
  1331. plot_jitter(ELA, 1, '#5E93CB')
  1332. plot_jitter(Non_ELA, 2, '#9A4942')
  1334. # Calculate standard errors for error bars
  1335. se_ELA = stats.sem(ELA)
  1336. se_Non_ELA = stats.sem(Non_ELA)
  1338. # Mark means with larger diamond shapes (s=150) and thicker error bars
  1339. mean_ELA = np.mean(ELA)
  1340. mean_Non_ELA = np.mean(Non_ELA)
  1342. # Plot error bars with thicker lines (linewidth=3)
  1343. plt.errorbar([1], [mean_ELA], #yerr=se_ELA,
  1344. fmt='none', ecolor='black', elinewidth=3, capsize=12, capthick=3)
  1345. plt.errorbar([2], [mean_Non_ELA], #yerr=se_Non_ELA,
  1346. fmt='none', ecolor='black', elinewidth=3, capsize=12, capthick=3)
  1348. # 显著增大均值标记的大小到s=300
  1349. plt.scatter(1, mean_ELA, color='white', s=250, edgecolor='black', # 增大到300
  1350. marker='D', zorder=10, linewidth=3, label='Mean') # 同时加粗边框到3
  1351. plt.scatter(2, mean_Non_ELA, color='white', s=250, edgecolor='black',
  1352. marker='D', zorder=10, linewidth=3)
  1354. # Add significance annotation with thicker line
  1355. y_max = max(np.max(ELA), np.max(Non_ELA)) + 0.1
  1356. plt.plot([1, 1, 2, 2], [y_max-0.05, y_max, y_max, y_max-0.05],
  1357. color='black', lw=3) # Thicker significance line
  1359. if p_value < 0.001:
  1360. sig_text = '***'
  1361. elif p_value < 0.01:
  1362. sig_text = '**'
  1363. elif p_value < 0.05:
  1364. sig_text = '*'
  1365. else:
  1366. sig_text = 'n.s.'
  1368. plt.text(1.5, y_max, sig_text,
  1369. ha='center', va='bottom',
  1370. fontsize=48, fontweight='bold')
  1372. # Customize plot with larger fonts
  1373. plt.xticks(fontsize=48, fontweight='bold')
  1374. plt.yticks(fontsize=48, fontweight='bold')
  1376. # Adjust spines (thicker)
  1377. ax = plt.gca()
  1378. ax.spines['top'].set_visible(False)
  1379. ax.spines['right'].set_visible(False)
  1380. ax.spines['bottom'].set_linewidth(3)
  1381. ax.spines['left'].set_linewidth(3)
  1382. ax.spines['bottom'].set_color('black')
  1383. ax.spines['left'].set_color('black')
  1385. plt.tight_layout()
  1386. plt.savefig(output_path + 'ELA_vs_Non_ELA_boxplot_enhanced.png',
  1387. dpi=300, bbox_inches='tight', transparent=True)
  1388. plt.show()
  1390. # Print model summary
  1391. print("Mixed Effects Model Results:")
  1392. print(model.summary())
  1393. print(f"\nGroup difference p-value: {p_value:.4f}")
  1395. # %%
  1396. import numpy as np
  1397. import os
  1398. import matplotlib.pyplot as plt
  1399. import seaborn as sns
  1400. import pandas as pd
  1401. import statsmodels.api as sm
  1402. import statsmodels.formula.api as smf
  1403. from scipy.stats import zscore
  1404. from statsmodels.stats.multitest import multipletests # For p-value correction
  1406. # Set paths
  1407. input_path = 'J:/HCP_analysis/deaf/10_rest_predict_task/'
  1408. output_path = 'J:/HCP_analysis/deaf/10_rest_predict_task/results_lan_DNS/'
  1409. os.makedirs(output_path, exist_ok=True)
  1411. # Load and standardize data
  1412. shap_values_all = np.load(input_path + 'shap_feature_importance.npy') # shape: (91, 4, 10)
  1413. shap_values_all = zscore(shap_values_all, axis=2) # Standardize within each dimension
  1415. # Define group sizes
  1416. n_DNS = 27
  1417. n_Non_DNS = 30
  1418. n_subjects = n_DNS + n_Non_DNS
  1419. n_blocks = 4
  1420. n_dims = 10
  1422. # Prepare DataFrame for mixed model
  1423. data_list = []
  1424. for subj in range(n_subjects):
  1425. group = 'DNS' if subj < n_DNS else 'DLS'
  1426. actual_subj_idx = 34 + subj # Adjust index since DNS group starts at position 34
  1427. for block in range(n_blocks):
  1428. for dim in range(n_dims):
  1429. data_list.append({
  1430. 'subject': f'subj_{actual_subj_idx}',
  1431. 'group': group,
  1432. 'block': block, # Changed to numeric
  1433. 'dimension': dim,
  1434. 'shap_value': shap_values_all[actual_subj_idx, block, dim]
  1435. })
  1437. df = pd.DataFrame(data_list)
  1439. # Ensure group is categorical and set reference level
  1440. df['group'] = pd.Categorical(df['group'], categories=['DLS', 'DNS'])
  1441. df['block'] = pd.Categorical(df['block']) # Make block categorical
  1443. # Initialize storage for results
  1444. lmm_results = []
  1446. # Fit mixed models for each SHAP dimension
  1447. for dim in range(n_dims):
  1448. dim_df = df[df['dimension'] == dim].copy()
  1450. try:
  1451. model = smf.mixedlm(
  1452. "shap_value ~ C(group, Treatment('DLS'))",
  1453. data=dim_df,
  1454. groups=dim_df["subject"],
  1455. re_formula="~1"
  1456. )
  1458. result = model.fit(reml=True)
  1460. coef = result.fe_params["C(group, Treatment('DLS'))[T.DNS]"]
  1461. pval = result.pvalues["C(group, Treatment('DLS'))[T.DNS]"]
  1462. t_stat = result.tvalues["C(group, Treatment('DLS'))[T.DNS]"]
  1464. lmm_results.append({
  1465. 'dimension': dim,
  1466. 'coefficient': coef,
  1467. 't_statistic': t_stat,
  1468. 'p_value': pval,
  1469. 'significant': pval < 0.05,
  1470. 'random_effect_variance': result.cov_re.iloc[0, 0]
  1471. })
  1472. except Exception as e:
  1473. print(f"Failed to fit model for dimension {dim}: {str(e)}")
  1474. lmm_results.append({
  1475. 'dimension': dim,
  1476. 'coefficient': np.nan,
  1477. 't_statistic': np.nan,
  1478. 'p_value': np.nan,
  1479. 'significant': False,
  1480. 'random_effect_variance': np.nan
  1481. })
  1483. # Convert to DataFrame
  1484. results_df = pd.DataFrame(lmm_results)
  1486. # Apply p-value correction (FDR)
  1487. pvals = results_df['p_value'].values
  1488. rejected, pvals_corrected, _, _ = multipletests(pvals, alpha=0.05, method='fdr_bh')
  1489. results_df['p_value_corrected'] = pvals_corrected
  1490. results_df['significant_corrected'] = rejected
  1492. # Save corrected results
  1493. results_df.to_csv(output_path + 'lmm_results_DNS_vs_DLS_corrected.csv', index=False)
  1495. # Visualization with gradient colors and black borders
  1496. plt.figure(figsize=(12, 6.5))
  1498. # Calculate group means (averaged across blocks)
  1499. # 用“被试层面均值”（跨 4 个 block 平均）来计算组均值与SEM，避免SEM过小
  1500. df_subj = df.groupby(['group', 'subject', 'dimension'], as_index=False)['shap_value'].mean()
  1502. group_means = df_subj.groupby(['group', 'dimension'])['shap_value'].mean().unstack()
  1503. group_sems = df_subj.groupby(['group', 'dimension'])['shap_value'].sem().unstack()
  1505. # Create color gradients for each group
  1506. dns_base = np.array([0.3686, 0.5765, 0.7961]) # base blue
  1507. dls_base = np.array([0.6039, 0.2863, 0.2588]) # base red
  1509. dns_gradient = [dns_base * (1 - i*0.08) + np.array([1, 1, 1]) * (i*0.08) for i in range(10)]
  1510. dls_gradient = [dls_base * (1 - i*0.08) + np.array([1, 1, 1]) * (i*0.08) for i in range(10)]
  1512. # ===================== 新增：计算“每个被试在每个维度”的单个点（跨block平均） =====================
  1513. # 这就是你说的 shap_feature_importance.npy 里对应条件的“被试点”
  1514. subj_means = df.groupby(['group', 'subject', 'dimension'])['shap_value'].mean()
  1515. dns_subjects = subj_means.loc['DNS'].index.get_level_values('subject').unique()
  1516. dls_subjects = subj_means.loc['DLS'].index.get_level_values('subject').unique()
  1518. np.random.seed(0) # 仅用于抖动可复现（不想固定可删）
  1519. # =============================================================================================
  1521. # Plot each dimension with its own color from the gradient and black borders
  1522. for dim in range(10):
  1523. # ===================== 新增：叠加“所有被试点”（颜色与条件/维度一致，透明不遮挡误差线） =====================
  1524. # DNS 所有被试
  1525. y_dns_pts = np.array([subj_means.loc[('DNS', s, dim)] for s in dns_subjects])
  1526. x_dns_pts = dim + np.random.uniform(-0.18, -0.06, size=len(y_dns_pts))
  1527. plt.scatter(x_dns_pts, y_dns_pts,
  1528. s=85, color=dns_gradient[dim], alpha=0.88,
  1529. edgecolors='white', zorder=1)
  1531. # DLS 所有被试
  1532. y_dls_pts = np.array([subj_means.loc[('DLS', s, dim)] for s in dls_subjects])
  1533. x_dls_pts = dim + np.random.uniform(0.06, 0.18, size=len(y_dls_pts))
  1534. plt.scatter(x_dls_pts, y_dls_pts,
  1535. s=85, color=dls_gradient[dim], alpha=0.88,
  1536. edgecolors='white', zorder=1)
  1537. # ====================================================================================================
  1539. # DNS plot with black borders
  1540. plt.errorbar(dim, group_means.loc['DNS', dim],
  1541. yerr=group_sems.loc['DNS', dim],
  1542. fmt='o', zorder=5, capsize=6, linewidth=3, markersize=20,
  1543. color=dns_gradient[dim], ecolor='black',
  1544. elinewidth=2,
  1545. markeredgecolor='black', markeredgewidth=2,
  1546. label='DNS' if dim == 0 else "")
  1548. # DLS plot with black borders
  1549. plt.errorbar(dim, group_means.loc['DLS', dim],
  1550. yerr=group_sems.loc['DLS', dim],
  1551. fmt='s', zorder=5, capsize=6, linewidth=3, markersize=20,
  1552. color=dls_gradient[dim], ecolor='black',
  1553. elinewidth=2,
  1554. markeredgecolor='black', markeredgewidth=2,
  1555. label='DLS' if dim == 0 else "")
  1558. # Connect the dots with lines in the gradient colors
  1559. plt.plot(range(10), group_means.loc['DNS'],
  1560. color='#5E93CB', linewidth=3, alpha=1)
  1561. plt.plot(range(10), group_means.loc['DLS'],
  1562. color='#9A4942', linewidth=3, alpha=1)
  1564. # Mark significant differences (corrected p < 0.05)
  1565. for i in range(10):
  1566. if results_df['significant_corrected'].iloc[i]:
  1567. y_pos = max(group_means.loc['DNS', i], group_means.loc['DLS', i]) + 0.03
  1568. plt.scatter(i, y_pos + 0.2, marker='*', color='black', s=600, zorder=10)
  1570. plt.ylabel('Standardized SHAP Value', fontsize=29, fontweight='bold', labelpad=10)
  1571. plt.xticks(range(10), [f'G{i+1}' for i in range(10)], fontsize=30, fontweight='bold')
  1572. plt.yticks(fontsize=30, fontweight='bold')
  1574. # Create custom legend with black borders
  1575. from matplotlib.lines import Line2D
  1576. legend_elements = [
  1577. Line2D([0], [0], marker='o', color='w', label='DNS',
  1578. markerfacecolor='#5E93CB', markersize=25,
  1579. markeredgecolor='black', markeredgewidth=2),
  1580. Line2D([0], [0], marker='s', color='w', label='DLS',
  1581. markerfacecolor='#9A4942', markersize=25,
  1582. markeredgecolor='black', markeredgewidth=2)
  1583. ]
  1584. plt.legend(handles=legend_elements, fontsize=26,
  1585. prop={'weight': 'bold', 'size': 26},
  1586. framealpha=1, loc='upper right')
  1588. # Customize spines
  1589. ax = plt.gca()
  1590. ax.spines['top'].set_visible(False)
  1591. ax.spines['right'].set_visible(False)
  1592. ax.spines['bottom'].set_linewidth(4)
  1593. ax.spines['left'].set_linewidth(4)
  1594. ax.spines['bottom'].set_color('black')
  1595. ax.spines['left'].set_color('black')
  1597. plt.tight_layout()
  1598. plt.savefig(output_path + 'shap_comparison_gradient_DNS_DLS_bordered_corrected.png',
  1599. dpi=300, bbox_inches='tight', transparent=True)
  1600. plt.close()
  1602. print("Analysis complete. Results saved to:", output_path)
  1604. # %%

S05_the_rsFC_gradient_rest_predict_task_activation.ipynb, no license · at the source

Overview

Authors: Linyan Liu1,2, Chan Tang1,2, Jie Chen1,2, Yang Lei1,2, Yanyan Li3,4, Guosheng Ding1,2
ORCID iDs: Linyan Liu, Chan Tang, Yanyan Li, Guosheng Ding
  1. Faculty of Psychology, Beijing Normal University,Beijing, China
  2. State Key Laboratory of Cognitive Neuroscience and Learning & IDG/McGovern Institute for Brain Research, Beijing Normal University,Beijing, China
  3. State Key Laboratory of Cognitive Science and Mental Health, Institute of Psychology, Chinese Academy of Sciences,Beijing, China
  4. Department of Psychology, University of Chinese Academy of Sciences,Beijing, China
Institutions: Beijing Normal University (China); Chinese Academy of Sciences (China); University of Chinese Academy of Sciences (China)
Journal: Communications biology, volume 9, issue 1, article 703
Dates: received 1 September 2025; accepted 9 March 2026; published online 30 March 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s42003-026-09900-8 · PMID 41912615 · PMCID PMC13199478 · OpenAlex W7143385668
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), other condition (population), cognitive (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging
Keywords: Language, Sensory processing
MeSH: Auditory Perception*, Brain*, Deafness*, Language*, Neuronal Plasticity*, Sensory Deprivation*, Acoustic Stimulation, Auditory Cortex, Brain Mapping, Female, Humans, Magnetic Resonance Imaging, Male (* major topic)
Topic: Language Development and Disorders (Developmental and Educational Psychology, Psychology), according to OpenAlex
Funding: National Natural Science Foundation of China (32471100, 32371141)
Citations: cited by 1 paper (Europe PMC); 61 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 10 matches between paragraphs and lines of code.

OSF ebsvz

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: Jupyter (6)
Size: 31 files, 6 scripts
Software Heritage: not checked
Found in: “Code availability”
Holds: 6 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: BrainSpace (6 files), h5py (6 files), Matplotlib (6 files), neuromaps (6 files), NiBabel (6 files), Nilearn (6 files), NumPy (6 files), pandas (6 files), scikit-learn (6 files), SciPy (6 files), seaborn (6 files), statsmodels (6 files), SHAP (4 files), Pingouin (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
6 files
At the source: osf.io/ebsvz/

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s42003-026-09900-8.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 6 scripts, each with its path and the digest of its content;
  • 10 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s42003-026-09900-8.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 2 keywords, 13 MeSH terms, 1 funder, 61 references.

Cite

This paper

Liu, L., Tang, C., Chen, J., Lei, Y., Li, Y., & Ding, G. (2026). A multidimensional framework for dissociating the neuroplasticity of auditory and early language deprivation. Communications biology, 9(1), 703. https://doi.org/10.1038/s42003-026-09900-8

BibTeX

@article{liu2026multidimensional,
author = {Liu, Linyan and Tang, Chan and Chen, Jie and Lei, Yang and Li, Yanyan and Ding, Guosheng},
title = {{A multidimensional framework for dissociating the neuroplasticity of auditory and early language deprivation}},
journal = {Communications biology},
year = {2026},
month = mar,
volume = {9},
number = {1},
pages = {703},
publisher = {Nature Publishing Group},
issn = {2399-3642},
doi = {10.1038/s42003-026-09900-8},
url = {https://doi.org/10.1038/s42003-026-09900-8},
pmid = {41912615},
pmcid = {PMC13199478}
}

RIS

TY - JOUR
AU - Liu, Linyan
AU - Tang, Chan
AU - Chen, Jie
AU - Lei, Yang
AU - Li, Yanyan
AU - Ding, Guosheng
TI - A multidimensional framework for dissociating the neuroplasticity of auditory and early language deprivation
T2 - Communications biology
J2 - Commun Biol
PY - 2026
DA - 2026/03/30
VL - 9
IS - 1
SP - 703
SN - 2399-3642
PB - Nature Publishing Group
DO - 10.1038/s42003-026-09900-8
UR - https://doi.org/10.1038/s42003-026-09900-8
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42003-026-09900-8",
"type": "article-journal",
"title": "A multidimensional framework for dissociating the neuroplasticity of auditory and early language deprivation",
"container-title": "Communications biology",
"author": [
{
"family": "Liu",
"given": "Linyan"
},
{
"family": "Tang",
"given": "Chan"
},
{
"family": "Chen",
"given": "Jie"
},
{
"family": "Lei",
"given": "Yang"
},
{
"family": "Li",
"given": "Yanyan"
},
{
"family": "Ding",
"given": "Guosheng"
}
],
"container-title-short": "Commun Biol",
"volume": "9",
"issue": "1",
"page": "703",
"DOI": "10.1038/s42003-026-09900-8",
"PMID": "41912615",
"PMCID": "PMC13199478",
"ISSN": "2399-3642",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s42003-026-09900-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
30
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.21203/rs.3.rs-9326213/v1 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: neuromaps, Nilearn, h5py, 7 other tools, cognitive, 7 references
[2] doi:10.64898/2026.03.09.710558 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: neuromaps, Nilearn, h5py, 7 other tools, cognitive, 7 references
[3] doi:10.7554/elife.96944 [code]
Unraveling the impact of congenital deafness on individual brain organization.
Journal: eLife
In common: other condition, 11 references
[4] doi:10.1162/netn_a_00382 [code]
Comparison of different group-level templates in gradient-based multimodal connectivity analysis.
Journal: Network neuroscience (Cambridge, Mass.)
In common: BrainSpace, Nilearn, h5py, 8 other tools, 6 references
[5] doi:10.1162/netn_a_00402 [code]
The impact of functional correlations on task information coding.
Journal: Network neuroscience (Cambridge, Mass.)
In common: neuromaps, Nilearn, h5py, 8 other tools, 6 references
[6] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: neuromaps, BrainSpace, Pingouin, 8 other tools, 4 references
[7] doi:10.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: neuromaps, BrainSpace, Nilearn, 9 other tools, 4 references
[8] doi:10.1162/imag_a_00474 [code]
Exploring neuroendocrine influences on the sensorimotor-association axis in a female and a male individual.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: neuromaps, BrainSpace, Nilearn, 8 other tools, 5 references
[9] doi:10.1038/s41467-025-66783-9 [code]
Evidence for divergent cortical organisation in Parkinson's disease and Lewy Body Dementia.
Journal: Nature communications
In common: neuromaps, BrainSpace, Pingouin, 8 other tools, 4 references
[10] doi:10.1038/s41467-024-51942-1 [code]
Sex differences in functional cortical organization reflect differences in network topology rather than cortical morphometry.
Journal: Nature communications
In common: neuromaps, BrainSpace, Pingouin, 9 other tools, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8.svg)](https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8.svg
   :target: https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8"><img src="https://openscicode.org/badge/doi/10.1038/s42003-026-09900-8.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.