OSCR

Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease.

Code ↔ Paper

12 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 12 matches
  1. [1] § Methods › Dimensional symptom analysis ↔ Fingerprints_Cocuzza.py, lines 698–750 · score 0.92 · score_samples, log likelihood, nbClust, distance matrix, linkage, scikit-learn
  2. [2] § Methods › Data collection › Study outline ↔ Data_Splitting_Cocuzza.py, lines 213–275 · score 0.91 · depression NOS, social anxiety, eating disorder, panic disorder, OCD, cyclothymia
  3. [3] § Methods › Dimensional symptom analysis ↔ Fingerprints_Cocuzza.py, lines 89–222 · score 0.83 · childhood trauma questionnaire, raw score, redundancies, subscale, curated, inflated
  4. [4] § Methods › Dimensional symptom analysis ↔ Fingerprints_Cocuzza.py, lines 332–415 · score 0.83 · Anderson Darling, Kolmogorov Smirnov, Shapiro Wilk, Gaussianity, zero, transformations
  5. [5] § Methods › Dimensional symptom analysis ↔ Fingerprints_Cocuzza.py, lines 473–597 · score 0.76 · miss forest, min max normalization, Random forest, imputed, space, binarizing
  6. [6] § Methods › Brain network dynamics analysis › Non-negative matrix factorization (NMF) ↔ NMF_Cocuzza.py, lines 223–274 · score 0.75 · feasibility, scikit-learn, computational, divergence, solver, evenly
  7. [7] § Results › Brain network dynamics are flattened in psychiatric illness ↔ NMF_Cocuzza.py, lines 513–581 · score 0.73 · encoding matrix, features matrix, brain regions, brain network, tuned, beta
  8. [8] § Results › Brain network dynamics are flattened in psychiatric illness ↔ NMF_Cocuzza.py, lines 167–221 · score 0.71 · upper triangle, configuration matrix, network assignment, fMRI, parcel, thresholded
  9. [9] § Methods › Data collection › Participants ↔ Data_Splitting_Cocuzza.py, lines 213–275 · score 0.67 · eating disorder, panic disorder, phobia, ADHD, anxiety, depression
  10. [10] § Methods › Brain network dynamics analysis › Non-negative matrix factorization (NMF) ↔ NMF_Cocuzza.py, lines 513–581 · score 0.62 · encoding matrix, brain regions, brain network, components, weights, coefficients
  11. [11] § Methods › Data collection › Participants ↔ NMF_Cocuzza.py, lines 167–221 · score 0.51 · fMRI, functional network, age, Connectome, validation, Brain
  12. [12] § Results › Symptom fingerprints: core dimensions of functioning are exhibited across health and psychiatric illness ↔ Fingerprints_Cocuzza.py, lines 698–750 · score 0.50 · log likelihood, Agglomerative, PCA, Ward, component, distance

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

