OSCR

In-scanner thoughts contribute to resting-state functional connectivity.

Code ↔ Paper

15 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 15 matches · 3 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › In-scanner experience reports analyses › Clustering ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 415–442 · score 0.94 · Gaussian Mixture Modeling, Python library scikit, soft clustering nature, spherical clusters, provides membership probabilities, neuroimaging literature
  2. [2] § Methods › In-scanner experience reports analyses › Dimensionality reduction: T-SNE ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 277–310 · score 0.91 · draw biplot arrows, NSE embedding space, SNE dimensions relate, linear function, maximal change, SNE coordinates
  3. [3] § Methods › In-scanner experience reports analyses › Dimensionality reduction: T-SNE ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 218–253 · score 0.85 · explored hyper parameter, embedding preserves local, SNE algorithm, t-SNE, trustworthiness, perplexity
  4. [4] § Methods › fMRI Data Pre-processing ↔ bash/S06_NuissanceRegression.sh, lines 86–173 · score 0.79 · nuisance regression, noise regressors, 0.15 Hz, derivative, filtering, AFNI
  5. [5] § Methods › In-scanner experience reports analyses › Data scaling ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 66–75 · score 0.79 · RobustScaler, inter quantile range, excessive influence, scikit learn, median, outliers
  6. [6] § Methods › In-scanner experience reports analyses › Clustering ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 415–442 · score 0.79 · low dimensional space, clustering technique, clustering algorithms, inner experience, separate scans, robustness
  7. [7] § Methods › In-scanner experience reports analyses › Outlier detection ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 85–101 · score 0.65 · squared Mahalanobis distance, Outlier detection, quantile, scan
  8. [8] § Methods › fMRI Data Pre-processing ↔ bash/S04_TransformToMNI.pass01.sh, the whole file · a weak match · score 0.63 · MNI space, motion correction, Pre processing, AFNI, scans
  9. [9] § Methods › Literature search › Connectome-based predictive modeling studies ↔ notebooks/S16_CPM.Development.ipynb, lines 257–315 · score 0.60 · feature selection, cross validation, Fold, CPM, predicted, models
  10. [10] § Methods › Functional connectivity analyses - statistical differences across scan sets › Parcellation ↔ bash/S08_ExtractROIts.sh, the whole file · a weak match · score 0.57 · dNetCorr, pre processing, FC matrices, ROIs, AFNI, atlas
  11. [11] § Methods › Functional connectivity analyses – connectome-based predictive modeling ↔ notebooks/S16_CPM.Development.ipynb, lines 257–315 · score 0.56 · cross validation, gather, training, fold, CPM, metric
  12. [12] § Methods › fMRI Data Pre-processing ↔ bash/S04_TransformToMNI.pass01.sh, the whole file · a weak match · score 0.55 · ANTs, pre processing, MNI152, scan
  13. [13] § Methods › In-scanner experience reports analyses › Clustering ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 479–515 · score 0.53 · Adjusted Rand, subsamples, ARI, bootstrapping, clusters, SI
  14. [14] § Results › Self-report annotations of resting-state scans ↔ notebooks/S10_SNYCQ_TSNE_Clustering.ipynb, lines 479–515 · score 0.52 · adjusted Rand, subsampling, ARI, bootstrapping, fit, GMM
  15. [15] § Methods › Functional connectivity analyses - statistical differences across scan sets › Parcellation ↔ notebooks/S07_PrepareAtlas.ipynb, lines 1–32 · score 0.52 · Limbic network, subcortical regions, pallidum, AAL, putamen, ROIs

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,101 lines · 40 KB · CC0-1.0 · 8 matches

  1. # %% [markdown]
  2. # # Description: SNYCQ Initial Exploration, Clustering and TSNE
  3. #
  4. # This notebook contains the following analytical steps associated with the in-scanner experience data
  5. #
  6. # 1. Data scaling: this is accomplished using skicit-learn [```RobustScaler```](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.RobustScaler.html)
  7. #
  8. # 2. Outlier detection
  9. #
  10. # 3. Creates figure looking a potential correlations between sNYCQ items
  11. #
  12. # 4. Dimensionality reduction with T-SNE
  13. #
  14. # 5. Clustering analysis in original 11D space
  16. # %%
  17. from utils.basics import get_sbj_scan_list, RESOURCES_DINFO_DIR, RESOURCES_SNYCQ_DIR, ORIG_DEMO_PATH
  18. from utils.plotting import show_correlations_with_statistics
  19. import os.path as osp
  20. from scipy.stats import ttest_ind, mannwhitneyu, wilcoxon
  21. import hvplot.pandas
  22. import holoviews as hv
  23. import seaborn as sns
  24. from tqdm.notebook import tqdm
  26. import warnings
  27. warnings.simplefilter(action='ignore', category=FutureWarning)
  29. # %%
  30. import numpy as np
  31. import pandas as pd
  32. import matplotlib.pyplot as plt
  33. import plotly.graph_objects as go
  34. import panel as pn
  35. from scipy.stats import pearsonr
  37. from sklearn.preprocessing import RobustScaler
  38. from sklearn.covariance import MinCovDet
  39. from sklearn.mixture import GaussianMixture
  40. from sklearn.cluster import KMeans
  41. from sklearn.metrics import silhouette_score, adjusted_rand_score
  42. from sklearn.manifold import TSNE, trustworthiness
  43. from sklearn.linear_model import LinearRegression
  45. # %% [markdown]
  46. # Configurations for outlier detection, ambigous clusters and random seed
  48. # %%
  49. OUTLIER_Q = 0.997 # robust cutoff on squared Mahalanobis distances
  50. AMBIG_THRESH = 0.8 # ambiguous if max(prob) < 0.8
  51. RANDOM_STATE = 42
  53. # %% [markdown]
  54. # ***
  55. # # 1. Load In-scanner Experience Data (SNYCQ)
  56. #
  57. # We load this data only for the scans that have passed our QA for the imaging data
  59. # %%
  60. SBJs, SCANs, SNYCQ_wVigilance = get_sbj_scan_list(when='post_motion', return_snycq=True)
  61. SNYCQ = SNYCQ_wVigilance.drop('Vigilance',axis=1)
  62. Nscans, Nquestions = SNYCQ.shape
  63. print(SNYCQ.shape)
  64. SNYCQ_items = SNYCQ.columns
  66. # %% [markdown]
  67. # ***
  68. # # 2. Data Scaling
  69. # Answers to all questions, except wakefulness, were scaled using scikit-learn ```RobustScaler```. This scaler object removes the median and scales the data by the inter-quantile range. This form of scaling was performed to avoid excessive influence of outliers in the scaling process.
  71. # %%
  72. X_raw_df = SNYCQ[SNYCQ_items].replace([np.inf, -np.inf], np.nan).dropna(axis=0)
  73. idx = X_raw_df.index
  74. X_scaled = RobustScaler().fit_transform(X_raw_df.values)
  75. X_scaled_df = pd.DataFrame(X_scaled,index=idx, columns=X_raw_df.columns)
  77. # %% [markdown]
  78. # We look at the distributions of the sNYCQ data before and after scaling
  80. # %%
  81. plot_dist = X_raw_df.hvplot.hist(title='SNYC-Q pre scaling') + X_scaled_df.hvplot.hist(title='SNYC-Q post scaling', shared_axes=False)
  82. hv.save(plot_dist, osp.join('figures', 'FigureXX-SNYCQ_histograms_pre_post_scaling.html'))
  83. plot_dist
  85. # %% [markdown]
  86. # ![Distributions of SNYQC items before and after scaling](./figures/FigureXX-SNYCQ_histograms_pre_post_scaling.png)
  88. # %% [markdown]
  89. # ***
  90. # # 3. Outlier Detection
  91. #
  92. # We relied on the Mahalanobis distance of each scan to the sample's mean in order to detect outlier scans. We set the threshold to the top 3% quantile
  94. # %%
  95. mcd = MinCovDet().fit(X_scaled_df.values)
  96. # squared Mahalanobis distances for the training set
  97. md2 = mcd.mahalanobis(X_scaled_df.values) if hasattr(mcd, "mahalanobis") else mcd.dist_
  99. # Threshold
  100. thr = np.quantile(md2, OUTLIER_Q)
  101. keep = md2 < thr
  103. # %%
  104. md2_df = pd.DataFrame(md2, index=idx)
  105. md2_df.hvplot(hover_cols=['Subject','Run'], title='Mahalanobis ditance',ylabel='Mahalanobis distance') *hv.HLine(thr).opts(line_width=0.5, line_dash='dashed', line_color='k')
  107. # %% [markdown]
  108. # Now that we have identified two outlier scans, we create a new version of the data where those scans have been removed (e.g., ```X_$$$$_kept```)
  110. # %%
  111. idx_kept = idx[keep]
  112. X_scaled_kept_df = X_scaled_df.loc[idx_kept] # Scaled data for scans not marked as outliers
  113. X_raw_kept_df = X_raw_df.loc[idx_kept] # Original data for scans not marked as outliers
  114. print(f"[Outliers] Flagged {(~keep).sum()} / {len(md2)}; keeping {keep.sum()}")
  116. # %% [markdown]
  117. # # 4. Correlation between SNYCQ items
  118. #
  119. # To explore the structure of in-scanner experience reports, we first computed the Pearson’s correlation between the 11 in-scanner experience items.
  121. # %%
  122. # Compute correlation matrix
  123. X_raw_kept_corr_df = X_raw_kept_df.corr()
  125. # %%
  126. # Estimate P-value matrix (Pearson)
  127. cols = X_raw_kept_df.columns
  128. pval_df = pd.DataFrame(np.zeros((len(cols), len(cols))),
  129. index=cols, columns=cols)
  131. for i, c1 in enumerate(cols):
  132. for j, c2 in enumerate(cols):
  133. if i <= j:
  134. r, p = pearsonr(X_raw_kept_df[c1], X_raw_kept_df[c2])
  135. pval_df.loc[c1, c2] = p
  136. pval_df.loc[c2, c1] = p
  138. # %%
  139. # Get clustering order from seaborn
  140. clustergrid = sns.clustermap(X_raw_kept_corr_df)
  141. plt.close()
  143. row_order = clustergrid.dendrogram_row.reordered_ind
  144. ordered = X_raw_kept_corr_df.index[row_order]
  146. corr_ord = X_raw_kept_corr_df.round(2).loc[ordered, ordered]
  147. pval_ord = pval_df.loc[ordered, ordered]
  149. # %%
  150. # Make sure names exist (used by show_results)
  151. corr_ord.index.name = 'index'
  152. corr_ord.columns.name = 'col'
  153. corr_ord.name = 'corr'
  155. pval_ord.index.name = 'index'
  156. pval_ord.columns.name = 'col'
  157. pval_ord.name = 'pval'
  159. # %%
  160. # Calcualte the number of unique entries to do Bonferroni correction
  161. n_comps = corr_ord.shape[0]*(corr_ord.shape[0]-1)/2
  163. # %%
  164. # Plot: correlation heatmap with bold black outline for pBonf < 0.05
  165. plot = show_correlations_with_statistics(
  166. data_val=corr_ord,
  167. data_pval=pval_ord,
  168. pval_thr=0.05/n_comps,
  169. clabel='Pearson correlation',
  170. height=700,
  171. width=800,
  172. cmap='RdBu_r',
  173. fontscale=1.5,
  174. clim=(-0.7, 0.7)
  175. )
  177. # %% [markdown]
  178. # ![Figure 01 - Panel B](./figures/Figure01_B-SNYCQcorr.png)
  179. #
  180. # ### Saving Publication ready format data and figure panel
  181. #
  182. # We will save now the panel as HTML in a way that we can then export to PDF via the "Print..." option of the browser
  184. # %%
  185. import holoviews as hv
  187. from bokeh.io import save
  188. from bokeh.models.plots import Plot
  189. from bokeh.resources import INLINE
  191. hv.extension("bokeh")
  193. def svg_backend(hv_plot, element):
  194. hv_plot.state.output_backend = "svg"
  196. svg_plot = plot.opts(hooks=[svg_backend])
  198. bokeh_obj = hv.render(svg_plot, backend="bokeh")
  200. # Extra safety for layouts / overlays
  201. if isinstance(bokeh_obj, Plot):
  202. bokeh_obj.output_backend = "svg"
  204. for p in bokeh_obj.select({"type": Plot}):
  205. p.output_backend = "svg"
  207. save(
  208. bokeh_obj,
  209. filename="./figures/Figure01_B-SNYCQcorr.html",
  210. resources=INLINE,
  211. title="Figure01_B",
  212. )
  214. # %%
  215. corr_ord.to_csv('./source_data_files/figure_01_b.csv',float_format='%.2f')
  216. print('Saving source data for Figure 01B: SNYCQ correlation matrix with significant correlations outlined in black [./source_data_files/figure_01_b.csv]')
  218. # %% [markdown]
  219. # ***
  220. # # 4. Dimensionality Reduction with T-SNE
  221. # ## 4.1. Hyper-parameter Optimization
  222. # Two key hyper-parameters of the T-SNE algorithm are dimensionality and perplexity. We relied on [trustworthiness](https://scikit-learn.org/stable/modules/generated/sklearn.manifold.trustworthiness.html)—an estimate of how well a given embedding preserves local distances—to select these two hyper-parameters on a data-driven manner. The explored hyper-parameter space was:
  223. #
  224. # | Hyper-parameter | Space Explored |
  225. # |:----------------|:---------------|
  226. # |dimensionality| {1,2,3}|
  227. # | perplexity | {5,9,10,13,15,18,20,30,40,50} |
  229. # %%
  230. def choose_tsne_perplexity(X, random_state=RANDOM_STATE,n_components=[1,2,3]):
  231. n = X.shape[0]
  232. # Theoretical upper bound: perplexity < (n - 1); practical upper ~ n/3
  233. max_perp = max(5, min(60, (n - 1) // 3))
  234. base = min(max_perp, max(5, int(round(n / 50)))) # ~ n/50 baseline, clamped [5, 60]
  235. # Candidate set around 'base' plus some classics
  236. cands = sorted(set([
  237. 5, 10, 15, 20, 30, 40, 50,
  238. base, int(base*1.5), int(base*2)
  239. ]))
  240. cands = [p for p in cands if 5 <= p <= max_perp and p < (n - 1)]
  241. scores = []
  242. for nc in tqdm(n_components, position=0, desc='Num Components'):
  243. for p in tqdm(cands, position=1, desc='Perplexity', leave=False):
  244. tsne_tmp = TSNE(n_components=nc, perplexity=p, learning_rate='auto', init='pca',
  245. n_iter_without_progress=600, early_exaggeration=12.0, random_state=random_state)
  246. Ytmp = tsne_tmp.fit_transform(X)
  247. tw = trustworthiness(X, Ytmp, n_neighbors=10, metric='euclidean')
  248. scores.append((nc,p, tw))
  249. # pick the perplexity with max trustworthiness; break ties by preferring mid-range
  250. scores.sort(key=lambda t: (-t[2], abs(t[0] - 30)))
  251. return scores[0][0], scores[0][1], pd.DataFrame(scores, columns=["num_components","perplexity", "trustworthiness"])
  253. best_nc, best_perp, tw_table = choose_tsne_perplexity(X_scaled_kept_df.values, RANDOM_STATE)
  255. # %% [markdown]
  256. # Now we show what were the values that maximized the trustworthiness of the T-SNE embeddings
  258. # %%
  259. print(f"[t-SNE tuning] Selected perplexity = {best_perp}")
  260. print(f"[t-SNE tuning] Selected dimensionaity = {best_nc}")
  261. print(f"[t-SNE tuning] Trustworthiness = {tw_table['trustworthiness'].max()}")
  263. # %% [markdown]
  264. # ## 4.2. Compute T-SNE with optimal dimensionality and perplexity
  266. # %%
  267. tsne = TSNE(
  268. n_components=best_nc, perplexity=best_perp, learning_rate='auto', init='pca',
  269. n_iter_without_progress=1000, early_exaggeration=12.0, random_state=RANDOM_STATE
  270. )
  271. Y = tsne.fit_transform(X_scaled_kept_df.values)
  272. if best_nc == 2:
  273. emb = pd.DataFrame(Y, index=idx_kept, columns=["TSNE1", "TSNE2"])
  274. else:
  275. emb = pd.DataFrame(Y, index=idx_kept, columns=["TSNE1", "TSNE2","TSNE3"])
  277. # %% [markdown]
  278. # ## 4.3. Compute Biplot Arrows inidicating directions of maximal variance for each SNYCQ item
  279. #
  280. # To aid with the interpretation of how T-SNE dimensions relate to the original items in the SNYC survey, we modeled each of the SNYCQ items as a linear function of the three T-SNE coordinates. We then used these coefficients to draw biplot arrows depicting the directions of maximal change for each SNYCQ item in the 3D T-NSE embedding space.
  282. # %%
  283. def compute_biplot_arrows(emb_df, features_df, feature_names):
  284. """
  285. For each feature f: fit f ~ a*X + b*Y on the 2D embedding; return direction (a,b) and R^2.
  286. features_df should be standardized/robust-scaled (we use Xr_kept).
  287. """
  288. Xy = emb_df.values # shape (n,2)
  289. n_components = Xy.shape[1]
  290. arrows = []
  291. reg = LinearRegression()
  292. for j, f in enumerate(feature_names):
  293. y = features_df.iloc[:, j].values
  294. reg.fit(Xy, y)
  295. if n_components == 2:
  296. a, b = reg.coef_
  297. r2 = reg.score(Xy, y)
  298. arrows.append((f, a, b, r2))
  299. elif n_components == 3:
  300. a,b,c = reg.coef_
  301. r2 = reg.score(Xy,y)
  302. arrows.append((f,a,b,c,r2))
  303. else:
  304. print('++ ERROR: This function only works with 2D and 3D embeddings')
  305. return None
  306. if n_components==2:
  307. arrows_df = pd.DataFrame(arrows, columns=["feature", "beta_x", "beta_y", "R2"]).sort_values("R2", ascending=False)
  308. elif n_components == 3:
  309. arrows_df = pd.DataFrame(arrows, columns=["feature", "beta_x", "beta_y", "beta_z","R2"]).sort_values("R2", ascending=False)
  310. return arrows_df
  312. # %%
  313. tsne_arrows = compute_biplot_arrows(emb,X_raw_kept_df,SNYCQ_items)
  314. # Scaling so that they are clearly visible in the plot
  315. tsne_arrows['beta_x'] = tsne_arrows['beta_x'] * 5
  316. tsne_arrows['beta_y'] = tsne_arrows['beta_y'] * 5
  317. tsne_arrows['beta_z'] = tsne_arrows['beta_z'] * 5
  319. # %% [markdown]
  320. # ## 4.4. Plot the T-SNE embedding colored by one SNYCQ item
  322. # %%
  323. # Create a new DF with both th original values and the dimensions in T-SNE (for plotting purposes)
  324. emb_plus = pd.concat([emb, X_raw_kept_df], axis=1)
  326. # %%
  327. tsne_camera_object = dict(
  328. center=dict(x=6.661338147750939e-16, y=-1.942890293094024e-16, z=8.881784197001252e-16),
  329. eye=dict(x=0.4505864572659788, y=1.5176226342505497, z=-1.6428294069570146),
  330. up=dict(x=-0.7264163680411282, y=0.05191238226958195, z=0.6852914451596732)
  331. )
  332. # Create the scatter trace
  333. scatter_trace = go.Scatter3d(
  334. x=emb_plus['TSNE1'],
  335. y=emb_plus['TSNE2'],
  336. z=emb_plus['TSNE3'],
  337. mode='markers',
  338. marker=dict(
  339. size=5, # Fixed small size
  340. color=emb_plus['People'], # Color mapped to 'People'
  341. colorscale='viridis', # Viridis colormap
  342. opacity=0.8,
  343. colorbar=dict(title='People') # Optional colorbar
  344. ),
  345. text=[
  346. "<br>".join(f"{col}: {row[col]}" for col in emb_plus.columns if col not in ['TSNE1', 'TSNE2', 'TSNE3'])
  347. for _, row in emb_plus.iterrows()
  348. ],
  349. hoverinfo='text'
  350. )
  352. # === 2. Arrows from (0, 0, 0) to (beta_x, beta_y, beta_z) ===
  354. arrow_traces = []
  356. for i, row in tsne_arrows.iterrows():
  357. # Line (arrow body)
  358. if row['feature'] == 'People':
  359. arrow_color = 'red'
  360. else:
  361. arrow_color = 'black'
  362. arrow_trace = go.Scatter3d(
  363. x=[0, row['beta_x']],
  364. y=[0, row['beta_y']],
  365. z=[0, row['beta_z']],
  366. mode='lines',
  367. line=dict(
  368. color=arrow_color,
  369. width=max(1, row['R2'] * 10), # Scale width by R2
  370. ),
  371. showlegend=False
  372. )
  374. # Text label at arrow tip
  375. label_trace = go.Scatter3d(
  376. x=[row['beta_x']],
  377. y=[row['beta_y']],
  378. z=[row['beta_z']],
  379. mode='text',
  380. text=[row['feature']],
  381. textposition='top center',
  382. showlegend=False
  383. )
  385. arrow_traces.extend([arrow_trace, label_trace])
  388. # === 3. Combine all traces ===
  390. fig = go.Figure(data=[scatter_trace] + arrow_traces)
  392. # === 4. Update layout ===
  394. fig.update_layout(
  395. scene=dict(
  396. xaxis_title='TSNE1',
  397. yaxis_title='TSNE2',
  398. zaxis_title='TSNE3',
  399. aspectmode='data',
  400. camera=tsne_camera_object
  402. ),
  403. title='3D t-SNE Plot with Feature Arrows',
  404. margin=dict(l=0, r=0, b=0, t=40),
  405. height=700,
  406. showlegend=False
  407. )
  409. fig.show()
  412. # %%
  413. fig.write_html(osp.join('figures', 'FigureXX-SNYCQtsne_colorby_People.html'))
  415. # %% [markdown]
  416. # ![TSNE Map colored by the People question](./figures/FigureXX-SNYCQtsne_colorby_People.png)
  418. # %% [markdown]
  419. # ***
  420. # # 5. Clustering Analyses
  421. #
  422. # ## 5.1 Apply K-means and Gaussian Mixture Modeling to the data (K=2)
  423. #
  424. # To separate scans into two sets with well-differentiated inner-experience, we decided to apply two different clustering algorithms to the in-experience data (11 questions) following robust scaling.
  425. #
  426. # >NOTE: To be clear, clustering was performed in the original 11D space, not the low dimensional space generated by T-SNE.
  427. #
  428. # Working with two different clustering methods allows to evaluate the robustness of results against clustering technique.
  429. #
  430. # The two selected methods were [K-Means](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html) and [Gaussian Mixture Modeling](https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html); both as implemented in the python library scikit-learn.
  431. #
  432. # K-Means was chosen as a representative hard-clustering method often used in the neuroimaging literature.
  433. #
  434. # Gaussian Mixture Modeling was chosen because of its soft-clustering nature (i.e., it provides membership probabilities) and because it allows non-spherical clusters. In both instances, we set k=2.
  435. #
  436. # The next cell computes the GMM clustering for K = 2
  438. # %%
  439. K = 2
  440. gm = GaussianMixture(n_components=K, covariance_type="spherical", random_state=RANDOM_STATE).fit(X_scaled_kept_df.values)
  441. proba = gm.predict_proba(X_scaled_kept_df.values)
  442. labels_gmm = proba.argmax(axis=1)
  444. # %% [markdown]
  445. # Now we do the same using KMeans for comparison
  447. # %%
  448. km = KMeans(n_clusters=K, random_state=RANDOM_STATE, n_init=10).fit(X_scaled_kept_df.values)
  449. labels_km = km.predict(X_scaled_kept_df.values)
  451. # %% [markdown]
  452. # ## 5.2. Cluster method comparison with ARI and Silhouette Index
  453. #
  454. # We now check for the consistency of the clustering results across both methods using the Asjusted Rand Index (ARI), and also for their quality, separately, using the Silhouette Index (SI)
  455. #
  456. # 1. Compute the SI for each clustering result
  458. # %%
  459. sil_gmm = silhouette_score(X_scaled_kept_df.values, labels_gmm)
  460. sil_km = silhouette_score(X_scaled_kept_df.values, labels_km)
  462. # %% [markdown]
  463. # 2. Compute the ARI comparing both methods
  465. # %%
  466. ari_km_gmm = adjusted_rand_score(labels_km, labels_gmm)
  468. # %% [markdown]
  469. # 3. Print the computed statistics
  471. # %%
  472. print(f"[GMM k=2 sph] silhouette={sil_gmm:.3f}")
  473. print(f"[K-M k=2 ] silhouette={sil_km:.3f}")
  475. print(f"[Agreement] ARI(KMeans vs GMM) = {ari_km_gmm:.3f}")
  476. print(f"[Sizes] GMM clusters: {np.bincount(labels_gmm)}")
  479. # %% [markdown]
  480. # Based on the SI, we decided to move forward with the GMM solustion, which we explore in further detail.
  482. # %% [markdown]
  483. # ## 5.3. Bootstraing Analysis for GMM
  485. # %%
  486. def subsample_labels(model, X, frac=0.8, seed=0):
  487. rng = np.random.RandomState(seed)
  488. n = X.shape[0]
  489. take = np.sort(rng.choice(n, int(frac*n), replace=False)) # sort for convenience
  490. # clone model with same hyperparams
  491. m2 = type(model)(**model.get_params())
  492. if hasattr(m2, "random_state"):
  493. m2.random_state = seed
  494. # (KMeans has fit_predict; GMM needs fit + predict)
  495. if hasattr(m2, "fit_predict"):
  496. labels = m2.fit_predict(X[take])
  497. else:
  498. m2.fit(X[take]); labels = m2.predict(X[take])
  499. return take, labels
  501. def bootstrap_ari(model, X, n_runs=20, frac=0.8, base_seed=0):
  502. aris = []
  503. for s in range(n_runs):
  504. i1, l1 = subsample_labels(model, X, frac=frac, seed=base_seed + 2*s)
  505. i2, l2 = subsample_labels(model, X, frac=frac, seed=base_seed + 2*s + 1)
  506. # align to the same samples, same order
  507. common = np.intersect1d(i1, i2)
  508. # positions of the common indices in each subsample
  509. pos1 = np.searchsorted(i1, common)
  510. pos2 = np.searchsorted(i2, common)
  511. aris.append(adjusted_rand_score(l1[pos1], l2[pos2]))
  512. return float(np.mean(aris)), float(np.std(aris))
  514. gm_mean_ari, gm_sd_ari = bootstrap_ari(gm, X_scaled_kept_df.values, n_runs=20, frac=0.8, base_seed=100)
  515. print("Bootstrap ARI (GMM) mean±sd: %.2f +/- %.2f" %(gm_mean_ari, gm_sd_ari))
  517. # %% [markdown]
  518. # ***
  519. # ## 5.4. Detection of scans with ambigous cluster membership
  520. #
  521. # One additional bonus of GMM over K-means is that GMM not being a hard clustering algorithm, it outputs cluster membershup probabilities. We will use these to detect scans with ambiguous membership. Such scans will not be included in the population differences analyses later on.
  523. # %%
  524. proba_df = pd.DataFrame(proba, index=idx_kept, columns=['P(c1)','P(c2)'])
  525. proba_df.hvplot.hist('P(c1)', title='Probability Distribution for Membership in Cluster 1')
  527. # %%
  528. # store probabilities & an ambiguity flag
  529. p1 = proba[:, 1]
  530. ambiguity = np.maximum(1 - p1, p1) < (AMBIG_THRESH)
  532. # %% [markdown]
  533. # Save final cluster/sets labels in a new pandas dataframe: ```group_info_df```
  535. # %%
  536. # Store that information with clear labels in a pandas Dataframe
  537. group_info_df = pd.DataFrame(index=idx_kept, columns=['Set Label', 'Group Probability'])
  538. for i,scan in enumerate(idx_kept):
  539. if ambiguity[i]:
  540. group_info_df.loc[scan,"Set Label"] = "Ambiguous"
  541. group_info_df.loc[scan,"Group Probability"] = np.max([p1[i],1-p1[i]])
  542. else:
  543. if p1[i] > 1 - p1[i]:
  544. group_info_df.loc[scan,"Set Label"] = "Set B"
  545. group_info_df.loc[scan,"Group Probability"] = p1[i]
  546. else:
  547. group_info_df.loc[scan,"Set Label"] = "Set A"
  548. group_info_df.loc[scan,"Group Probability"] = 1 - p1[i]
  549. group_info_df['Set Label'].value_counts()
  551. # %% [markdown]
  552. # Add cluster membership information to the TSNE embedding information, and save two versions to disk: one with the original SNYCQ items, one with their scaled values.
  554. # %%
  555. emb_plus = pd.concat([emb_plus, group_info_df], axis=1)
  556. emb_plus.to_csv(osp.join(RESOURCES_SNYCQ_DIR, 'SNYCQ_tsne_embeddings_plus.csv'))
  557. print('++ INFO: Saved t-SNE embeddings with original features and group info to CSV: %s' % osp.join(RESOURCES_SNYCQ_DIR, 'SNYCQ_tsne_embeddings_plus.csv'))
  558. emb_plus_scaled = pd.concat([emb, X_scaled_kept_df, group_info_df], axis=1)
  559. emb_plus_scaled.to_csv(osp.join(RESOURCES_SNYCQ_DIR, 'SNYCQ_tsne_embeddings_plus_scaled.csv'))
  560. print('++ INFO: Saved t-SNE embeddings with scaled features and group info to CSV: %s' % osp.join(RESOURCES_SNYCQ_DIR, 'SNYCQ_tsne_embeddings_plus_scaled.csv'))
  562. # %% [markdown]
  563. # ## 5.5. Plot the T-SNE embedding again, but this time with scans colored according to set membership
  565. # %%
  566. tsne_camera_object = dict(
  567. center=dict(x=0.0, y=0.0, z=0.0),
  568. eye=dict(x=1.25, y=-1.25, z=1.25),
  569. up=dict(x=0.0, y=0.0, z=0.0)
  570. )
  571. scatter_traces = []
  573. group_colors = {
  574. "Set A": "#1f77b4", # light blue
  575. "Set B": "#ff7f0e", # orange
  576. "Ambiguous": "#ffffff" # white fill
  577. }
  579. for group, color in group_colors.items():
  580. group_data = emb_plus[emb_plus["Set Label"] == group]
  581. scatter_traces.append(
  582. go.Scatter3d(
  583. x=group_data["TSNE1"],
  584. y=group_data["TSNE2"],
  585. z=group_data["TSNE3"],
  586. mode='markers',
  587. name=group,
  588. marker=dict(
  589. size=5,
  590. color=color,
  591. opacity=0.9,
  592. line=dict(
  593. color='black' if group == "Ambiguous" else color,
  594. width=2 if group == "Ambiguous" else 0
  595. )
  596. ),
  597. text=[
  598. "<br>".join(
  599. f"{col}: {row[col]}"
  600. for col in emb_plus.columns if col not in ['TSNE1', 'TSNE2', 'TSNE3']
  601. )
  602. for _, row in group_data.iterrows()
  603. ],
  604. hoverinfo='text'
  605. )
  606. )
  607. # === 3. Combine all traces ===
  608. fig = go.Figure(data=scatter_traces + arrow_traces)
  610. fig.update_layout(
  611. scene=dict(
  612. xaxis_title='TSNE1',
  613. yaxis_title='TSNE2',
  614. zaxis_title='TSNE3',
  615. aspectmode='data',
  616. camera=tsne_camera_object
  617. ),
  618. legend=dict(
  619. title='Group:',
  620. itemsizing='constant',
  621. x=0.55, # move horizontally (1.0 is far right, 0.5 is middle)
  622. y=0.9, # move vertically (1.0 is top)
  623. xanchor='left',
  624. yanchor='top',
  625. bgcolor='rgba(255,255,255,0.7)', # optional background for readability
  626. bordercolor='rgba(0,0,0,0.2)',
  627. borderwidth=1
  628. ),
  629. title='3D t-SNE Plot with Feature Arrows',
  630. margin=dict(l=0, r=0, b=0, t=40),
  631. height=700,
  632. showlegend=True
  633. )
  635. fig.show()
  637. # %%
  638. fig.write_html(osp.join('figures', 'Figure01_K_SNYCQtsneWclusters.html'))
  640. # %% [markdown]
  641. # ![TSNE with clusters](./figures/Figure01_K-SNYCQtnseWclusters.png)
  643. # %% [markdown]
  644. # ### Save Publication Ready Figure Panel and Data source
  646. # %%
  647. fig.write_image(osp.join('figures', 'Figure01_K_SNYCQtsneWclusters.svg'))
  649. # %%
  650. emb_plus[['TSNE1', 'TSNE2', 'TSNE3','Set Label']].to_csv('./source_data_files/figure_01_k.csv', float_format='%.3f', index=False)
  652. # %% [markdown]
  653. # ***
  654. # ## 6. Explore how Sets A and B differ in terms of SNYCQ items, vigilance, head motion and basic demographics
  655. #
  656. # We will seek for potential differences across groups in the following variables:
  657. #
  658. # * All entries in the SNYCQ, including wakefulness
  659. # * Mean head motion
  660. # * Age distribution
  661. # * Gender distribution
  662. #
  663. # We will do this using Cohen's d, MannWhitney tests and Wilconxon tests
  665. # %% [markdown]
  666. # ## 6.1. Examination of differences in age distribution
  668. # %%
  669. # Load Demographic Data
  670. demographics = pd.read_csv(ORIG_DEMO_PATH, index_col=0,sep='\t')
  671. demographics = demographics.loc[list(SBJs)]
  672. # Load Demographic Data
  673. #demographics = pd.read_csv(osp.join(RESOURCES_SNYCQ_DIR,'participants_post_motion_QA.csv'), index_col=0)
  674. # Extract information about age
  675. age_per_scan = pd.DataFrame(index=idx_kept, columns=['Set Label','Age (5-year bins)'])
  676. for sbj,run in tqdm(idx_kept):
  677. age_range = demographics.loc[sbj,'age (5-year bins)']
  678. group_label = emb_plus.loc[(sbj,run),'Set Label']
  679. age_per_scan.loc[(sbj,run),'Age (5-year bins)'] = age_range
  680. age_per_scan.loc[(sbj,run),'Set Label'] = group_label
  681. # Remove entries for ambiguous scans
  682. age_per_scan = age_per_scan[age_per_scan['Set Label']!='Ambiguous']
  683. # Get counts of scas in each age range
  684. age_counts_per_group = age_per_scan.groupby('Set Label').value_counts()
  685. age_counts_per_group = age_counts_per_group.infer_objects()
  687. # Prepare Dataframe for plotting with hvplot
  688. age_counts_per_group = age_counts_per_group.reset_index()
  689. age_counts_per_group.columns = ['Group','Age Range','# Scans']
  690. age_counts_per_group.replace({'Set A':'A', 'Set B':'B'}, inplace=True)
  691. age_counts_per_group['color'] = '#ffffff'
  692. age_counts_per_group.loc[age_counts_per_group['Group']=='A','color'] = group_colors['Set A']
  693. age_counts_per_group.loc[age_counts_per_group['Group']=='B','color'] = group_colors['Set B']
  694. age_counts_per_group = age_counts_per_group.infer_objects()
  695. age_counts_per_group = age_counts_per_group.sort_values(by='Age Range', ascending=True)
  697. # %%
  698. A = age_counts_per_group.set_index(['Group','Age Range']).loc['A',:]['# Scans']
  699. B = age_counts_per_group.set_index(['Group','Age Range']).loc['B',:]['# Scans']
  700. W, w_p = wilcoxon(A,B, alternative='two-sided', method='exact')
  701. print('++ AGE ACROSS SETS: Wilcoxon = %.2f (p = %.2f)' % (W,w_p))
  704. # %%
  705. # Generate graph that will get later added to a Grid with information about all variables
  706. age_bar_plot = age_counts_per_group.hvplot.bar(x='Age Range',by='Group', alpha=0.5, xlabel='Age',cmap=["#1f77b4","#ff7f0e"]).opts(toolbar=None, xrotation=90, width=250, height=200, fontscale=1)
  708. # %% [markdown]
  709. # ![Age per set](./figures/Figure01_I-AgePerSet.png)
  711. # %% [markdown]
  712. # ### Save publication ready and source data for Figure 01I
  714. # %%
  715. age_counts_per_group.to_csv('./source_data_files/figure_01_i.csv', index=False)
  717. # %%
  718. svg_plot = age_bar_plot.opts(hooks=[svg_backend])
  720. bokeh_obj = hv.render(svg_plot, backend="bokeh")
  722. # Extra safety for layouts / overlays
  723. if isinstance(bokeh_obj, Plot):
  724. bokeh_obj.output_backend = "svg"
  726. for p in bokeh_obj.select({"type": Plot}):
  727. p.output_backend = "svg"
  729. save(
  730. bokeh_obj,
  731. filename="./figures/Figure01_I-AgePerSet.html",
  732. resources=INLINE,
  733. title="Figure01_I",
  734. )
  736. # %% [markdown]
  737. # ## 6.2. Examination of differneces in gender distribution
  739. # %%
  740. # Extract information about age
  741. sex_per_scan = pd.DataFrame(index=idx_kept, columns=['Set Label','Sex'])
  742. for sbj,run in tqdm(idx_kept):
  743. sex = demographics.loc[sbj,'gender']
  744. if sex == 'M':
  745. sex = 'Male'
  746. else:
  747. sex = 'Female'
  748. group_label = emb_plus.loc[(sbj,run),'Set Label']
  749. sex_per_scan.loc[(sbj,run),'Sex'] = sex
  750. sex_per_scan.loc[(sbj,run),'Set Label'] = group_label
  751. # Remove entries for ambiguous scans
  752. sex_per_scan = sex_per_scan[sex_per_scan['Set Label']!='Ambiguous']
  753. # Get counts of scas in each age range
  754. sex_counts_per_group = sex_per_scan.groupby('Set Label').value_counts()
  755. sex_counts_per_group = sex_counts_per_group.infer_objects()
  756. sex_counts_per_group
  758. # %%
  759. sex_bar_plot = sex_counts_per_group.hvplot.bar(stacked=True, xlabel='', legend='top_left', title='', ylabel='# Scans', color=['white','gray']).opts(toolbar=None, width=250, height=200, fontscale=1)
  760. #hv.save(sex_bar_plot, osp.join('figures', 'Figure01_J-SexPerSet.html'))
  762. # %% [markdown]
  763. # ![Sex per Set](./figures/Figure01_J-SexPerSet.png)
  764. #
  765. # ### Save Publication Ready and Data Source
  767. # %%
  768. svg_plot = sex_bar_plot.opts(hooks=[svg_backend])
  770. bokeh_obj = hv.render(svg_plot, backend="bokeh")
  772. # Extra safety for layouts / overlays
  773. if isinstance(bokeh_obj, Plot):
  774. bokeh_obj.output_backend = "svg"
  776. for p in bokeh_obj.select({"type": Plot}):
  777. p.output_backend = "svg"
  779. save(
  780. bokeh_obj,
  781. filename="./figures/Figure01_J-SexPerSet.html",
  782. resources=INLINE,
  783. title="Figure01_J",
  784. )
  786. # %%
  787. sex_counts_per_group.to_csv('./source_data_files/figure_01_j.csv')
  789. # %% [markdown]
  790. # ## 6.3. Examination of diffrences in head motion
  792. # %%
  793. # Load motion information for each scan
  794. mot_info = pd.read_csv(osp.join(RESOURCES_DINFO_DIR,'motion_confounds.csv'),index_col=['Subject','Run'])
  795. scans_in_A = emb_plus[emb_plus['Set Label'] == 'Set A'].index
  796. scans_in_B = emb_plus[emb_plus['Set Label'] == 'Set B'].index
  797. mot_A = mot_info.loc[scans_in_A,'Mean Rel Motion'].values
  798. mot_B = mot_info.loc[scans_in_B,'Mean Rel Motion'].values
  800. U, u_p = mannwhitneyu(mot_A,mot_B,alternative='two-sided')
  801. print('++ AGE ACROSS SETS: Mann-Whiteney U = %.2f (p = %.2f)' % (U,u_p))
  803. # %%
  804. mot_A = mot_info.loc[scans_in_A,'Mean Rel Motion']
  805. mot_B = mot_info.loc[scans_in_B,'Mean Rel Motion']
  806. overlay = mot_A.hvplot.hist(label='Set A', c=group_colors['Set A'], title='', width=250, height=200, alpha=0.5, shared_axes=False, bins=20, normed=True).opts(toolbar=None) * \
  807. mot_B.hvplot.hist(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False, bins=20, normed=True) * \
  808. mot_A.hvplot.kde(label='Set A', c=group_colors['Set A'], alpha=0.5, shared_axes=False) * \
  809. mot_B.hvplot.kde(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False)
  811. # %% [markdown]
  812. # ![Motion per set](./figures/Figure01_H-MotionPerSet.png)
  813. #
  814. # ### Save Publication ready figure and Source Data
  816. # %%
  817. svg_plot = overlay.opts(hooks=[svg_backend])
  819. bokeh_obj = hv.render(svg_plot, backend="bokeh")
  821. # Extra safety for layouts / overlays
  822. if isinstance(bokeh_obj, Plot):
  823. bokeh_obj.output_backend = "svg"
  825. for p in bokeh_obj.select({"type": Plot}):
  826. p.output_backend = "svg"
  828. save(
  829. bokeh_obj,
  830. filename="./figures/Figure01_H-MotionPerSet.html",
  831. resources=INLINE,
  832. title="Figure01_H",
  833. )
  835. # %%
  836. pd.concat([mot_A.reset_index(drop=True), mot_B.reset_index(drop=True)], axis=1, keys=['Set A','Set B']).to_csv('./source_data_files/figure_01_h.csv', index=False)
  838. # %% [markdown]
  839. # ## 6.4 Examination of Vigilance
  841. # %%
  842. vigilance = SNYCQ_wVigilance['Vigilance']
  843. vigilance.name = 'Wakefulness'
  844. scans_in_A = emb_plus[emb_plus['Set Label'] == 'Set A'].index
  845. scans_in_B = emb_plus[emb_plus['Set Label'] == 'Set B'].index
  846. vigilance_A = vigilance.loc[scans_in_A].values
  847. vigilance_B = vigilance.loc[scans_in_B].values
  849. U, u_p = mannwhitneyu(vigilance_A,vigilance_B,alternative='two-sided')
  850. print('++ VIGILANCE ACROSS SETS: Mann-Whiteney U = %.2f (p = %.2f)' % (U,u_p))
  852. # %%
  853. vigilance_A = vigilance.loc[scans_in_A]
  854. vigilance_B = vigilance.loc[scans_in_B]
  855. overlay = vigilance_A.hvplot.hist(label='Set A', c=group_colors['Set A'], title='', width=250, height=200, alpha=0.5, shared_axes=False, bins=[0,10,20,30,40,50,60,70,80,90,100], normed=True, legend=False).opts(toolbar=None) * \
  856. vigilance_B.hvplot.hist(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False, bins=[0,10,20,30,40,50,60,70,80,90,100], normed=True) * \
  857. vigilance_A.hvplot.kde(label='Set A', c=group_colors['Set A'], alpha=0.5, shared_axes=False) * \
  858. vigilance_B.hvplot.kde(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False)
  859. #hv.save(overlay, osp.join('figures', 'Figure01_G-VigilancePerSet.html'))
  861. # %% [markdown]
  862. # ![Wakefulness per set](./figures/Figure01_G-VigilancePerSet.png)
  863. # ### Save Publication Ready Panel and Source Data
  865. # %%
  866. svg_plot = overlay.opts(hooks=[svg_backend])
  868. bokeh_obj = hv.render(svg_plot, backend="bokeh")
  870. # Extra safety for layouts / overlays
  871. if isinstance(bokeh_obj, Plot):
  872. bokeh_obj.output_backend = "svg"
  874. for p in bokeh_obj.select({"type": Plot}):
  875. p.output_backend = "svg"
  877. save(
  878. bokeh_obj,
  879. filename="./figures/Figure01_G-VigilancePerSet.html",
  880. resources=INLINE,
  881. title="Figure01_G",
  882. )
  884. # %%
  885. pd.concat([vigilance_A.reset_index(drop=True), vigilance_B.reset_index(drop=True)], axis=1, keys=['Set A','Set B']).to_csv('./source_data_files/figure_01_g.csv', index=False)
  887. # %% [markdown]
  888. # ## 6.5. Examination of differences in SNYCQ items
  890. # %%
  891. def cohens_d(x0, x1):
  892. m0, m1 = np.nanmean(x0), np.nanmean(x1)
  893. s0, s1 = np.nanstd(x0, ddof=1), np.nanstd(x1, ddof=1)
  894. n0, n1 = np.sum(~np.isnan(x0)), np.sum(~np.isnan(x1))
  895. sp = np.sqrt(((n0-1)*s0**2 + (n1-1)*s1**2) / (n0 + n1 - 2))
  896. return (m1 - m0) / sp if sp > 0 else np.nan
  898. def calculate_stats_per_set(df, items, label_col='Set Label', method="bootstrap",
  899. B=2000, # bootstrap reps
  900. ci=95,random_state=123):
  902. rng = np.random.default_rng(random_state)
  903. alpha_low = (100 - ci) / 2.0
  904. alpha_high = 100 - alpha_low
  906. out = []
  907. for f in items:
  908. x0 = df.loc[df[label_col] == "Set A", f].dropna().values
  909. x1 = df.loc[df[label_col] == "Set B", f].dropna().values
  910. m0 = float(np.mean(x0)) if x0.size else np.nan
  911. m1 = float(np.mean(x1)) if x1.size else np.nan
  912. d = cohens_d(x0,x1)
  913. u, p = mannwhitneyu(x0,x1,alternative='two-sided')
  914. if method == "bootstrap":
  915. # nonparametric bootstrap of the mean per cluster
  916. if x0.size >= 2:
  917. boots0 = [np.mean(x0[rng.integers(0, x0.size, x0.size)]) for _ in range(B)]
  918. lo0, hi0 = np.percentile(boots0, [alpha_low, alpha_high])
  919. else:
  920. lo0 = hi0 = np.nan
  922. if x1.size >= 2:
  923. boots1 = [np.mean(x1[rng.integers(0, x1.size, x1.size)]) for _ in range(B)]
  924. lo1, hi1 = np.percentile(boots1, [alpha_low, alpha_high])
  925. else:
  926. lo1 = hi1 = np.nan
  928. elif method == "analytic":
  929. # mean ± z * SE; SE = s / sqrt(n)
  930. z = 1.96 if ci == 95 else None
  931. if z is None:
  932. from scipy.stats import norm
  933. z = float(norm.ppf(0.5 + ci/200.0))
  934. if x0.size >= 2:
  935. se0 = float(np.std(x0, ddof=1) / np.sqrt(x0.size))
  936. lo0, hi0 = m0 - z*se0, m0 + z*se0
  937. else:
  938. lo0 = hi0 = np.nan
  939. if x1.size >= 2:
  940. se1 = float(np.std(x1, ddof=1) / np.sqrt(x1.size))
  941. lo1, hi1 = m1 - z*se1, m1 + z*se1
  942. else:
  943. lo1 = hi1 = np.nan
  944. else:
  945. raise ValueError("method must be 'bootstrap' or 'analytic'.")
  947. out.append((f, d, u,p, m0, lo0, hi0, m1, lo1, hi1))
  949. out_df = pd.DataFrame(out, columns=["Item", "d(Set B - Set A)", "MW (U)","MW (p)","Set A (mean)", "lo (Set A)", "hi (Set A)", "Set B (mean)", "lo (Set B)", "hi (Set B)"]).set_index("Item")
  950. return out_df
  953. # %%
  954. non_ambiguous_scans = emb_plus[emb_plus['Set Label']!='Ambiguous'].index
  955. data_items = SNYCQ_wVigilance.loc[non_ambiguous_scans,:]
  956. data_labels = emb_plus.loc[non_ambiguous_scans,'Set Label']
  957. data = pd.concat([data_items,data_labels],axis=1)
  958. stats_per_set = calculate_stats_per_set(data,[c for c in data_items.columns if c !='Vigilance'],'Set Label')
  959. table_02 = stats_per_set.sort_values(by="d(Set B - Set A)", ascending=False).round(2)[['Set A (mean)','Set B (mean)','d(Set B - Set A)','MW (U)','MW (p)']]
  960. table_02
  962. # %%
  963. table_02.to_csv('./source_data_files/table_02.csv', float_format='%.2f')
  965. # %%
  966. layout = pn.GridBox(ncols=4)
  967. items_in_descending_d = stats_per_set.sort_values(by="d(Set B - Set A)", ascending=False).index
  968. set_a_scans = emb_plus[emb_plus['Set Label']=='Set A'].index
  969. set_b_scans = emb_plus[emb_plus['Set Label']=='Set B'].index
  970. for item in items_in_descending_d:
  971. bins = [0,10,20,30,40,50,60,70,80,90,100]
  972. this_pval = stats_per_set.loc[item,"MW (p)"]
  973. if this_pval >= 0.01:
  974. title = "Cohen d=%.2f | n.s." % stats_per_set.loc[item]["d(Set B - Set A)"]
  975. else:
  976. title = "Cohen d=%.2f | p < 0.01" % stats_per_set.loc[item]["d(Set B - Set A)"]
  977. overlay = SNYCQ.loc[set_a_scans,item].hvplot.hist(label='Set A', c=group_colors['Set A'], title=title, width=225, height=200, alpha=0.5, shared_axes=False, bins=bins, normed=True) * \
  978. SNYCQ.loc[set_b_scans,item].hvplot.hist(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False, bins=bins, normed=True) * \
  979. SNYCQ.loc[set_a_scans,item].hvplot.kde(label='Set A', c=group_colors['Set A'], title=title, width=300, height=200, alpha=0.5, shared_axes=False) * \
  980. SNYCQ.loc[set_b_scans,item].hvplot.kde(label='Set B', c=group_colors['Set B'], alpha=0.5, shared_axes=False)
  981. overlay = overlay.opts(show_legend=False,shared_axes=False, toolbar=None)
  982. layout.append(overlay)
  983. layout.save( osp.join('figures', 'Supplementary_Figure02.html'))
  985. # %%
  986. SNYCQ.loc[set_a_scans].to_csv('./source_data_files/suppfig_02_parta.csv',index=False)
  987. SNYCQ.loc[set_b_scans].to_csv('./source_data_files/suppfig_02_partb.csv',index=False)
  989. # %% [markdown]
  990. # ![Supplementary Figure 02](./figures/Supplementary_Figure02.png)
  992. # %% [markdown]
  993. # Provide the same information in more concise manner in the form of a radar plot
  995. # %%
  996. def plot_radar_means_with_ci(
  997. stats_df, # output of cluster_means_ci (indexed by feature)
  998. order=None, # optional ordering of features
  999. title="Per-cluster item means (with 95% CI)",
  1000. ci_label="95% CI"
  1001. ):
  1002. """
  1003. Draws radar with Cluster 0 & 1 mean lines and shaded CI bands.
  1004. """
  1005. # order features
  1006. if order is None:
  1007. features = list(stats_df.index)
  1008. else:
  1009. features = list(order)
  1010. stats_df = stats_df.loc[features]
  1012. # angles for axes + close the loop
  1013. angles = np.linspace(0, 2*np.pi, len(features), endpoint=False).tolist()
  1014. angles += angles[:1]
  1016. # extract arrays and close loops
  1017. m0 = stats_df["Set A (mean)"].to_numpy().tolist(); m0 += m0[:1]
  1018. m1 = stats_df["Set B (mean)"].to_numpy().tolist(); m1 += m1[:1]
  1019. lo0 = stats_df["lo (Set A)"].to_numpy().tolist(); lo0 += lo0[:1]
  1020. hi0 = stats_df["hi (Set A)"].to_numpy().tolist(); hi0 += hi0[:1]
  1021. lo1 = stats_df["lo (Set B)"].to_numpy().tolist(); lo1 += lo1[:1]
  1022. hi1 = stats_df["hi (Set B)"].to_numpy().tolist(); hi1 += hi1[:1]
  1024. # limits based on CI envelopes
  1025. all_vals = np.array((lo0 + hi0 + lo1 + hi1), dtype=float)
  1026. finite_vals = all_vals[np.isfinite(all_vals)]
  1027. if finite_vals.size:
  1028. rmin, rmax = float(np.nanmin(finite_vals)), float(np.nanmax(finite_vals))
  1029. else:
  1030. rmin, rmax = 0.0, 1.0
  1031. pad = max(1.0, 0.05 * (rmax - rmin))
  1033. # plot
  1034. fig = plt.figure(figsize=(7, 7))
  1035. ax = plt.subplot(111, polar=True)
  1037. # mean lines
  1038. l0, = ax.plot(angles, m0, linewidth=2, label="Set A")
  1039. l1, = ax.plot(angles, m1, linewidth=2, label="Set B")
  1041. # shaded CI polygons (build as upper path + reversed lower path)
  1042. # cluster 0
  1043. ang_np = np.array(angles)
  1044. poly0_ang = np.concatenate([ang_np, ang_np[::-1]])
  1045. poly0_rad = np.concatenate([np.array(hi0), np.array(lo0)[::-1]])
  1046. s0 = ax.fill(poly0_ang, poly0_rad, alpha=0.15, label=f"Set A {ci_label}")
  1048. # cluster 1
  1049. poly1_ang = np.concatenate([ang_np, ang_np[::-1]])
  1050. poly1_rad = np.concatenate([np.array(hi1), np.array(lo1)[::-1]])
  1051. s1 = ax.fill(poly1_ang, poly1_rad, alpha=0.15, label=f"Set B {ci_label}")
  1053. # axes & labels
  1054. ax.set_xticks(angles[:-1])
  1055. ax.set_xticklabels(features)
  1056. ax.set_ylim(rmin - pad, rmax + pad)
  1057. ax.set_title(title)
  1058. ax.legend(loc="upper right", bbox_to_anchor=(0.1, 1.0))
  1060. plt.tight_layout()
  1061. plt.show()
  1062. return fig
  1064. # %%
  1065. # Plot (keep your preferred order of spokes)
  1066. plot = plot_radar_means_with_ci(stats_per_set, order=list(items_in_descending_d),
  1067. title="SNYCQ per-cluster means with 95% CI")
  1069. # %% [markdown]
  1070. # ### Save Publication Ready Panel and Source Data for the Radar Plot
  1072. # %%
  1073. stats_per_set.to_csv('./source_data_files/figure_01_f.csv', float_format='%.2f')
  1075. # %%
  1076. plot.savefig('./figures/Figure01_F.svg', format='svg', bbox_inches='tight')
  1078. # %% [markdown]
  1079. # ***
  1080. #
  1081. # # Distribution of SNYCQ values (Supplementary Figure 1)
  1083. # %%
  1084. layout = None
  1085. for q in SNYCQ.columns:
  1086. plot = SNYCQ[q].reset_index(drop=True).hvplot.hist(bins=np.linspace(0,100,20), width=250, height=200, normed=True, ylabel='Density', fontsize=12) * SNYCQ[q].reset_index(drop=True).hvplot.kde()
  1087. if layout is None:
  1088. layout = plot
  1089. else:
  1090. layout = layout + plot
  1091. layout = layout.cols(3).opts(toolbar=None)
  1092. hv.save(layout, osp.join('figures', 'Supplementary_Figure01.html'))
  1094. # %% [markdown]
  1095. # ![Supplementary Figure 01](./figures/Supplementary_Figure01.png)
  1097. # %%
  1098. SNYCQ.to_csv('./source_data_files/suppfig_01.csv',index=False)
  1100. # %% [markdown]

S10_SNYCQ_TSNE_Clustering.ipynb at commit 2aa1bb1, under CC0-1.0 · at the source

Overview

Authors: Javier Gonzalez-Castillo1, Megan A. Spurney1,2, Ka Chun Lam3, Isabel S. Gephart1, Francisco Pereira3, Daniel A. Handwerker1, Julia W. Y. Kam4,5, Peter A. Bandettini1,6
ORCID iDs: Javier Gonzalez-Castillo, Ka Chun Lam, Daniel A. Handwerker, Julia W. Y. Kam, Peter A. Bandettini
  1. Section on Functional Imaging Methods, NIMH, NIH,Bethesda, MA USA
  2. Department of Psychology, Northwestern University,Evanston, IL USA
  3. Machine Learning Team, NIMH, NIH,Bethesda, MA USA
  4. Department of Psychology, University of Calgary,Calgary, AL Canada
  5. Hotchkiss Brain Institute, University of Calgary,Calgary, AL Canada
  6. Functional MRI Core, NIMH, NIH,Bethesda, MA USA
Institutions: National Institutes of Health (United States); Northwestern University (United States); University of Calgary (Canada)
Journal: Nature communications, volume 17, issue 1, article 8058
Dates: received 8 December 2025; accepted 16 June 2026; published online 27 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-74953-6 · PMID 42365014 · PMCID PMC13454257 · OpenAlex W7166358421
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), human (organism), clinical / translational (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging
Keywords: Consciousness, Cognitive neuroscience, Biomarkers
MeSH: Brain*, Magnetic Resonance Imaging*, Rest*, Thinking*, Adult, Brain Mapping, Female, Humans, Male, Retrospective Studies (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institute of Mental Health (ZIAMH002783, ZICMH002968, ZICMH002884)
Citations: not cited yet (Europe PMC); 99 references in the paper

Abstract

Resting-state fMRI (rsfMRI) scans—acquired in the absence of experimentally controlled stimuli or task demands—are widely used to identify aberrant patterns of functional connectivity (FC) in clinical populations. To minimize interpretational uncertainty, researchers routinely control for across-cohort disparities in age, gender, comorbidities, and head motion. Yet, studies rarely consider the possibility that systematic differences in inner experience (i.e., how subjects think and feel during the scan) directly affect FC measures. Here, using an rsfMRI dataset comprising 469 scans with retrospective experiential annotations, we show that summary descriptors of in-scanner experience are reproducible across visits and subject-specific, consistent with trait-like characteristics. We further show that widespread significant differences in FC are observed between scans that are associated with different reported experiential profiles, and that FC can predict specific experiential dimensions with performance comparable to that reported for demographic, cognitive, and clinical variables. Together, these findings highlight the key role that in-scanner experience should play when interpreting FC in the context of rsfMRI. Given that the available experiential measures are retrospective summaries, these results speak to stable experiential tendencies rather than potential moment-to-moment, state-dependent relationships between ongoing experience and concurrent brain activity.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

nimh-sfim/fc_introspection

License: CC0-1.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 2aa1bb1b653111e0486001780837e9734d90eaf3, 1 June 2026
Languages: Python (36), Jupyter (24), Shell (9), MATLAB (3)
Size: 326 files, 72 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 24 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (33 files), NumPy (17 files), AFNI (11 files), Matplotlib (10 files), seaborn (10 files), scikit-learn (5 files), SciPy (5 files), ANTs (4 files), xarray (4 files), FreeSurfer (2 files), FSL (2 files), Plotly (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
47 files
At the source: github.com/nimh-sfim/fc_introspection

Code availability

Code used in this study is available at https://github.com/nimh-sfim/fc_introspection (release tagged as v4.0.0).

Reproduced under the paper's license (CC BY), from the paper cited above.

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;
  • 45 scripts, each with its path and the digest of its content;
  • 15 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

Datasets cited

Data availability

This study is based on publicly available data as described in Mendes et al35 Imaging data can be accessed at ftp://ftp.gwdg.de/pub/misc/MPILeipzig_Mind-Brain-Body. Behavioral data is available at https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/VMJ6NV. Additional access locations can be found in the Mendes et al. publication. In addition, Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

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, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 3 keywords, 10 MeSH terms, 1 funder, 93 references.

Cite

This paper

Gonzalez-Castillo, J., Spurney, M. A., Lam, K. C., Gephart, I. S., Pereira, F., Handwerker, D. A., Kam, J. W. Y., & Bandettini, P. A. (2026). In-scanner thoughts contribute to resting-state functional connectivity. Nature communications, 17(1), 8058. https://doi.org/10.1038/s41467-026-74953-6

BibTeX

@article{gonzalezcastillo2026scanner,
author = {Gonzalez-Castillo, Javier and Spurney, Megan A. and Lam, Ka Chun and Gephart, Isabel S. and Pereira, Francisco and Handwerker, Daniel A. and Kam, Julia W. Y. and Bandettini, Peter A.},
title = {{In-scanner thoughts contribute to resting-state functional connectivity}},
journal = {Nature communications},
year = {2026},
month = jun,
volume = {17},
number = {1},
pages = {8058},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-74953-6},
url = {https://doi.org/10.1038/s41467-026-74953-6},
pmid = {42365014},
pmcid = {PMC13454257}
}

RIS

TY - JOUR
AU - Gonzalez-Castillo, Javier
AU - Spurney, Megan A.
AU - Lam, Ka Chun
AU - Gephart, Isabel S.
AU - Pereira, Francisco
AU - Handwerker, Daniel A.
AU - Kam, Julia W. Y.
AU - Bandettini, Peter A.
TI - In-scanner thoughts contribute to resting-state functional connectivity
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/06/27
VL - 17
IS - 1
SP - 8058
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-74953-6
UR - https://doi.org/10.1038/s41467-026-74953-6
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-74953-6",
"type": "article-journal",
"title": "In-scanner thoughts contribute to resting-state functional connectivity",
"container-title": "Nature communications",
"author": [
{
"family": "Gonzalez-Castillo",
"given": "Javier"
},
{
"family": "Spurney",
"given": "Megan A."
},
{
"family": "Lam",
"given": "Ka Chun"
},
{
"family": "Gephart",
"given": "Isabel S."
},
{
"family": "Pereira",
"given": "Francisco"
},
{
"family": "Handwerker",
"given": "Daniel A."
},
{
"family": "Kam",
"given": "Julia W. Y."
},
{
"family": "Bandettini",
"given": "Peter A."
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8058",
"DOI": "10.1038/s41467-026-74953-6",
"PMID": "42365014",
"PMCID": "PMC13454257",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-74953-6",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
27
]
]
}
}

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.7554/elife.104053 [code]
Brain cognition gaps reveal associations with dopamine and factors related to brain health through artificial intelligence prediction of functional connectome.
Journal: eLife
In common: scikit-learn, pandas, SciPy, 1 other tool, 10 references
[2] doi:10.1002/epi.70323 [code]
Individual-specific resting-state networks predict language dominance in drug-resistant epilepsy.
Journal: Epilepsia
In common: AFNI, ANTs, FreeSurfer, 5 other tools, fMRI, 4 references
[3] 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: ANTs, FreeSurfer, FSL, 7 other tools, fMRI, 3 references
[4] 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: ANTs, FreeSurfer, FSL, 7 other tools, fMRI, 3 references
[5] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: AFNI, ANTs, FreeSurfer, 5 other tools, 3 references
[6] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: AFNI, ANTs, FreeSurfer, 7 other tools, fMRI, 2 references
[7] doi:10.1016/j.neuron.2026.04.011 [code]
Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.
Journal: Neuron
In common: AFNI, ANTs, FreeSurfer, 5 other tools, fMRI, 3 references
[8] doi:10.1162/imag.a.1336 [code]
Brain network dynamics in the wake-sleep transition reorganize according to task engagement.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: 7 references, author Javier Gonzalez-Castillo
[9] doi:10.1038/s41398-026-04025-2 [code]
Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: AFNI, ANTs, FreeSurfer, 7 other tools, 2 references
[10] doi:10.1016/j.xcrm.2026.102943 [code]
Parent-of-origin effects in Alzheimer's liability dissociate neurocognitive and cardiovascular traits in at-risk individuals.
Journal: Cell reports. Medicine
In common: AFNI, ANTs, FreeSurfer, 8 other tools, clinical / translational

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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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