Python · 776 lines · 47 KB · GPL-3.0 · 5 matches

  1. # C.V. Cocuzza, 2025. Symptom fingerprinting code for "Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease". Cocuzza et al., 2025.
  3. # The steps of this pipeline are organized into python functions that are possible to adapt to other datasets,
  4. # but the TCP dataset is quite unique (at the time of preparing this manuscript) in that it collected many more
  5. # behavioral, cognitive, and clinical measures (i.e., self-report surveys, cognitive batteries, clinical instruments, etc.)
  6. # than is typical. Thus, we've uploaded (to GitHub) de-identified behavioral data files to make implementing the functions here easier
  7. # (and hopefully easier to adapt to other datasets). The participant IDs match those openly available on OpenNeuro and NDA for
  8. # TCP dataset, but all PHI has been removed. Example usage is recommended throughout along with explanatory notes.
  10. # NOTE: the data in all_TCP_behav_concatenated_mapped_domains.pkl was pre-labelled and verified by coauthors
  11. # (see Methods of manuscript; and large Supplemental Table that also reflects this data)
  13. # WORKFLOW (corresponds to modular functions below; so see each function for more detailed notes):
  14. # 1) Organize and curate data: fingerprints_curate_Cocuzza()
  15. # 2) Transform measures toward Gaussianity: fingerprints_transform_Cocuzza(), cnd_project_imputation_R_pipe.Rmd
  16. # ** NOTE: this step is multi-part and partially involves RStudio. See cnd_project_imputation_R_pipe.Rmd
  17. # 3) Binarize 0-inflated variables: fingerprints_binarize_Cocuzza()
  18. # 4) Impute missing values/scores: fingerprints_imputation_Cocuzza()
  19. # 5) Estimate distance matrix (individual differences correlations of normalized and transformed measures): fingerprints_distance_Cocuzza()
  20. # 6) Input distance matrix into agglomerative hiearchical clustering to partition data into empirically-driven dimensions of functioning: fingerprints_clustering_Cocuzza()
  21. # ** NOTE: this step partially involves RStudio, see cnd_project_clustering_checks_pipe.Rmd
  22. # 7) Perform PCA on each cluster of measures (normalized and imputed): fingerprints_PCA_Cocuzza()
  24. # NOTE: steps 8 and 9 are TBA to GitHub.
  25. # 8) Use PCA results (i.e., weight by factor loadings) and pre-labellings to name clusters: fingerprints_cluster_naming_Cocuzza()
  26. # 9) Quantify multi-dimensional symptom profiles (i.e., symptom fingerprints) for each participant: fingerprints_profiles_Cocuzza()
  28. # NOTE: we were able to take many of the steps here *because* of the Transdiagnostic Connectome Project (TCP) dataset
  29. # having a relatively *very large* set of behavioral, cognitive, and clinical measures to use. That is to say,
  30. # if we excluded a given measure, it was very likely that there were still other measures to cover that proportion of variance in functioning.
  31. # Even with all the above steps, the final set of behavioral data included 110 variables.
  32. # This may not be the case (and is likely not) for other datasets, thus these steps should be adapted with caution and thorough consideration.
  34. '''
  35. # Example code to load behavioral data shared in GitHub repo
  36. # (copy and paste this section into jupyter notebook, python script, etc.)
  37. # dataHere_All can be used in first function below; see notes for all functions for usages (these are modular and in-order).
  39. import pandas as pd
  40. import numpy as np
  41. import pickle
  43. dirHere = '/directory/path/to/where/you/saved/scripts/and/data/' # CHANGE
  45. keysAll = np.load(dirHere + 'all_TCP_behav_concatenated_keys.npy',allow_pickle=True).copy()
  47. dataHere_All = pd.read_csv(dirHere + 'all_TCP_behav_concatenated.csv')
  49. with open(dirHere + 'all_TCP_behav_concatenated_mapped_domains.pkl', 'rb') as f:
  50. behavTestLegend = pickle.load(f)
  52. numMeasures = dataHere_All.shape[1]-1 # 1st column is subject ID
  53. print(f"Initial data dimensions: N = {dataHere_All.shape[0]} participants and {numMeasures} behavioral/clinical/cognitive measures.")
  54. '''
  56. # Other notes:
  58. # Using missForest to impute missing values in observed data
  59. # NOTE: can use missCompare in R to see if missForest is still best solution (TBA), but for now going straight to missForest
  60. # See here: https://github.com/rawat126/Feature-Engineering/blob/c523d48154d7808ee26214e3c70066a0981d634a/Miss%20Forest%20Imputation/Miss_forest_imputation.ipynb#L6
  61. # And here: https://ragvenderrawat.medium.com/miss-forest-imputaion-the-best-way-to-handle-missing-data-feature-engineering-techniques-2e6922e5cecb
  62. # Also see: https://scikit-learn.org/stable/modules/impute.html
  63. # Also see: https://ragvenderrawat.medium.com/miss-forest-imputaion-the-best-way-to-handle-missing-data-feature-engineering-techniques-2e6922e5cecb
  64. # See: https://www.analyticsvidhya.com/blog/2022/05/handling-missing-values-with-random-forest/
  66. ####################################################################################
  67. # IMPORTS (some may need to be installed)
  68. import sys
  69. import os
  70. #import subprocess
  71. #import time
  73. import numpy as np
  75. import pandas as pd
  76. import pickle
  77. import h5py as h5
  79. import scipy
  80. #import scipy.io as spio
  81. import scipy.stats as stats
  82. #from scipy import signal
  83. import sklearn
  84. import statistics as stat
  86. from sklearn.ensemble import RandomForestClassifier # for categorical data
  87. from sklearn.ensemble import RandomForestRegressor # for continuous data
  89. ####################################################################################
  90. def fingerprints_curate_Cocuzza(dataHere_All,keysAll,saveResults=False,saveDir=None,percentThresh=25,removeRedunancies=True,removeSparseSubjs=False,percentThresh_Subjs=50):
  92. '''
  93. Data organizing & curating: trim measures based on a series of related conditions:
  94. a) remove measures with > threshold % of subjects missing
  95. b) remove redundant meaures (e.g., having both t-scores and raw scores) -- this was manually coded for TCP and agreed upon by coauthors (and based on literature/best practices)
  96. c) remove subjects above a certain threshold of sparseness. This was done in the TCP data release manuscript but we did *not* do it here because downstream steps tended to exclude those anyway (see imputation function for how we accounted for this)
  98. INPUTS:
  99. dataHere_All: see top of script for how to load example data; from csv file with all behavioral data concatenated into dataframe
  100. keysAll: see top of script for how to load example data; from npy file with all keys for dataHere_All
  101. saveResults: Boolean w/ Default False; whether or not to also save out the curated dataframe (into .csv)
  102. saveDir: if saveResults=True; a directory path (include entire string) to save out results (csv) from curation steps
  103. percentThresh: subject-missingness threshold for measure removal, Default is 25
  104. removeRedunancies: Boolean w/ Default True; whether or not to remove redundant measures (recommended)
  105. removeSparseSubjs: Boolean w/ Default False; whether or not to remove sparse subjects (not recommended for this study)
  106. percentThresh_Subjs: If removeSparseSubjs is True (not recommended), this is the sparsity threshold to use. 50% (Default) is relatively lenient.
  108. OUTPUTS:
  109. dataHere_TrimFinal: new dataframe with measures trimmed; returned and also saved to [saveDir + 'dataHere_TrimFinal_BeforeGaussTest.csv'] when saveResults=True
  110. '''
  112. # a) Remove measures with more than a given % of subjects missing
  114. numMeasures = dataHere_All.shape[1]-1 # ignore 1st column, which is subject ID
  115. allKeys = dataHere_All.keys().to_numpy().copy()
  117. # Loop over measures and find those with >percentThresh% of subjects' missing data
  118. measuresToRemove = []
  119. measuresToRemove_Ixs = []
  120. for measureIx in range(numMeasures):
  121. thisMeasureStr = keysAll[1:][measureIx] # skip 1st column (subject ID)
  123. measureVec_AllSubjs = dataHere_All[thisMeasureStr].to_numpy().copy()
  124. numMissingSubjs = np.where(np.isnan(measureVec_AllSubjs))[0].shape[0]
  125. percentMissingSubjs = np.floor((numMissingSubjs / dataHere_All.shape[0]) * 100)
  127. if percentMissingSubjs > percentThresh:
  128. print(f"{thisMeasureStr} has {percentMissingSubjs}% subjects missing, flagging it for removal...")
  129. measuresToRemove.append(thisMeasureStr)
  130. measuresToRemove_Ixs.append(measureIx)
  132. measuresToRemove = np.asarray(measuresToRemove)
  133. measuresToRemove_Ixs = np.asarray(measuresToRemove_Ixs)
  135. dataHere_Trimmed = dataHere_All.copy()
  136. dataHere_Trimmed = dataHere_Trimmed.drop(measuresToRemove,axis=1)
  138. keys_Trimmed = allKeys[1:].copy()
  139. keys_Trimmed = np.delete(keys_Trimmed,measuresToRemove_Ixs)
  141. print(f"\nTrimming the dataset by {measuresToRemove.shape[0]} behavioral/clinical/cognitive measures that had > {percentThresh}% of subjects missing. "+
  142. f"New number of measures = {dataHere_Trimmed.shape[1]-1}.\n")
  144. # b) data curating continued from above: trim redundant variables
  146. # NOTE: important to avoid results being overly driven by 1 scale; but I'm making it optional
  147. # NOTE: these were decided by hand because it typically involved reviewing the available literature and finding the field standard/preference for that measure
  149. # NOTES on what was trimmed:
  150. # POMS, SHIPLEY: keep t-scores
  151. # YMRS: keep sub-scales and remove total (sub-scales likely offer more diverse coverage of behavioral/cognitive constructs / explain more variance)
  152. # BAPQ, BIS, MSPSS, TCI, POMS, COGFQ, RSRI: keep sub-scales and remove total
  153. # MADRS: keep total and remove subscales because the subscales are 0 inflated and would end up being binarized (which isn't ideal) later anyway
  154. # LRIFT: keep total and remove subscales because the subscales are 0 inflated and would end up being binarized (which isn't ideal) later anyway
  155. # PUM: keep total-via-average and remove total-via-sum
  156. # CTQ: removing to keep as a demographic variable/covariate (childhood trauma questionnaire is unique in that it's based on past experiences)
  157. # CRT: these variables end up being binarized (not ideal) and overfit the clustering solutions downstream
  159. dataHere_Trimmed_NonRedundant = dataHere_Trimmed.copy()
  160. if removeRedunancies:
  161. # This list corresponds to the notes above (re: "what was trimmed" and why)
  162. keysToRM = ['lrift_employment','lrift_household', 'lrift_school', 'lrift_work','lrift_spouse_relate', 'lrift_children_relate','lrift_relatives_relate',
  163. 'lrift_friends_relate','lrift_relationships', 'lrift_satisfaction', 'lrift_recreation',
  164. 'madrs_reduced_appetite', 'madrs_apparent_sadness','madrs_concentration_issues', 'madrs_inability_feel','madrs_lassitude',
  165. 'madrs_pessimism', 'madrs_reported_sadness','madrs_reduced_sleep', 'madrs_suicidal_thoughts','madrs_inner_tension',
  166. 'ctq_denial_validity', 'ctq_emo_abuse', 'ctq_emo_neglect','ctq_phys_abuse', 'ctq_phys_neglect', 'ctq_sex_abuse',
  167. 'POMS_anger','POMS_confusion','POMS_depression','POMS_fatigue','POMS_tension','POMS_vigour','POMS_total_raw',
  168. 'SHIPLEY_vocab','ymrs_total','bapq_total_avg','bis_total_sum','mspss_total','tci_total_sum','cogfq_total',
  169. 'crt_correct_first_3','crt_intuitive_first_3', 'crt_correct_all_5','crt_intuitive_all_5','rsri_total_average','pum_total_sum']
  171. dataHere_Trimmed_NonRedundant = dataHere_Trimmed_NonRedundant.drop(keysToRM,axis=1)
  172. print(f"Removing {len(keysToRM)} behavioral/clinical/cognitive measures that were redundant and/or circular with other measures. "+
  173. f"New number of measures = {dataHere_Trimmed_NonRedundant.shape[1]-1}.\n")
  175. # c) remove subjects with > a certain % of missing variables
  177. # NOTE: for the purposes of CND project we did NOT do this, but keeping here for other projects/uses/datasets. Reasons:
  178. # (1) subjects will be removed later anyway (only certain functional runs were of-interest for NMF/dynamics pipeline)
  179. # (2) we need as many subjects as possible for train / test / validation splits,
  180. # (3) so instead we opted for a well-validated imputation procedure (see later function)
  182. dataHere_TrimFinal = dataHere_Trimmed_NonRedundant.copy()
  184. if removeSparseSubjs:
  186. subjKeysAll = dataHere_Trimmed_NonRedundant['subjectkey'].to_numpy().copy()
  187. keysHere = dataHere_Trimmed_NonRedundant.keys().to_numpy().copy()
  189. subjsToRemove_IDs = []
  190. subjsToRemove_Ixs = []
  191. for subjIx in range(dataHere_Trimmed_NonRedundant.shape[0]):
  192. subjID = subjKeysAll[subjIx]
  193. subjVec = dataHere_Trimmed_NonRedundant[dataHere_Trimmed_NonRedundant['subjectkey']==subjID][keysHere[1:]].values[0,:].copy()
  194. numNaNs = np.where(np.isnan(subjVec))[0].shape[0]
  195. percNaNs = np.floor(((numNaNs/subjVec.shape[0])*100))
  196. if percNaNs > percentThresh_Subjs:
  197. subjsToRemove_IDs.append(subjID)
  198. subjsToRemove_Ixs.append(subjIx)
  199. print(f"Flagging subject {subjID} for having {percNaNs}% missing behavioral/clinical/cognitive measures...")
  201. subjsToRemove_IDs = np.asarray(subjsToRemove_IDs)
  202. subjsToRemove_Ixs = np.asarray(subjsToRemove_Ixs)
  204. dataHere_TrimFinal = dataHere_TrimFinal.drop(subjsToRemove_Ixs,axis=0)
  206. print(f"\nTrimming the dataset by {subjsToRemove_Ixs.shape[0]} subjects that were missing > {percentThresh_Subjs}% of behavioral/clinical/cognitive measures. "+
  207. f"New number of subjects = {dataHere_TrimFinal.shape[0]}.\n")
  209. # QA:
  210. subjIxVec_All = np.arange(dataHere_Trimmed_NonRedundant.shape[0])
  211. subjIxVec_All = np.delete(subjIxVec_All,subjsToRemove_Ixs)
  212. if not np.all(subjKeysAll[subjIxVec_All] == dataHere_TrimFinal['subjectkey'].to_numpy()):
  213. print(f"WARNING: subject ID order is off, please check.")
  215. elif not removeSparseSubjs:
  216. print(f"\nNOT trimming the dataset by subjects missing behavioral/clinical/cognitive measures. Number of subjects = {dataHere_TrimFinal.shape[0]}.\n")
  218. # SAVE
  219. if saveResults:
  220. dataHere_TrimFinal.to_csv(saveDir + 'dataHere_TrimFinal_BeforeGaussTest.csv')
  222. return dataHere_TrimFinal
  224. ####################################################################################
  225. def fingerprints_transform_Cocuzza(dataHere_TrimFinal,normData=False,useZ=False,saveResults=False,saveDir=None):
  226. '''
  227. The goal of this step is to normalize variables that are deemed non-Gaussian. NOTE that this was aided by an RStudio toolkit called bestNormalize, with the following workflow:
  229. 1) Run fingerprints_curate_Cocuzza() [function above] to curate which variables will be part of the final set for analyses.
  230. 2) Run those variables through this function fingerprints_transform_Cocuzza() to identify which need to be normalized --> save the identified non-Gaussian variables.
  231. 3) Load non-Gaussian variables identified in step 2 into RStudio and run bestNormalize (see cnd_project_imputation_R_pipe.Rmd).
  232. 4) Re-run this function fingerprints_transform_Cocuzza() to see if any of the variables are still being flagged as non-Gaussian (likely a small proportion of original set of variables)
  233. 5) The variables that are still non-Gaussian will need to be manually inspected. The vast majority were 0-inflated (or otherwise point-inflated), in which case we binarized (see next function, fingerprints_binarize_Cocuzza()), but a small number were reasonably close to having a normal distribution and were kept as-is. It may also be reasonable to trim those last remaining variables from further analyses, depending on how problematic the distribution is for later steps (e.g., strongly multimodal distributions might need further consideration)
  235. INPUTS:
  236. dataHere_TrimFinal: output dataframe of fingerprints_curate_Cocuzza(); NOTE if saveResults=True in fingerprints_curate_Cocuzza(), you could load this csv then use as input (or some variant thereof)
  237. normData: Boolean w/ Default False. Whether or not to pre-normalize with z-score or min-max.
  238. useZ: Boolean w/ Default False. Only when normData=True. If True, pre-normalize with z-score; if False, pre-normalize with min-max method.
  239. saveResults: Boolean w/ Default False (but we recommend setting to True); whether or not to also save out the curated dataframe (into .csv)
  240. saveDir: if saveResults=True; a directory path (include entire string) to save out results (csv) from transformation steps
  242. OUTPUTS:
  243. dataHere_CND_ToTransform: a dataframe of the variables you should pass to R for transformation; if you do a 2nd pass of this function (see notes above), the resulting dataframe from R will be the input dataHere_TrimFinal (make sure normData is False)
  244. dataHere_TrimFinal_Normed: only if normData=True; an intermedite output dataframe; all the curated data before identifying which variables need to go to R for transforming; i.e., all data in the input dataHere_TrimFinal plus min-max normalization or z-scoring (based on useZ input)
  246. '''
  248. # Normalize data: optional
  249. if normData:
  250. dataHere_TrimFinal_Normed = dataHere_TrimFinal.copy()
  252. if useZ:
  253. # NOTE: indexed at 1 to skip the first column (subject IDs), but may need to adjust to 2 if another, extra index column was added in the last step
  254. for keyIx in range(1,dataHere_TrimFinal.keys().to_numpy().shape[0]):
  255. keyHere = dataHere_TrimFinal.keys().to_numpy()[keyIx]
  256. vecHere = dataHere_TrimFinal[keyHere].to_numpy().copy()
  257. zVec = stats.zscore(vecHere,nan_policy='omit').copy()
  258. dataHere_TrimFinal_Normed[keyHere] = zVec.copy()
  260. # SAVE:
  261. if saveResults:
  262. dataHere_TrimFinal_Normed.to_csv(saveDir + 'dataHere_TrimFinal_Normed_Z.csv')
  264. elif not useZ:
  265. # NOTE: indexed at 1 to skip the first column (subject IDs), but may need to adjust to 2 if another, extra index column was added in the last step
  266. for keyIx in range(1,dataHere_TrimFinal.keys().to_numpy().shape[0]):
  267. keyHere = dataHere_TrimFinal.keys().to_numpy()[keyIx]
  268. vecHere = dataHere_TrimFinal[keyHere].to_numpy().copy()
  269. normVec = minMaxNorm(vecHere,minVal=0,maxVal=1).copy()
  270. dataHere_TrimFinal_Normed[keyHere] = normVec.copy()
  272. # SAVE:
  273. if saveResults:
  274. dataHere_TrimFinal_Normed.to_csv(saveDir + 'dataHere_TrimFinal_Normed_MinMax.csv')
  276. elif not normData:
  277. # SAVE:
  278. if saveResults:
  279. dataHere_TrimFinal.to_csv(saveDir + 'dataHere_TrimFinal_NotNormed.csv')
  282. # Identify non-gaussian data for transformation
  283. # Here, I'm running some tests to identify which variables need transforming, I'll then save those out --> go to R --> load them back in
  284. # See cnd_project_imputation_R_pipe.rmd
  285. if normData:
  286. keysHere_All = dataHere_TrimFinal_Normed.keys().to_numpy().copy()
  287. dataHere_TestGauss_Pass1 = dataHere_TrimFinal_Normed.copy()
  289. elif not normData:
  290. keysHere_All = dataHere_TrimFinal.keys().to_numpy().copy()
  291. dataHere_TestGauss_Pass1 = dataHere_TrimFinal.copy()
  293. # Test whether scale's distribution differs from normal distribution (D’Agostino’s K^2 test)
  294. keyIxs_NonNormal = []
  295. keyIxs_NonNormal_Cusp = [] # i.e., metrics that only just fail the test, so inspect further
  296. for keyIx in range(1,keysHere_All.shape[0]):
  297. keyHere = keysHere_All[keyIx]
  298. dataVecHere = dataHere_TestGauss_Pass1[keyHere].to_numpy().copy()
  299. kStat,pVal = stats.normaltest(dataVecHere,nan_policy='omit')
  300. if pVal<0.05:
  301. #print(f"{keyHere} (py index {keyIx}) is likely not normal")
  302. keyIxs_NonNormal.append(keyIx)
  303. if pVal > 0.01:
  304. print(f"{keyHere} (py index {keyIx}) may be on the cusp, please check")
  305. keyIxs_NonNormal_Cusp.append(keyIx)
  307. keyIxs_NonNormal = np.asarray(keyIxs_NonNormal)
  308. keyIxs_NonNormal_Cusp = np.asarray(keyIxs_NonNormal_Cusp)
  309. print(f"\nOf {keysHere_All.shape[0]-1} measures, {keyIxs_NonNormal.shape[0]} were identified as not matching normal distribution\n")
  311. # Test Gaussianity with Shapiro-Wilk test
  312. keyIxs_ShapiroWilk = []
  313. keyIxs_ShapiroWilk_Cusp = [] # i.e., metrics that only just fail the test, so inspect further
  314. for keyIx in range(1,keysHere_All.shape[0]):
  315. keyHere = keysHere_All[keyIx]
  316. dataVecHere = dataHere_TestGauss_Pass1[keyHere].to_numpy().copy()
  317. nanIxs = np.where(np.isnan(dataVecHere))[0]
  318. if nanIxs.shape[0]!=0:
  319. dataVecHere = np.delete(dataVecHere,nanIxs)
  320. kStat,pVal = stats.shapiro(dataVecHere)
  321. if pVal<0.05:
  322. #print(f"{keyHere} (py index {keyIx}) is likely not normal")
  323. keyIxs_ShapiroWilk.append(keyIx)
  324. if pVal > 0.01:
  325. print(f"{keyHere} (py index {keyIx}) may be on the cusp, please check")
  326. keyIxs_ShapiroWilk_Cusp.append(keyIx)
  328. keyIxs_ShapiroWilk = np.asarray(keyIxs_ShapiroWilk)
  329. keyIxs_ShapiroWilk_Cusp = np.asarray(keyIxs_ShapiroWilk_Cusp)
  330. print(f"\nOf {keysHere_All.shape[0]-1} measures, {keyIxs_ShapiroWilk.shape[0]} were identified as non-Gaussian with Shapiro-Wilk\n")
  332. # Test normality with Anderson-Darling test
  333. keyIxs_AD = []
  334. keyIxs_AD_Cusp = [] # i.e., metrics that only just fail the test, so inspect further
  335. for keyIx in range(1,keysHere_All.shape[0]):
  336. keyHere = keysHere_All[keyIx]
  337. dataVecHere = dataHere_TestGauss_Pass1[keyHere].to_numpy().copy()
  338. nanIxs = np.where(np.isnan(dataVecHere))[0]
  339. if nanIxs.shape[0]!=0:
  340. dataVecHere = np.delete(dataVecHere,nanIxs)
  341. testModel = stats.anderson(dataVecHere,dist='norm')
  342. if testModel.statistic>testModel.critical_values[2]: # see testModel.significance_level for associated alphas for each critical value:
  343. keyIxs_AD.append(keyIx)
  344. if testModel.statistic<testModel.critical_values[4]: # see testModel.significance_level for associated alphas for each critical value
  345. print(f"{keyHere} (py index {keyIx}) may be on the cusp, please check")
  346. keyIxs_AD_Cusp.append(keyIx)
  348. keyIxs_AD = np.asarray(keyIxs_AD)
  349. keyIxs_AD_Cusp = np.asarray(keyIxs_AD_Cusp)
  350. print(f"\nOf {keysHere_All.shape[0]-1} measures, {keyIxs_AD.shape[0]} were identified as non-normal with Anderson-Darling\n")
  352. # Test normality with Kolmogorov-Smirnov test
  353. keyIxs_KS = []
  354. keyIxs_KS_Cusp = [] # i.e., metrics that only just fail the test, so inspect further
  355. for keyIx in range(1,keysHere_All.shape[0]):
  356. keyHere = keysHere_All[keyIx]
  357. dataVecHere = dataHere_TestGauss_Pass1[keyHere].to_numpy().copy()
  358. nanIxs = np.where(np.isnan(dataVecHere))[0]
  359. if nanIxs.shape[0]!=0:
  360. dataVecHere = np.delete(dataVecHere,nanIxs)
  361. kStat,pVal = stats.kstest(dataVecHere,'norm')
  362. if pVal<0.05:
  363. #print(f"{keyHere} (py index {keyIx}) is likely not normal")
  364. keyIxs_KS.append(keyIx)
  365. if pVal > 0.01:
  366. print(f"{keyHere} (py index {keyIx}) may be on the cusp, please check")
  367. keyIxs_KS_Cusp.append(keyIx)
  369. keyIxs_KS = np.asarray(keyIxs_KS)
  370. keyIxs_KS_Cusp = np.asarray(keyIxs_KS_Cusp)
  371. print(f"\nOf {keysHere_All.shape[0]-1} measures, {keyIxs_KS.shape[0]} were identified as non-normal with Kolmogorov-Smirnov\n")
  373. # Identify which measures/variables were identified across a few tests for normality/Gaussianity
  374. uniqueVec,countsVec = np.unique(np.concatenate((keyIxs_NonNormal,keyIxs_ShapiroWilk,keyIxs_AD,keyIxs_KS)),return_counts=True)
  375. nonNormalWinnerIxs = np.where(countsVec>=2)[0].copy()
  376. nonNormalWinners = uniqueVec[nonNormalWinnerIxs]
  377. #print(f"\nAcross all 4 tests, {nonNormalWinners.shape[0]} measures were identified as "+
  378. # f"non-normal for >=2 tests: \n{keysHere_All[nonNormalWinners]} \n(py indices: {nonNormalWinners})\n")
  379. print(f"\nAcross all 4 tests, {nonNormalWinners.shape[0]} measures were identified as non-normal for >=2 tests:\n")
  381. # And identify which metrics were on the cusp across a few tests
  382. uniqueVec,countsVec = np.unique(np.concatenate((keyIxs_NonNormal_Cusp,keyIxs_ShapiroWilk_Cusp,keyIxs_AD_Cusp,keyIxs_KS_Cusp)),return_counts=True)
  383. cuspWinnerIxs = np.where(countsVec>=3)[0].copy() # lower number = less strict; higher number = more strict; max = 4 (b/c that's the number of tests I used)
  384. if cuspWinnerIxs.shape[0]!=0:
  385. cuspWinners = uniqueVec[cuspWinnerIxs]
  386. print(f"\nAcross all 4 tests, the following {cuspWinners.shape[0]} measures were on the cusp "+
  387. f"for being non-normal for >=3 tests: \n{keysHere_All[cuspWinners]} \n(py indices: {cuspWinners})\n")
  389. cuspWinnerIxs_Adj = np.zeros_like(cuspWinnerIxs)
  390. for ixHere in range(cuspWinnerIxs.shape[0]):
  391. cuspWinnerIxs_Adj[ixHere] = np.where(nonNormalWinners==cuspWinners[ixHere])[0]
  393. # And generate the final list to send to R
  394. nonNormalWinners = np.delete(nonNormalWinners,cuspWinnerIxs_Adj)
  395. keysToTransform = np.concatenate((np.asarray(['subjectkey']),keysHere_All[nonNormalWinners]))
  397. # Curate dataframe and save
  398. dataHere_CND_ToTransform = dataHere_TestGauss_Pass1.copy()
  399. dataHere_CND_ToTransform = dataHere_CND_ToTransform[keysToTransform]
  401. if saveResults:
  402. if normData:
  403. if useZ:
  404. saveFile = saveDir + 'dataHere_TrimFinal_Normed_Z_ToTransform.csv'
  405. elif not useZ:
  406. saveFile = saveDir + 'dataHere_TrimFinal_Normed_MinMax_ToTransform.csv'
  407. elif not normData:
  408. saveFile = saveDir + 'dataHere_TrimFinal_NotNormed_ToTransform.csv'
  410. dataHere_CND_ToTransform.to_csv(saveFile)
  412. if normData:
  413. return dataHere_CND_ToTransform, dataHere_TrimFinal_Normed
  414. elif not normData:
  415. return dataHere_CND_ToTransform
  417. ####################################################################################
  418. def fingerprints_binarize_Cocuzza(dataHere_ToBinarize,normData=False,useZ=False,saveResults=False,saveDir=None):
  419. '''
  420. The goal of this step is to binarize measures that have been identified as 0-inflated (after 2 runs through fingerprints_transform_Cocuzza() above and RStudio's bestNormalize).
  421. NOTE that there are a few possible ways to deal with point-inflated data (e.g., poisson transforms; gamma models; etc.), but we opted for binarization given that it's straightforward.
  423. INPUTS:
  424. dataHere_ToBinarize: dataframe; where general formatting is the same as other functions, but the variables (columns) here are those identified for binarization
  425. normData: whether or not data was normed in prior transformation step/function; note this is only relevant for strings used to save results if saveResults=True
  426. useZ: if normData=True, whether or not data was normalized with z-scoring in prior transformation step/function; note this is only relevant for strings used to save results if saveResults=True
  427. saveResults: Boolean w/ Default False (but we recommend setting to True); whether or not to also save out the results dataframe (into .csv)
  428. saveDir: if saveResults=True; a directory path (include entire string) to save out results (csv) from binarization steps
  430. OUTPUTS:
  431. dataHere_Transformed_Binarized: dataframe with measures/variables binarized
  432. '''
  434. # Binarize highly non-gaussians (0-inflated etc.): binarize based on mode
  436. # Initialize dictionary to make new dataframe
  437. dataHere_Transformed_Pass2_DICT = {'subjectkey':dataHere_ToBinarize['subjectkey'].to_numpy()}
  439. for metricIx in range(1,keysHere_All.shape[0]):
  440. thisKey = keysHere_All[metricIx]
  442. vecHere = dataHere_ToBinarize[thisKey].to_numpy()
  443. modeVal,countVal = stats.mode(vecHere)
  445. nanIxs = np.where(np.isnan(vecHere))[0] # hold out nans temporarily
  446. vecHere[nanIxs] = 0
  448. modeIxs = np.where(np.round(vecHere,4)==np.round(modeVal,4))[0]
  449. vecHere_Bin = np.zeros((vecHere.shape[0]))
  450. vecHere_Bin[modeIxs] = 1
  452. vecHere_Bin[nanIxs] = np.nan # add nans back in
  454. # NOTE: can use sklearn (see line below), but I think it does the opposite mapping from what we want (1's given to the non-mode)
  455. # vecHere_Bin = sklearn.preprocessing.binarize(vecHere.reshape(-1,1),threshold=modeVal)[:,0].copy()
  457. dataHere_Transformed_Pass2_DICT[thisKey] = vecHere_Bin.copy()
  459. dataHere_Transformed_Binarized = pd.DataFrame(dataHere_Transformed_Pass2_DICT).copy()
  461. # SAVE
  462. if saveResults:
  463. if normData:
  464. if useZ:
  465. dataHere_Transformed_Binarized.to_csv(saveDir + 'dataHere_TrimFinal_Normed_Z_Transformed_Pass2.csv')
  466. elif not useZ:
  467. dataHere_Transformed_Binarized.to_csv(saveDir + 'dataHere_TrimFinal_Normed_MinMax_Transformed_Pass2.csv')
  468. elif not normData:
  469. dataHere_Transformed_Binarized.to_csv(saveDir + 'dataHere_TrimFinal_NotNormed_Transformed_Pass2.csv')
  471. return dataHere_Transformed_Binarized
  473. ####################################################################################
  474. def fingerprints_imputation_Cocuzza(dataHere_Collated_ToImpute,percentTest=20,reNormData=True,useZ=False,saveResults=False,saveDir=None):
  476. '''
  477. The goal of this step is to impute missing data points (example: a given participant did not complete a survey) using the MissForest algorithm. See Chopra et al. 2025, Sci Data, for extensive testing on the TCP dataset, which suggested that missforest yielded the most robust imputation results.
  478. NOTE: in the manuscript, we used a conservative train/test/validation split of the data. It is up to the user to divvy up their data appropriately.
  480. INPUTS:
  481. dataHere_Collated_ToImpute: this is a dataframe where all of the results of the prior steps (curating, transforming (twice), and binarizing) are all collated back together (i.e., the columns are stacked back into 1 dataframe)
  482. percentTest: percent of the data to set as the test set. Default is 20. e.g., if 20 is given, this is an 80/20 train/test split.
  483. reNormData: whether or not to re-normalize the data; the bestNormalize transformation process in R appears to put some (but not all) variables in a different space/range; and this is particularly important for min-max norming (i.e., putting data back into range of 0-1); Also should consider using this when prior steps/functions had normData=False. The idea here is that the prior transformations were on raw scores, then here normalize just before imputation.
  484. useZ: Boolean w/ Default False. Only when normData=True. If True, pre-normalize with z-score; if False, pre-normalize with min-max method.
  485. saveResults: Boolean w/ Default False (but we recommend setting to True); whether or not to also save out the curated dataframe (into .csv)
  486. saveDir: if saveResults=True; a directory path (include entire string) to save out results (csv) from transformation steps
  488. OUTPUTS:
  489. dataHere_Imputed: dataframe where formatting matches inputted dataHere_Collated_ToImpute, but missing values have been imputed
  490. '''
  492. # Note: dataHere_Collated_ToImpute was pre-computed (and saved, if saveResults=True) above; user needs to collate those dataframes before using this function
  493. if reNormData:
  494. dataHere_ToImpute_Final = dataHere_Collated_ToImpute.copy()
  495. if not useZ:
  496. for keyIx in range(1,dataHere_Collated_ToImpute.keys().to_numpy().shape[0]):
  497. keyStrHere = dataHere_Collated_ToImpute.keys().to_numpy()[keyIx]
  498. vecHere = dataHere_Collated_ToImpute[keyStrHere].to_numpy().copy()
  499. vecHere_Normed = minMaxNorm(vecHere).copy()
  500. dataHere_ToImpute_Final[keyStrHere] = vecHere_Normed.copy()
  502. elif not reNormData:
  503. dataHere_ToImpute_Final = dataHere_Collated_ToImpute.copy()
  505. # MissForest with sklearn
  507. subjIDs_DF = dataHere_ToImpute_Final['subjectkey'].to_numpy().copy()
  508. dataHere_ToImpute_Final_TrainTest = dataHere_ToImpute_Final.copy()
  509. nSubjsCND = dataHere_ToImpute_Final.shape[0]
  510. nSubjsTrainTest = dataHere_ToImpute_Final_TrainTest.shape[0]
  511. percTrainTest = (nSubjsTrainTest * 100)/nSubjsCND
  513. print(f"Final dataset for imputation: {nSubjsTrainTest} subjects (of N={nSubjsCND} total) by {dataHere_ToImpute_Final_TrainTest.shape[1]} measures. "+
  514. f"{percTrainTest}% subjects (n = {nSubjsTrainTest}) will be used for train/test")
  516. percentTrain = 100 - percentTest
  517. nSubjs_ImputePipe = dataHere_ToImpute_Final_TrainTest.shape[0]
  518. numTrain = int(np.round((percentTrain * nSubjs_ImputePipe)/100))
  519. numTest = int(np.round((percentTest * nSubjs_ImputePipe)/100))
  520. if numTest + numTrain != nSubjs_ImputePipe:
  521. print(f"train and test subject numbers do not equal full N, please check.")
  523. print(f"\nFor missForest: training set: {percentTrain}% of subjects (n = {numTrain}). Test set: {percentTest}% of subjects (n = {numTest}). (Total N={nSubjs_ImputePipe}).")
  525. # Loop over scales; NOTE: can loop over shuffles and take average (for continuous data) or mode (categorical data)
  526. subjVec_Shuffled = np.arange(nSubjs_ImputePipe)
  527. np.random.shuffle(subjVec_Shuffled)
  528. subjsTrain = subjVec_Shuffled[:numTrain].copy()
  529. subjsTest = subjVec_Shuffled[numTrain:].copy()
  531. keysHere_All = dataHere_ToImpute_Final_TrainTest.keys().to_numpy()[1:].copy()
  533. dataHere_Imputed_DICT = {'subjectkey':dataHere_ToImpute_Final_TrainTest['subjectkey'].to_numpy()}
  534. for keyIx in range(keysHere_All.shape[0]):
  535. keyStr = keysHere_All[keyIx]
  536. allOtherKeys = keysHere_All.copy()
  537. allOtherKeys = np.delete(allOtherKeys,keyIx)
  539. yData_Orig = pd.DataFrame({keyStr:dataHere_ToImpute_Final_TrainTest[keyStr].to_numpy()})[keyStr].copy()
  540. nanIxs = np.where(np.isnan(yData_Orig))[0]
  542. xData = dataHere_ToImpute_Final_TrainTest[allOtherKeys].values.copy()
  544. if nanIxs.shape[0]!=0:
  545. #print(f"Predicting {keyStr}...")
  547. if keyStr in keysHere_Trans2_Adj: # Categorical (binarized)
  548. yData = mode_imputation(yData_Orig).to_numpy().copy() # for categorical data
  550. xData_Train = xData[subjsTrain,:].copy()
  551. xData_Test = xData[subjsTest,:].copy()
  553. yData_Train = yData[subjsTrain].copy()
  554. yData_Test = yData[subjsTest].copy()
  556. modelHere = RandomForestClassifier()
  557. modelHere.fit(xData,yData)
  559. elif keyStr not in keysHere_Trans2_Adj: # Continuous
  560. #yData = mean_median_imputation(yData_Orig,kind='mean').to_numpy().copy() # for continuous data
  561. yData = mean_median_imputation(yData_Orig,kind='median').to_numpy().copy() # for continuous data
  563. xData_Train = xData[subjsTrain,:].copy()
  564. xData_Test = xData[subjsTest,:].copy()
  566. yData_Train = yData[subjsTrain].copy()
  567. yData_Test = yData[subjsTest].copy()
  569. modelHere = RandomForestRegressor()
  570. modelHere.fit(xData,yData)
  572. yPred = modelHere.predict(xData_Test[:nanIxs.shape[0],:]) # values
  573. yData_Imputed = yData_Orig.copy()
  574. yData_Imputed[nanIxs] = yPred
  576. else:
  577. print(f"{keyStr} is not missing data, skipping...")
  578. yData_Imputed = yData_Orig.copy()
  580. dataHere_Imputed_DICT[keyStr] = yData_Imputed.copy()
  582. dataHere_Imputed = pd.DataFrame(dataHere_Imputed_DICT)
  584. nanIxsDF = dataHere_Imputed.isnull().sum()[dataHere_Imputed.isnull().sum() != 0].index # use ==0 to get not-nans
  585. if nanIxsDF.shape[0] != 0:
  586. print(f"There is still missing data points, please check.")
  588. if saveResults:
  589. if reNormData:
  590. if useZ:
  591. dataHere_Imputed.to_csv(saveDir + 'dataHere_Normed_Z_Transformed_Imputed.csv')
  592. elif not useZ:
  593. dataHere_Imputed.to_csv(saveDir + 'dataHere_Normed_MinMax_Transformed_Imputed.csv')
  594. elif not reNormData:
  595. dataHere_Imputed.to_csv(saveDir + 'dataHere_NotNormed_Imputed.csv')
  597. return dataHere_Imputed
  599. ####################################################################################
  600. def fingerprints_distance_Cocuzza(dataHere_Imputed,startIx_CSV=2):
  601. '''
  602. The goal of this step is to estimate individual differences correlation. Here, each variable-variable pair (i.e., each pair of behavioral measures) is correlated across subjects. So the resulting correlation score indicates the extent that individual differences (across those 2 given measures) are similar. NOTE: we used spearman rho because of it's inference for rank correlation that doesn't assume specific distributions a priori (i.e., is less parametric than r)
  604. INPUTS:
  605. dataHere_Imputed: output dataframe from prior step/function (imputation step) above.
  606. startIx_CSV: an artifact of pandas is that sometimes an indexing column is added. You always want this to be at least 1 (because we want to skip the true 1st colum which is subject IDs), but maybe increase to 2 (default) or even 3 etc. if more columns were added to left of that subjectkey column.
  608. OUTPUTS:
  609. behavCorrs_IndivDiffs: an adjacency matrix (square, symmetric). x/y dimensions correspond to behavioral variables in input dataframe.
  610. '''
  612. # Correlation matrix: individual differences across subjects
  613. numVarsCND = dataHere_Imputed.shape[1]-startIx_CSV
  614. behavCorrs_IndivDiffs = np.zeros((numVarsCND,numVarsCND))
  615. for metricIx in range(numVarsCND):
  616. thisTargetKey = dataHere_Imputed.keys().to_numpy()[startIx_CSV:][metricIx]
  618. for metricIx_Next in range(numVarsCND):
  619. thisSourceKey = dataHere_Imputed.keys().to_numpy()[startIx_CSV:][metricIx_Next]
  621. r,p = stats.spearmanr(dataHere_Imputed[thisTargetKey].to_numpy(),dataHere_Imputed[thisSourceKey].to_numpy(),nan_policy='omit')
  622. behavCorrs_IndivDiffs[metricIx,metricIx_Next] = r
  624. np.fill_diagonal(behavCorrs_IndivDiffs,0)
  626. return behavCorrs_IndivDiffs
  628. ####################################################################################
  629. def fingerprints_clustering_Cocuzza(behavCorrs_IndivDiffs,dataHere_Imputed,nClusters=4,startIx_CSV=2):
  630. '''
  631. The goal of this step is to perform heirarchical agglomerative clustering on individual differences correlations. See cnd_project_clustering_checks_pipe.Rmd for example usage of. Rpackage that helps optimize number of clusters to use.
  633. INPUTS:
  634. behavCorrs_IndivDiffs: output matrix from step/function above (distance step).
  635. dataHere_Imputed: output dataframe from step/function above (imputation step).
  636. nClusters: Number of clusters to use; Default = 4; based on R package nbClust.
  637. startIx_CSV: an artifact of pandas is that sometimes an indexing column is added. You always want this to be at least 1 (because we want to skip the true 1st colum which is subject IDs), but maybe increase to 2 (default) or even 3 etc. if more columns were added to left of that subjectkey column.
  639. OUTPUTS:
  640. clusterBoundaries: array of size nClusters x 3: column 1 = start ix, column 2 = stop ix (py --> +1), column 3 = number of features in cluster.
  641. clusterOrder: variable that assigns each variable from the input dataframe to a cluster.
  642. '''
  644. originalOrder = dataHere_Imputed.keys().to_numpy()[startIx_CSV:].copy()
  645. scoreNames = dataHere_Imputed.keys().to_numpy()[startIx_CSV:].copy()
  647. linkagMatrix = scipy.cluster.hierarchy.ward(behavCorrs_IndivDiffs) # same as colLink
  648. cutTree = scipy.cluster.hierarchy.cut_tree(linkagMatrix, n_clusters = nClusters).flatten() # note: n_clusters can be 2D, will return 2D group membership
  649. percentileList = pd.DataFrame({'column_names':scoreNames,'cluster_membership':cutTree})
  650. clusterList = percentileList.sort_values(by='cluster_membership')
  652. clusterBoundaries = np.zeros((nClusters,3)) # Col 1 = start ix, Col 2 = stop ix (py --> +1), Col 3 = number of features in cluster
  653. clusterOrder = np.zeros((originalOrder.shape[0])) # Use this as a sorting/indexing variable
  654. for clusterIx in range(nClusters):
  655. if not useHandOrder:
  656. clusterSetIxs = np.where(clusterList['cluster_membership'].to_numpy() == clusterIx)[0]
  657. elif useHandOrder:
  658. clusterSetIxs = np.where(clusterListAdj['cluster_membership'].to_numpy() == clusterIx)[0] # NOTE: hand-ordered above
  659. clusterBoundaries[clusterIx,0] = clusterSetIxs[0]
  660. clusterBoundaries[clusterIx,1] = clusterSetIxs[-1] + 1
  661. clusterBoundaries[clusterIx,2] = clusterSetIxs.shape[0]
  663. if not useHandOrder:
  664. variablesThisCluster = clusterList['column_names'].to_numpy().copy()
  665. elif useHandOrder:
  666. variablesThisCluster = clusterListAdj['column_names'].to_numpy().copy() # NOTE: hand-ordered above
  668. print(f"\nCluster {clusterIx+1}:\n{clusterList[clusterList['cluster_membership']==clusterIx]['column_names'].to_numpy()}")
  670. for varIx in range(variablesThisCluster.shape[0]):
  671. thisVar = variablesThisCluster[varIx]
  672. origIx = np.where(originalOrder==thisVar)[0]
  673. if not useHandOrder:
  674. newIx = np.where(clusterList['column_names'].to_numpy()==thisVar)[0]
  675. elif useHandOrder:
  676. newIx = np.where(clusterListAdj['column_names'].to_numpy()==thisVar)[0] # NOTE: hand-ordered above
  677. #clusterOrder[varIx] = origIx
  678. clusterOrder[newIx] = origIx
  680. # NOTE: can verify with:
  681. # np.array_equal(originalOrder[clusterOrder.astype(int)],clusterList['column_names'].to_numpy())
  683. # NOTE: an alternative with sklearn; results are effectively equivalent to scipy method above:
  684. #numClustersStr = 4
  685. # Appending distance between clusters at every iteration in for loop
  686. #hClustering_Distances_All = []
  687. #hClustering_CH_Index_All = []
  688. #for numClusters in range(2,numVarsCND):
  689. # hClustering = sklearn.cluster.AgglomerativeClustering(n_clusters=numClusters, metric='euclidean', linkage='ward', compute_distances=True).fit(behavCorrs_IndivDiffs)
  690. # hClustering_Distances_All.append(hClustering.distances_)
  691. # hClustering_CH_Index_All.append(sklearn.metrics.calinski_harabasz_score(behavCorrs_IndivDiffs,hClustering.labels_))
  692. #hClustering_CH_Index_All = np.asarray(hClustering_CH_Index_All)
  693. #hClustering_Distances = hClustering_Distances_All[0].copy() # note: all sub-arrays appear to be equal? why do loop? this appears to be equal to colLink
  694. #xNumClusters = np.array([indexHere for indexHere in range(numVarsCND,1,-1)]) # number of clusters in descending oder
  696. return clusterBoundaries, clusterOrder
  698. ####################################################################################
  699. def fingerprints_PCA_Cocuzza(behavCorrs_IndivDiffs,dataHere_Imputed,numClusters=4, startIx_CSV=2,printTestName=True,numComps=1,useAbs=Fale):
  700. '''
  701. After clustering has been performed (see above), run PCA on the original data (i.e.., not the distance matrix), one cluster at a time. In the manuscript we focused on the 1st PC for each cluster.
  703. INPUTS:
  704. behavCorrs_IndivDiffs: output matrix from step/function above (distance step).
  705. dataHere_Imputed: output dataframe from step/function above (imputation step).
  706. numClusters: Number of clusters to use; Default = 4; based on R package nbClust.
  707. startIx_CSV: an artifact of pandas is that sometimes an indexing column is added. You always want this to be at least 1 (because we want to skip the true 1st colum which is subject IDs), but maybe increase to 2 (default) or even 3 etc. if more columns were added to left of that subjectkey column.
  708. printTestName: Default True; whether to include test names in the print
  709. numComps: Default 1; For PCA, either use a number (less than # features) or None
  710. useAbs: Default False; Use absolute values of PCA scores for visualizing
  713. OUTPUTS:
  714. pcaScores_All: array of size number of subjects x number of clusters. This gives log-likelihood of each subject to express each cluster.
  715. loadings_All: dictionary where keys = cluster IDs; in each, an array of factor loadings corresponding to measures in that cluster
  717. '''
  719. modelHere = sklearn.cluster.AgglomerativeClustering(n_clusters=numClusters, metric='euclidean', linkage='ward', compute_distances=True)
  720. modelHere.fit_predict(behavCorrs_IndivDiffs)
  721. labelsHere = modelHere.labels_
  723. pcaScores_All = np.zeros((dataHere_Imputed.shape[0],numClusters))
  724. loadings_All = {} # keys given by cluster, e.g., 'cluster_1', etc.
  725. for labelIx in range(numClusters):
  726. tempList = []
  727. clusterIxsHere = np.where(labelsHere==labelIx)[0]
  728. testNamesHere = dataHere_Imputed.keys().to_numpy()[startIx_CSV:][clusterIxsHere]
  729. dataHere_Cluster = dataHere_Imputed[testNamesHere].copy()
  730. if numComps is not None:
  731. modelHere = sklearn.decomposition.PCA(n_components=numComps)
  732. else:
  733. modelHere = sklearn.decomposition.PCA()
  734. modelHere.fit(dataHere_Cluster)
  735. modelHere_Transform = modelHere.fit_transform(dataHere_Cluster).copy() # samples (subjects) x components
  737. # NOTE: other variables are here that aren't returned; feel free to edit to save out/return for other research studies
  738. compsHere = modelHere.components_.copy() # components x features (psych measures)
  739. explVarHere = modelHere.explained_variance_.copy() # components
  740. explVarRatioHere = modelHere.explained_variance_ratio_.copy() # components
  741. inverseTransformHere = modelHere.inverse_transform(modelHere_Transform).copy() # samples (subjects) x features (psych measures) -> original space
  742. pcaScoreHere = modelHere.score_samples(dataHere_Cluster).copy() # samples (subjects) -> NOTE: this is the log-likelihood of each sample (subject)
  743. loadingsHere = compsHere.T * np.sqrt(explVarHere) # features (psych measures) x components
  745. print(f"Cluster {labelIx+1}: the 1st PC explains {np.round(explVarRatioHere[0]*100,2)}% of the variance (of {np.round(np.nansum(explVarRatioHere)*100,2)}% total).")
  747. pcaScores_All[:,labelIx] = pcaScoreHere.copy()
  748. loadings_All['cluster_'+str(labelIx+1)] = loadingsHere.copy()
  750. return pcaScores_All, loadings_All
  752. ####################################################################################
  753. # Helper function to perform min-max normalization. Also known as min-max scaling and feature scaling
  754. def minMaxNorm(data,minVal=0,maxVal=1):
  755. normedData = minVal + (((data-np.nanmin(data))*(maxVal-minVal))/(np.nanmax(data)-np.nanmin(data)))
  756. return normedData
  758. ####################################################################################
  759. # Helper function for imputation
  760. def add_label(data,attr,name_notnan = 'Training', name_nan = 'Predict'):
  761. null_pos = data[attr][data[attr].isnull()].index
  762. data['Label'] = name_notnan
  763. data['Label'].iloc[null_pos] = name_nan
  765. ####################################################################################
  766. # Helper function for imputation
  767. def mean_median_imputation(data, kind = 'mean'):
  768. if kind=='mean':
  769. return data.fillna(np.mean(data.dropna()))
  770. elif kind == 'median':
  771. return data.fillna(np.median(data.dropna()))
  773. ####################################################################################
  774. # Helper function for imputation
  775. def mode_imputation(data):
  776. return data.fillna(stat.mode(data))

Fingerprints_Cocuzza.py at commit 939f37d, under GPL-3.0 · at the source

Overview

Authors: Carrisa V Cocuzza1,2, Sidhant Chopra3,4, Ashlea Segal5, Loïc Labache1,2, Rowena Chin2, Kaley Joss1, Avram J Holmes1
ORCID iDs: Carrisa V Cocuzza, Loïc Labache, Rowena Chin, Avram J Holmes
  1. Department of Psychiatry, Brain Health Institute, Rutgers University, Piscataway, NJ USA
  2. Department of Psychology, Yale University, New Haven, CT USA
  3. Orygen, Melbourne, VIC Australia
  4. Center for Youth Mental Health, University of Melbourne, Melbourne, VIC Australia
  5. Department of Neuroscience, Wu Tsai Institute, Yale University, New Haven, CT USA
Institutions: Rutgers, The State University of New Jersey (United States); Yale University (United States); The University of Melbourne (Australia); Orygen (Australia); Brown University (United States)
Journal: Nature communications, volume 17, issue 1, article 6678
Dates: received 25 June 2025; accepted 23 June 2026; published online 21 July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-75585-6 · PMID 42481484 · PMCID PMC13388696 · OpenAlex W4410826654
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.05.23.655864
Status: code verified
Categories: human (organism), depression (population), schizophrenia / psychosis (population)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Graphs, fMRI & imaging
Keywords: Network models, Dynamical systems, Anxiety, Psychosis, Depression
MeSH: Brain*, Mental Disorders*, Nerve Net*, Adult, Brain Mapping, Case-Control Studies, Cognition, Connectome, Female, Humans, Magnetic Resonance Imaging, Male, Young Adult (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: U.S. Department of Health &amp; Human Services | NIH | National Institute of Mental Health (NIMH) (R01MH120080, R01MH123245); NIMH NIH HHS (RF1 MH123245, R01 MH120080)
Citations: cited by 2 papers (Europe PMC); 179 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.

Repositories

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

HolmesLab/ClinicalNetDynamics

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 939f37da13ad1e1c2e1eeb3d4671251cb325545f, 26 September 2026
Languages: Python (4), R (2)
Size: 11 files, 6 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 2 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (4 files), pandas (2 files), scikit-learn (2 files), tidyverse (2 files), h5py (1 file), igraph (1 file), Matplotlib (1 file), psych (1 file), SciPy (1 file), seaborn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
8 files
At the source: github.com/HolmesLab/ClinicalNetDynamics

Zenodo 20074876

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (4 files), pandas (2 files), scikit-learn (2 files), tidyverse (2 files), h5py (1 file), igraph (1 file), Matplotlib (1 file), psych (1 file), SciPy (1 file), seaborn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
8 files
At the source:

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/s41467-026-75585-6.

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 12 scripts, each with its path and the digest of its content;
  • 12 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 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:

Read it in the paper: doi.org/10.1038/s41467-026-75585-6.

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, 7 authors, 5 keywords, 13 MeSH terms, 2 funders, 168 references.

Cite

This paper

Cocuzza, C. V., Chopra, S., Segal, A., Labache, L., Chin, R., Joss, K., & Holmes, A. J. (2026). Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease. Nature communications, 17(1), 6678. https://doi.org/10.1038/s41467-026-75585-6

BibTeX

@article{cocuzza2026brain,
author = {Cocuzza, Carrisa V and Chopra, Sidhant and Segal, Ashlea and Labache, Loïc and Chin, Rowena and Joss, Kaley and Holmes, Avram J},
title = {{Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {6678},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75585-6},
url = {https://doi.org/10.1038/s41467-026-75585-6},
pmid = {42481484},
pmcid = {PMC13388696}
}

RIS

TY - JOUR
AU - Cocuzza, Carrisa V
AU - Chopra, Sidhant
AU - Segal, Ashlea
AU - Labache, Loïc
AU - Chin, Rowena
AU - Joss, Kaley
AU - Holmes, Avram J
TI - Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/21
VL - 17
IS - 1
SP - 6678
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75585-6
UR - https://doi.org/10.1038/s41467-026-75585-6
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75585-6",
"type": "article-journal",
"title": "Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease",
"container-title": "Nature communications",
"author": [
{
"family": "Cocuzza",
"given": "Carrisa V"
},
{
"family": "Chopra",
"given": "Sidhant"
},
{
"family": "Segal",
"given": "Ashlea"
},
{
"family": "Labache",
"given": "Loïc"
},
{
"family": "Chin",
"given": "Rowena"
},
{
"family": "Joss",
"given": "Kaley"
},
{
"family": "Holmes",
"given": "Avram J"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "6678",
"DOI": "10.1038/s41467-026-75585-6",
"PMID": "42481484",
"PMCID": "PMC13388696",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75585-6",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
21
]
]
}
}

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.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: h5py, seaborn, scikit-learn, 4 other tools, 13 references
[2] 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: h5py, seaborn, scikit-learn, 4 other tools, 10 references
[3] doi:10.1038/s41467-026-75745-8 [code]
A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks.
Journal: Nature communications
In common: h5py, seaborn, scikit-learn, 4 other tools, 8 references
[4] doi:10.1162/netn.a.570 [code]
Higher-order statistics for constructing centered edge functional connectivity.
Journal: Network neuroscience (Cambridge, Mass.)
In common: igraph, seaborn, scikit-learn, 4 other tools, 9 references
[5] doi:10.1038/s41467-026-72931-6 [code]
Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 8 references
[6] doi:10.1038/s41467-026-73153-6 [code]
Latent neural architecture organising shared aesthetic evaluations of visual artworks.
Journal: Nature communications
In common: h5py, seaborn, scikit-learn, 4 other tools, 8 references
[7] doi:10.1038/s41593-026-02359-0 [code]
The cross-site reproducibility of MRI morphometric phenotypes in psychiatric disorders.
Journal: Nature neuroscience
In common: scikit-learn, pandas, SciPy, 2 other tools, schizophrenia / psychosis, 7 references
[8] 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: h5py, seaborn, scikit-learn, 4 other tools, depression, 6 references
[9] doi:10.1038/s44220-026-00680-y [code]
The neuroimaging correlates of depression established across six large-scale population datasets.
Journal: Nature. Mental health
In common: seaborn, tidyverse, scikit-learn, 4 other tools, depression, 6 references
[10] doi:10.7554/elife.110294 [code]
Arousal modulates functional connectivity through structured and hemispherically asymmetric community architecture during wakefulness.
Journal: eLife
In common: seaborn, scikit-learn, pandas, 3 other tools, 7 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.

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.