OSCR

Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease.

Code ↔ Paper

14 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 14 matches · 4 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Clustering Analysis ↔ ClusteringUtilz.py, lines 446–516 · score 0.83 · Calinski Harabasz, Davies Bouldin, Silhouette score, DBI, CH, validity
  2. [2] § Results › Neurophysiological Profiles ↔ src/processEEGWithFOOOF2.m, the whole file · a weak match · score 0.79 · beta power, alpha power, gamma power, relative power, frequency bands, absolute
  3. [3] § Methods › Statistical Testing ↔ ClusteringUtilz.py, lines 310–442 · score 0.79 · chi squared, way ANOVA, FDR, Tukey, variables, numerical
  4. [4] § Methods › Feature Engineering ↔ FeatureExtraction.py, lines 593–670 · score 0.77 · Lempel Ziv complexity, permutation entropy, PermEn, LZC, mapping, temporal
  5. [5] § Methods › EEG Acquisition and Preprocessing ↔ src/remove_chans.m, the whole file · a weak match · score 0.77 · bandpass filtering, channel rejection, EEGLAB, eyes, MATLAB, 40 Hz
  6. [6] § Methods › Statistical Testing ↔ Resting State - Classic - DBSCAN.ipynb, lines 47–71 · score 0.74 · gait speed, disease duration, MoCA, CTT1, CTT2, age
  7. [7] § Methods › Feature Engineering ↔ src/processEEGWithFOOOF2.m, the whole file · a weak match · score 0.72 · relative power, frequency bands, spectrum, gamma, algorithm, domains
  8. [8] § Methods › EEG Acquisition and Preprocessing ↔ src/flag_and_remove_artifacts.m, the whole file · a weak match · score 0.67 · independent component, EEGLAB, artifact, eyes, ICA, channel
  9. [9] § Results › Feature Importance ↔ FeatureExtraction.py, lines 593–670 · score 0.59 · Lempel Ziv complexity, permutation entropy, PermEn, LZC
  10. [10] § Results › Clustering Performance ↔ ClusteringUtilz.py, lines 446–516 · score 0.56 · CH scores, validity scores, DBI, Silhouette, noise, DBSCAN
  11. [11] § Results › Clinical Profiles ↔ Resting State - Classic - DBSCAN.ipynb, lines 47–71 · score 0.56 · gait speed, disease duration, age, III, LEDD, UPDRS
  12. [12] § Results › Clinical Profiles ↔ Resting State - Classic - GMM.ipynb, lines 46–62 · score 0.56 · gait speed, disease duration, age, III, LEDD, UPDRS
  13. [13] § Methods › Clinical Profile ↔ Resting State - Classic - GMM.ipynb, lines 46–62 · score 0.56 · disease duration, MoCA, age, III, LEDD, UPDRS
  14. [14] § Methods › Dimensionality Reduction ↔ Resting State - Classic - DBSCAN.ipynb, lines 99–122 · score 0.51 · latent dimensionality, UMAP, PCA, component

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 · 1,118 lines · 51 KB · MIT · 3 matches

  1. import os
  2. import mne
  3. import scipy
  4. import random
  5. import joblib
  6. import numpy as np
  7. import pandas as pd
  8. from tqdm import tqdm
  9. import seaborn as sns
  10. import umap.umap_ as umap
  11. import matplotlib.pyplot as plt
  12. from itertools import combinations
  14. from sklearn.impute import SimpleImputer
  15. from sklearn.mixture import GaussianMixture
  16. from sklearn.model_selection import train_test_split, KFold, LeaveOneOut, BaseCrossValidator
  17. from sklearn.base import clone
  18. from sklearn.metrics import (
  19. calinski_harabasz_score,
  20. davies_bouldin_score,
  21. silhouette_score,
  22. adjusted_rand_score
  23. )
  25. from statsmodels.stats.multitest import multipletests
  26. from statsmodels.stats.multicomp import pairwise_tukeyhsd
  27. from scipy.stats import f_oneway, chi2_contingency
  30. SEED = 42
  33. class LeaveOneOutTrainEval(BaseCrossValidator):
  34. """
  35. A Leave-One-Out (LOO) like CrossValidator object that fits and scores on each N-1 samples in the data.
  36. """
  37. def get_n_splits(self, X=None, y=None, groups=None):
  38. return len(X)
  40. def split(self, X, y=None, groups=None):
  41. n = len(X)
  42. idx = np.arange(n)
  43. for holdout in range(n):
  44. train_idx = np.delete(idx, holdout)
  45. yield train_idx, train_idx
  48. class IOUtilz:
  49. @staticmethod
  50. def read_clinical_features(path: str) -> pd.DataFrame:
  51. """
  52. Reads the data from the clinical features .xlsx into a dataframe.
  53. :param path: a path to the .xlsx file of the clinical data.
  54. """
  55. data_file = pd.ExcelFile(path)
  56. cl_data = pd.concat([data_file.parse(sheet_name=sheet) for sheet in data_file.sheet_names])
  57. cl_data['IsHC'] = cl_data['Subject'].str.startswith('HC')
  58. cl_data.loc[~cl_data['IsHC'], 'Group'] = cl_data.loc[~cl_data['IsHC'], 'Group'].fillna('Idiopathic')
  59. return cl_data.reset_index(drop=True)
  62. class PreprocessUtilz:
  63. @staticmethod
  64. def normalize_col(col: pd.Series) -> float:
  65. """
  66. Performs z-score normalization of a given column (pandas series).
  67. """
  68. return (col - col.mean()) / col.std()
  70. def normalize_group_means(means: np.ndarray | pd.Series, ses: np.ndarray | pd.Series = None,
  71. mode: str = None) -> tuple[np.ndarray | pd.Series]:
  72. """
  73. Normalize the means and standard errors (ses) of a given feature across groups (e.g. clusters).
  74. :param means: a sequence representing the mean of a given feature for each group.
  75. :param ses: a sequence representing the standard error of a given feature for each group.
  76. :param mode: the name of the normalization mode, should be one of 'minmax', 'z-score'.
  77. """
  78. if mode is None:
  79. return means, ses
  81. if mode.lower() == "minmax":
  82. fmin = means.min(axis=0)
  83. rng = (means.max(axis=0) - fmin)
  84. rng = rng.replace(0, 1.0)
  85. means_n = (means - fmin) / rng
  86. ses_n = None if ses is None else ses / rng
  87. return means_n, ses_n
  89. if mode.lower() == "z-score":
  90. mu = means.mean(axis=0)
  91. sigma = means.std(axis=0).replace(0, 1.0)
  92. means_n = (means - mu) / sigma
  93. ses_n = None if ses is None else ses / sigma
  94. return means_n, ses_n
  96. raise ValueError("mode must be None, 'minmax', or 'z-score'")
  98. @staticmethod
  99. def _preprocess_pipeline(pipe, data: np.ndarray):
  100. """
  101. Transforms the input data along all pipeline steps excluding the last one (e.g. clustering model).
  102. :param pipe: a Pipeline object representing the model.
  103. :param data: an input numpy array for the model.
  104. """
  105. return pipe[:-1].transform(data)
  107. @staticmethod
  108. def set_random_seeds(models: list, seed: int = None):
  109. """
  110. Sets a random seed to all steps in a the given pipelines.
  111. :param models: a list of pipeline objects.
  112. :param seed: a random seed to set. Default is None.
  113. """
  114. for m in models:
  115. params_to_set = {}
  117. for name, step in m.steps:
  118. if hasattr(step, 'random_state'):
  119. param_key = f"{name}__random_state"
  120. params_to_set[param_key] = seed
  122. if params_to_set:
  123. m.set_params(**params_to_set)
  125. @staticmethod
  126. def get_ylim(data: pd.DataFrame, feature_name: str, groupby: list[str], offset: bool = True
  127. , error: str = 'se') -> (float, float):
  128. """
  129. Calculates scalar limits for a given feature_name in the data using groupby columns.
  130. :param data: a dataframe with the relevant columns.
  131. :param feature_name: the feature name to get scalar limits for.
  132. :param groupby: a list of column names according to which to group the data and calculate the limits.
  133. :param offset: if True, widen the scalar limits in both directions by a standard deviation offset.
  134. :param error:
  135. :return ylim: a tuple (low, high) of the scalar limits for the given feature.
  136. """
  137. groupby = data.groupby(groupby)[feature_name]
  138. means = groupby.mean()
  139. ylim = means.min(), means.max()
  141. if offset:
  142. if error == 'se':
  143. offset = groupby.sem().max()
  144. else:
  145. offset = groupby.std().max()
  146. ylim = ylim[0] - 2*offset, ylim[1] + 2*offset
  148. return ylim
  151. class TopoPlot:
  152. def __init__(self, data: pd.DataFrame, group_col: str, channel_col: str, feature_names: list[str],
  153. display_names: list[str] = None, montage: str = 'standard_1020', cmap: str = 'viridis',
  154. dims: (float, float) = (3, 4), vlim: (float, float) = None, title: str = None):
  155. """
  156. :param data: a dataframe with the raw feature table.
  157. :param group_col: the name of a column in the feature table to split the figures by.
  158. :param channel_col: the name of a column in the feature table which corresponds to the channel name \ ID.
  159. :param feature_names: a list of feature name to analyze topologically.
  160. :param display_names: the display names for the given features (labels of color bars).
  161. :param montage: a name of an MNE-EEG montage for the topological figure.
  162. :param cmap: a seaborn colormap name.
  163. :param dims: a tuple of dimensions for each figure.
  164. :param title: a customized title for the figure.
  165. """
  166. self.features = data.copy()
  167. self.feature_names = feature_names
  169. if display_names is None:
  170. self.display_names = [f.replace('_', ' ').strip() for f in feature_names]
  171. else:
  172. self.display_names = display_names
  174. self.channel_col = channel_col
  175. self.group_col = group_col
  176. self.groups = self.features[self.group_col].unique()
  177. self.n_groups = len(self.groups)
  179. self.eeg_info = TopoPlot.get_eeg_info(
  180. electrodes=self.features[self.channel_col].unique()
  181. , montage=montage
  182. )
  184. self.cmap = cmap
  185. self.dims = dims
  187. self.vlim = vlim
  188. self.title = title
  190. self.plot_topomap_mat()
  192. @staticmethod
  193. def get_eeg_info(electrodes: np.ndarray | list, montage: str):
  194. """
  195. Creates an MNE info object for downstream MNE topological plots creation.
  196. :param electrodes: a sequence of electrode names for the given configuration.
  197. :param montage: the desired configuration of electrodes on the scalp.
  198. """
  199. ch_names = list(electrodes)
  200. info = mne.create_info(ch_names=ch_names, sfreq=1, ch_types='eeg')
  201. montage = mne.channels.make_standard_montage(montage)
  202. info.set_montage(montage, on_missing='ignore')
  203. return info
  205. def plot_topomap(self, features: pd.DataFrame, feature_name: str, ax, title: str = '', vlim: (float, float) = None):
  206. """
  207. Creates figure with the topological distribution of the given feature name over the scalp using '.montage'.
  208. :param features: a transformed feature table (dataframe) with features are columns.
  209. :param feature_name: the feature name to plot its topological distribution.
  210. :param channel_col:
  211. :param ax: a matplotlib axes to draw the plot within.
  212. :param title: a title for the given axes.
  213. :param vlim: a tuple of floats (min, max) representing the color bar limits.
  214. If not given, estimated from values using PreprocessUtilz.get_ylim.
  215. """
  216. if vlim is None:
  217. vlim = PreprocessUtilz.get_ylim(features
  218. , feature_name
  219. , self.channel_col
  220. , offset=False
  221. , error='se')
  223. values = features.groupby(self.channel_col)[feature_name].mean().values
  224. im, cn = mne.viz.plot_topomap(values
  225. , self.eeg_info
  226. , axes=ax
  227. , cmap=self.cmap
  228. , vlim=vlim
  229. , show=False
  230. )
  231. ax.set_title(title, fontsize=13, y=1.02)
  232. return im
  234. def plot_topomap_row(self, feature_name: str, display_name: str, group_names: bool = True, axes = None):
  235. """
  236. Creates figure with the topological distribution of the given feature name over the scalp using '.montage'.
  237. The figure consists of a 1D row of subplots, one per each group, defined using '.group_col'.
  238. :param feature_name: the feature name to plot its topological distribution.
  239. :param display_name: the display name for the given feature (label of the color bar).
  240. :param group_names: If true, adds groups names as titles to the subplots. Default is True.
  241. :param axes: a 1D matplotlib axes for plot creation. Default is None.
  242. """
  243. if self.vlim:
  244. vlim = self.vlim
  245. else:
  246. vlim = PreprocessUtilz.get_ylim(self.features
  247. , feature_name
  248. , [self.group_col, self.channel_col]
  249. , offset=False
  250. , error='se')
  252. if axes is None:
  253. fig, axes = plt.subplots(1, self.n_groups
  254. , figsize=(self.n_groups * self.dims[0], self.dims[1] * 0.64)
  255. , constrained_layout=True)
  256. else:
  257. fig = np.ravel(axes)[0].figure
  259. for i, group in enumerate(self.groups):
  260. group_mask = self.features[self.group_col] == group
  261. if group_names:
  262. im = self.plot_topomap(self.features.loc[group_mask]
  263. , feature_name
  264. , ax=axes[i]
  265. , title=group
  266. , vlim=vlim
  267. )
  268. else:
  269. im = self.plot_topomap(self.features.loc[group_mask]
  270. , feature_name
  271. , ax=axes[i]
  272. , vlim=vlim
  273. )
  275. # Add a single colorbar to the right
  276. cbar = fig.colorbar(im, ax=axes, orientation='vertical', fraction=0.03, pad=0.04)
  277. cbar.set_label(display_name, fontsize=11)
  279. return fig
  281. def plot_topomap_mat(self):
  282. """
  283. Creates figure with the topological distribution of the given feature name over the scalp using '.montage'.
  284. The figure consists of a 2D matrix of subplots: one column per each group, defined using '.group_col',
  285. and one row for each feature, defined using '.feature_names'.
  286. """
  287. n_features = len(self.feature_names)
  289. if n_features > 1:
  290. fig, axes = plt.subplots(n_features, self.n_groups,
  291. figsize=(self.n_groups * self.dims[0], n_features * self.dims[1] * 0.64))
  293. for i, fname in enumerate(self.feature_names):
  294. self.plot_topomap_row(feature_name=fname,
  295. display_name=self.display_names[i],
  296. group_names=(i == 0),
  297. axes=axes[i])
  299. else: ## only 1 feature
  300. fig = self.plot_topomap_row(feature_name=self.feature_names[0], display_name=self.display_names[0])
  302. if self.title:
  303. fig.suptitle(self.title, fontsize=14, y=0.96)
  304. else:
  305. fig.suptitle(f'Topological Distributions Between Groups', fontsize=14, y=0.96)
  307. plt.show()
  310. class StatsUtilz:
  311. @staticmethod
  312. def anova(data: pd.DataFrame, feature_name: str, group_col: str = 'clusterID', print_line: bool = True) -> dict:
  313. """
  314. One-way ANOVA across clusters for a single numeric feature.
  315. Reports F, p, df1, df2, eta-squared (η²).
  316. :param data: a dataframe with the actual numeric data and splitter columns.
  317. :param feature_name: a name of a column within data on which to perform the analysis.
  318. :param group_col: a name of a column within data for group splitting (e.g. 'clusterID').
  319. :param print_line: if True, prints all statistics in a comprehensive line.
  320. """
  321. groups = []
  322. levels = data[group_col].dropna().unique()
  323. for lev in levels:
  324. groups.append(data.loc[data[group_col] == lev, feature_name].dropna().values)
  326. groups = [g for g in groups if len(g) > 0]
  328. F, p = f_oneway(*groups)
  330. k = len(groups)
  331. n = sum(len(g) for g in groups)
  332. df1 = k - 1
  333. df2 = n - k
  335. eta_sq = (F * df1) / (F * df1 + df2) if (F * df1 + df2) > 0 else np.nan
  337. if print_line:
  338. print(f"{feature_name}: F({df1},{df2}) = {F:.3f}, p = {p:.6f}, η² = {eta_sq:.3f}")
  340. return {
  341. "feature": feature_name,
  342. "k_groups": k,
  343. "n_total": n,
  344. "F": F,
  345. "p": p,
  346. "df1": df1,
  347. "df2": df2,
  348. "eta_sq": eta_sq
  349. }
  351. @staticmethod
  352. def pairwise_posthoc(data: pd.DataFrame, feature_name: str, group_col: str = 'clusterID', print_line: bool = True) -> str:
  353. """
  354. Conducts a tukey posthoc analysis test, assuming a significant one-way ANOVA.
  355. :param data: a dataframe with the actual numeric data and splitter columns.
  356. :param feature_name: a name of a column within data on which to perform the analysis.
  357. :param group_col: a name of a column within data for group splitting (e.g. 'clusterID').
  358. :param print_line: if True, prints all statistics in a comprehensive line.
  359. """
  360. data_c = data[[feature_name, group_col]].dropna().copy()
  361. posthoc = pairwise_tukeyhsd(endog=data_c[feature_name], groups=data_c[group_col], alpha=0.05)
  363. summary = posthoc.summary()
  364. if print_line:
  365. print(summary)
  367. return summary
  369. def chi_squared(data: pd.DataFrame, feature_name: str, group_col: str = 'clusterID', id_col: str = 'Subject',
  370. print_line: bool = True) -> dict:
  371. """
  372. Chi-square test of independence between group_col (e.g., 'clusterID') and a categorical feature.
  373. Builds a contingency table using unique counts of id_col per cell (robust to duplicates).
  374. Reports χ², df, p and Cramer's V.
  375. :param data: a dataframe with the actual numeric data and splitter columns.
  376. :param feature_name: a name of a column within data on which to perform the analysis.
  377. :param group_col: a name of a column within data for group splitting (e.g. 'clusterID').
  378. :param id_col: a column name within data of unique IDs or keys for each sample (e.g. 'Subject').
  379. :param print_line: if True, prints all statistics in a comprehensive line.
  380. """
  381. ct = (data.groupby([group_col, feature_name])[id_col]
  382. .nunique()
  383. .unstack(fill_value=0))
  385. ct = ct.loc[(ct.sum(axis=1) > 5), (ct.sum(axis=0) > 5)]
  387. chi2, p, dof, exp = chi2_contingency(ct.values, correction=False)
  389. n = ct.values.sum()
  390. r, c = ct.shape
  391. k = min(r, c)
  392. cramers_v = np.sqrt(chi2 / (n * (k - 1))) if k > 1 else np.nan
  394. if print_line:
  395. print(f"{feature_name}: χ²({dof}) = {chi2:.3f}, p = {p:.6f}, V = {cramers_v:.3f}")
  397. return {
  398. "feature": feature_name,
  399. "table": ct,
  400. "chi2": chi2,
  401. "df": dof,
  402. "p": p,
  403. "cramers_v": cramers_v,
  404. "n_total": int(n),
  405. "shape": (r, c)
  406. }
  408. def correct_ps(Ps: np.ndarray | list) -> np.ndarray:
  409. """
  410. Performs FDR-BH statistical correction for multiple comparisons.
  411. :param Ps: a sequence of p-values to correct.
  412. """
  413. _, pvals_corrected, _, _ = multipletests(Ps, alpha=0.05, method='fdr_bh')
  414. return pvals_corrected
  416. def calc_stats(data: pd.DataFrame, feature_families: dict):
  417. """
  418. Calculating statistical test score including p-values and effect-size measures for the given features.
  419. In addition, computes corrected p-values for multiple comparisons within each feature family using the FDR-BH method.
  420. """
  421. all_frames = []
  422. effect_sizes = []
  423. for family, features in feature_families.items():
  424. Ps = []
  425. for f in features:
  426. if pd.api.types.is_numeric_dtype(data[f]):
  427. res = StatsUtilz.anova(data, f)
  428. effect_sizes.append(res['eta_sq'])
  429. else:
  430. res = StatsUtilz.chi_squared(data, group_col='clusterID', feature_name=f)
  431. effect_sizes.append(res['cramers_v'])
  432. Ps.append(res['p'])
  434. _, corrected, _, _ = multipletests(Ps, alpha=0.05, method='fdr_bh')
  436. all_frames.append(pd.DataFrame(
  437. {'Variable': features, 'RawP': Ps, 'CorrectedP': corrected}
  438. ))
  440. stats = pd.concat(all_frames)
  441. stats['EffectSize'] = effect_sizes
  442. return stats
  445. class ClusterEval:
  446. @staticmethod
  447. def get_intrinsic_metrics(data: np.ndarray, models: list, model_names: list[str]) -> pd.DataFrame:
  448. """
  449. Get intrinsic clustering metrics to evaluate the model.
  450. The metrics calculated are: silhouette score, Calniski-Harbasz index (CH) and Davies-Bouldin index (DBI).
  451. In addition, for the sake of DBSCAN evaluation, the function computes validity scores - defined by
  452. wether there are more than 1 non-noise labels. For other models, you may ignore that.
  453. :return: a dataframe with the intrinsic metrics calculated for each model.
  454. """
  455. n_models = len(models)
  456. loo = LeaveOneOut()
  458. test_shape = (n_models, data.shape[0])
  459. silhouette_scores = np.zeros(test_shape)
  460. CH_scores = np.zeros(test_shape)
  461. DBI_scores = np.zeros(test_shape)
  462. validity_scores = np.zeros(test_shape)
  464. for i, model in enumerate(models):
  465. for train_idx, j in loo.split(data):
  466. sample = data[train_idx]
  467. labels = model.fit_predict(sample)
  469. non_noise_mask = labels != -1
  470. non_noise_labels = labels[non_noise_mask]
  471. non_noise_data = sample[non_noise_mask]
  473. if np.unique(non_noise_labels).size > 1: # more than 1 cluster
  474. prep_sample = PreprocessUtilz._preprocess_pipeline(model, non_noise_data)
  475. silhouette_scores[i, j] = silhouette_score(prep_sample, non_noise_labels)
  476. CH_scores[i, j] = calinski_harabasz_score(prep_sample, non_noise_labels)
  477. DBI_scores[i, j] = davies_bouldin_score(prep_sample, non_noise_labels)
  478. validity_scores[i, j] = 1
  480. else: # 1 cluster or None
  481. silhouette_scores[i, j] = np.nan
  482. CH_scores[i, j] = np.nan
  483. DBI_scores[i, j] = np.nan
  485. silhouette_means = np.nanmean(silhouette_scores, axis=1)
  486. CH_means = np.nanmean(CH_scores, axis=1)
  487. DBI_means = np.nanmean(DBI_scores, axis=1)
  488. validity_means = np.nanmean(validity_scores, axis=1)
  490. silhouette_stds = np.nanstd(silhouette_scores, axis=1)
  491. CH_stds = np.nanstd(CH_scores, axis=1)
  492. DBI_stds = np.nanstd(DBI_scores, axis=1)
  493. validity_stds = np.nanstd(validity_scores, axis=1)
  495. # Creating a DataFrame with the metrics
  496. metric_names = ['silhouette_score', 'CH_score', 'DBI_score', 'Validity']
  497. if model_names is None:
  498. model_names = list(range(n_models))
  500. intrinsic_metrics = pd.DataFrame({
  501. 'ModelName': len(metric_names) * model_names
  502. , 'Metric': np.concatenate([n_models * [mname] for mname in metric_names])
  503. , 'Mean': np.concatenate([silhouette_means, CH_means, DBI_means, validity_means])
  504. , 'STD': np.concatenate([silhouette_stds, CH_stds, DBI_stds, validity_stds])
  505. })
  507. pivoted = intrinsic_metrics.melt(
  508. id_vars=['ModelName', 'Metric'],
  509. var_name='Statistic',
  510. value_name='Value'
  511. ).pivot_table(
  512. index='Metric',
  513. columns=['ModelName', 'Statistic'],
  514. values='Value'
  515. )
  516. return pivoted
  518. @staticmethod
  519. def test_stability(data: np.ndarray, model_path: str, n_trials: int = 100, random_seed: bool = True) -> np.ndarray:
  520. """
  521. Test within-model stability between clustering predictions using the adjusted rand index (ARI)
  522. in a leave-one-out paradigm, with concurrent validation trials.
  523. :param data: an input numpy array for the model.
  524. :param model_path: a path to .pkl files with the clustering model.
  525. This function assumes that the model is a Pipeline object.
  526. :param n_trials: the number of vaidation trials to perform.
  527. :param random_seed: if True, randomizes a new random seed for each iteration.
  528. If False, sets all random seeds to be None.
  529. """
  530. model_name = model_path.split(os.sep)[-1].split('.')[0]
  531. model_1 = joblib.load(model_path)
  532. model_2 = joblib.load(model_path)
  534. # Resetting existing random seeds in each step for each model
  535. if not random_seed:
  536. PreprocessUtilz.set_random_seeds([model_1, model_2], seed=None)
  538. n = data.shape[0]
  539. loo = LeaveOneOut()
  540. ARIs = np.zeros((n_trials, n))
  542. # Iterating on subjects in a Leave-One-Out paradigm
  543. for train_idx, test_idx in tqdm(loo.split(data)):
  544. sample = data[train_idx]
  546. # Iterating n_trials times
  547. for t in range(n_trials):
  548. # Specifing distinct random seeds for each model
  549. if random_seed:
  550. seed_1 = random.randint(1, 1e6)
  551. PreprocessUtilz.set_random_seeds([model_1], seed=seed_1)
  553. seed_2 = random.randint(1, 1e6)
  554. PreprocessUtilz.set_random_seeds([model_2], seed=seed_2)
  556. ARIs[t][test_idx] = adjusted_rand_score(model_1.fit_predict(sample), model_2.fit_predict(sample))
  558. return ARIs
  560. @staticmethod
  561. def test_agreement(data: np.ndarray, model_paths: list[str], n_trials: int = 100, random_seed: bool = True) -> dict:
  562. """
  563. Test pair-wise model agreement between clustering predictions using the adjusted rand index (ARI)
  564. in a leave-one-out paradigm, with concurrent validation trials.
  565. :param data: an input numpy array for the model.
  566. :param model_paths: a list of paths to .pkl files with the clustering models.
  567. This function assumes that each model is a Pipeline object.
  568. :param n_trials: the number of vaidation trials to perform.
  569. :param random_seed: if True, randomizes a new random seed for each iteration.
  570. If False, sets all random seeds to be None.
  571. """
  572. model_names = [path.split(os.sep)[-1].split('.')[0] for path in model_paths]
  573. models = [joblib.load(path) for path in model_paths]
  575. # Resetting existing random seeds in each step for each model
  576. if not random_seed:
  577. PreprocessUtilz.set_random_seeds(models, seed=None)
  579. n = data.shape[0]
  580. loo = LeaveOneOut()
  581. model_pairs = list(combinations(model_names, 2))
  582. pair_names = ["_".join(pair) for pair in model_pairs]
  583. ARIs = {pair : np.zeros((n_trials, n)) for pair in pair_names}
  585. # Iterating on subjects in a Leave-One-Out paradigm
  586. for train_idx, test_idx in tqdm(loo.split(data)):
  587. sample = data[train_idx]
  589. # Iterating n_trials times
  590. for t in range(n_trials):
  591. if random_seed:
  592. seed = random.randint(1, 1e6)
  593. PreprocessUtilz.set_random_seeds(models, seed=seed)
  595. labels = {name: None for name in model_names}
  597. # Predicting cluster labels for each model
  598. for i, m in enumerate(models):
  599. name = model_names[i]
  600. labels[name] = m.fit_predict(sample)
  602. # Computing pair-wise ARI
  603. for i, (a, b) in enumerate(model_pairs):
  604. ARIs[pair_names[i]][t][test_idx] += adjusted_rand_score(labels[a], labels[b])
  606. return ARIs
  608. @staticmethod
  609. def radar_plot(data: pd.DataFrame, metrics: list[str], group_col: str = 'clusterID', title: str = None,
  610. normalize: str = None, ylabel: str = None, figsize: (int, int) = (5, 6), yticks: np.ndarray = None,
  611. metric_names: list[str] = None, show_err: bool = True, spin = np.pi / 2):
  612. """
  613. A general function to create a radar plot including several feature.
  614. :param data: a dataframe with the feature values.
  615. :param metrics: a list of column names within data correspoding to the relevant features for the plot.
  616. :param group_col: the name of the group column within clusters. Default is 'clusterID'.
  617. :param title: a title for the figure. Default is None.
  618. :param normalize: if True, normalizing the values within each features using Z-score before creating the plot.
  619. :param ylabel: a label for the y-axis. Default is None.
  620. :param figsize: the figure size. Default is (5, 6).
  621. :param yticks: an array of specific ticks for the y-axis. Default is None.
  622. :param metric_names: a list of names correspoding to metrics to display in the plot.
  623. :param show_err: if True, draws the standard error bar for each group. Default is True.
  624. :parm spin: a rotation factor for the plot in radians. Default is np.pi / 2.
  625. """
  626. groupby = data.groupby([group_col])[metrics]
  627. means = groupby.mean().sort_index()
  628. ses = groupby.sem().sort_index()
  630. means, ses = PreprocessUtilz.normalize_group_means(means=means, ses=ses, mode=normalize)
  632. angles = np.linspace(0, 2*np.pi, len(metrics), endpoint=False)
  633. angles = np.concatenate([angles, angles[:1]])
  635. fig = plt.figure(figsize=figsize)
  636. ax = plt.subplot(111, polar=True) # << this line makes it a radar chart
  637. ax.set_xticks(angles[:-1])
  638. ax.set_xticklabels(metrics if metric_names is None else metric_names, fontsize=14)
  639. ax.tick_params(axis="x", pad=15) # push feature labels outward a bit
  641. if yticks is not None:
  642. ax.set_yticks(yticks)
  644. ax.yaxis.grid(True, alpha=0.6)
  646. if ylabel is None:
  647. ylabel = 'Value' if normalize is None else normalize
  648. ax.set_ylabel(ylabel, labelpad=50, fontsize=11)
  649. ax.yaxis.set_label_coords(0.5, -0.2)
  650. ax.yaxis.label.set_rotation(0)
  652. for cid, row in means.iterrows():
  653. vals = row.values
  654. curve = np.concatenate([vals, vals[:1]])
  655. ax.plot(angles, curve, linewidth=2, label=cid)
  657. if show_err:
  658. se = ses.loc[cid].values
  659. lower = np.concatenate([vals - se, [vals[0] - se[0]]])
  660. upper = np.concatenate([vals + se, [vals[0] + se[0]]])
  661. ax.fill_between(angles, lower, upper, alpha=0.15)
  663. ax.legend(loc="upper left", bbox_to_anchor=(1.05, 1.05), title="Cluster")
  664. ax.set_theta_offset(spin)
  665. ax.set_theta_direction(-1)
  666. ax.set_title(title, pad=16, fontsize=14, loc='center')
  667. plt.show()
  669. @staticmethod
  670. def plot_multiple_features(features: pd.DataFrame, clusters: pd.DataFrame, raw_feature_names: list[str], disp_feature_names: list[str] = None,
  671. xlabel: str = 'Feature Name', ylabel: str = 'Value', cluster_name: str = 'K-Means', title: str = 'Feature Values',
  672. plot_type: str = 'line', ylim = None, group_col: str = 'clusterID'):
  673. """
  674. A general function to plot inter-group differences of several features, between behavioral conditions.
  675. :param features: a dataframe with the feature values for each subject.
  676. :param clusters: a dataframe with the cluster assignment for each subject.
  677. :param raw_feature_names: a list of feature names corresponding to colums in features.
  678. :param disp_feature_names: a list of feature names for display in the figure, corresponding to raw_feature_names.
  679. If None, raw_feature_names are used.
  680. :param xlabel: a label for the x-axis. Default is 'Feature Name'.
  681. :param ylabel: a label for the y-axis. Default is 'Value'.
  682. :param cluster_name: the name of the clustering model, or otherwise grouping scheme for the subjects (e.g. K-Means).
  683. :param title: a title prefix for the figure. Full title will be '{title} between {cluster_name} Clusters Across Behavioral States'.
  684. :param plot_type: a name of a plot type, should be one of ['bar', 'line'].
  685. :param ylim: a tuple of limits for the y-axis in the figure. If None, limits are set automatically via seaborn.
  686. :param group_col: the name of the group column within clusters. Default is 'clusterID'.
  687. """
  688. if disp_feature_names is None:
  689. disp_feature_names = raw_feature_names
  691. selected_features = features.loc[features['FeatureName'].isin(raw_feature_names)]
  692. agged_features = selected_features.groupby(['Subject', 'Condition', 'FeatureName'])[['RawValue']].mean().reset_index()
  694. agged_features = pd.merge(agged_features, clusters, how='inner', on='Subject')
  696. agged_features['FeatureName'] = pd.Categorical(
  697. agged_features['FeatureName'],
  698. categories=raw_feature_names,
  699. ordered=True
  700. )
  702. fig, axes = plt.subplots(1, 2, figsize=(14, 5))
  704. if plot_type == 'line':
  705. sns.lineplot(agged_features.loc[agged_features['Condition'] == 'sit'],
  706. x='FeatureName', y='RawValue', hue=group_col, ax=axes[0])
  708. sns.lineplot(agged_features.loc[agged_features['Condition'] == 'walk'],
  709. x='FeatureName', y='RawValue', hue=group_col, ax=axes[1])
  710. elif plot_type == 'bar':
  711. sns.barplot(agged_features.loc[agged_features['Condition'] == 'sit'],
  712. x='FeatureName', y='RawValue', hue=group_col, ax=axes[0], capsize=0.2)
  714. sns.barplot(agged_features.loc[agged_features['Condition'] == 'walk'],
  715. x='FeatureName', y='RawValue', hue=group_col, ax=axes[1], capsize=0.2)
  717. axes[0].set_title('Resting-State')
  718. axes[1].set_title('Active-Walking')
  720. for i in range(2):
  721. axes[i].set_xlabel(xlabel)
  722. axes[i].set_ylabel(ylabel)
  723. axes[i].set_xticklabels(disp_feature_names)
  725. if ylim is not None:
  726. axes[i].set_ylim(ylim)
  729. plt.suptitle(f'{title} between {cluster_name} Clusters Across Behavioral States', fontsize=14, y=1.02)
  730. plt.show()
  732. @staticmethod
  733. def plot_single_feature(features: pd.DataFrame, clusters: pd.DataFrame, feature_name: str,
  734. ylabel: str = 'Value', cluster_name: str = 'K-Means', title: str = 'Feature Values',
  735. plot_type: str = 'line', ylim = None, group_col: str = 'clusterID'):
  736. """
  737. A general function to plot inter-group differences of a single feature, between behavioral conditions.
  738. :param features: a dataframe with the feature values for each subject.
  739. :param clusters: a dataframe with the cluster assignment for each subject.
  740. :param feature_name: a name of a feature corresponding to colums in features.
  741. :param ylabel: a label for the y-axis. Default is 'Value'.
  742. :param cluster_name: the name of the clustering model, or otherwise grouping scheme for the subjects (e.g. K-Means).
  743. :param title: a title prefix for the figure. Full title will be '{title} between {cluster_name} Clusters Across Behavioral States'.
  744. :param plot_type: a name of a plot type, should be one of ['bar', 'line'].
  745. :param ylim: a tuple of limits for the y-axis in the figure. If None, limits are set automatically via seaborn.
  746. :param group_col: the name of the group column within clusters. Default is 'clusterID'.
  747. """
  748. selected_features = features.loc[features['FeatureName'] == feature_name]
  749. agged_features = selected_features.groupby(['Subject', 'Condition'])[['RawValue']].mean().reset_index()
  750. agged_features = pd.merge(agged_features, clusters, how='inner', on='Subject')
  752. if plot_type == 'line':
  753. ax = sns.lineplot(agged_features, x='Condition', y='RawValue', hue=group_col)
  755. elif plot_type == 'bar':
  756. ax = sns.barplot(agged_features, x='Condition', y='RawValue', hue=group_col, capsize=0.2)
  758. plt.ylabel(ylabel)
  760. if ylim is not None:
  761. plt.ylim(ylim)
  763. ax.set_xticklabels(['Resting-State', 'Active-Walking'])
  765. plt.title(f'{title} between {cluster_name} Clusters Across Behavioral States', fontsize=12, y=1.02)
  766. plt.show()
  769. class ClusterReport:
  770. def __init__(self, models: list[GaussianMixture], model_names: list[str], PD_data: np.ndarray, PD_subjects: np.ndarray
  771. , HC_data: np.ndarray, cl_data_path: str, cl_fnames: list[str], k: int = 5, resample_rate: float = 0.8):
  772. """
  773. An object for fast evaluation of Parkinson's disease clustering models using conventional performance indices and clinical assessments.
  774. :param models: a list of unfitted sklearn clustering models which one wish to evaluate.
  775. :param model_names: labels corresponding to the models, for referral to results.
  776. :param PD_data: a list of feature matrices of the Parkinson's Disease (PD) subjects, for each model.
  777. :param PD_subjects: an array of the subjectIDs for the Parkinson's Disease (PD) subjects.
  778. :param HC_data: a list of feature matrices of the Parkinson's Disease (PD) subjects, for each model.
  779. :param cl_data_path: a path to the .xlsx file of the clinical data.
  780. :param cl_fnames: a list of feature names from the clinical data to include in evaluation (e.g. LEDD, UPDRS).
  781. """
  782. self.models = models
  783. self.model_names = model_names
  784. self.n_models = len(models)
  786. self.PD_data = PD_data
  787. self.PD_subjects = PD_subjects
  788. self.HC_data = HC_data
  790. self.cl_data: pd.DataFrame = IOUtilz.read_clinical_features(cl_data_path)
  791. self.cl_fnames = cl_fnames
  793. self.loo = LeaveOneOut()
  795. def set_models(self, models: list, model_names: list[str]):
  796. """
  797. Set the model attributes for object reuse purposes.
  798. """
  799. self.models = models
  800. self.model_names = model_names
  801. self.n_models = len(models)
  803. def set_data(self, PD_data: np.ndarray, HC_data: np.ndarray):
  804. """
  805. Set the data attributes for object reuse purposes.
  806. """
  807. self.PD_data = PD_data
  808. self.HC_data = HC_data
  810. def _fit_models(self):
  811. """
  812. Fit all models to all Parkinson's Disease (PD) data.
  813. """
  814. for i in range(self.n_models):
  815. self.models[i].fit(self.PD_data)
  817. def _get_intrinsic_metrics(self, data: np.ndarray) -> pd.DataFrame:
  818. """
  819. Get intrinsic clustering metrics to evaluate the model.
  820. The metrics calculated are: silhouette score, Calniski-Harbasz index (CH) and Davies-Bouldin index (DBI).
  821. In addition, for the sake of DBSCAN evaluation, the function computes validity scores - defined by
  822. wether there are more than 1 non-noise labels. For other models, you may ignore that.
  823. :return: a dataframe with the intrinsic metrics calculated for each model.
  824. """
  825. return ClusterEval.get_intrinsic_metrics(data=self.PD_data,
  826. models=self.models,
  827. model_names=self.model_names)
  829. def _get_clinical_means(self, fnames: list[str]) -> pd.DataFrame:
  830. """
  831. Get the means of the clinical measures for each cluster in each model.
  832. Assumes models are already fitted to the data.
  833. :param fnames: a list of feature names from the clinical data to include in calculation.
  834. :return: a dataframe with the mean of each clinical measure in fnames, per cluster and model.
  835. """
  836. # Slicing only PD subjects found in subjects
  837. PD_subject_mask = self.cl_data['Subject'].isin(self.PD_subjects)
  838. cl_features = self.cl_data.loc[PD_subject_mask, fnames]
  839. cl_feature_mat = cl_features.values
  840. cl_feature_mat[np.isnan(cl_feature_mat)] = 0
  842. # Slicing only PD subjects found in clinical data
  843. subject_mask = np.where(np.isin(self.PD_subjects, self.cl_data['Subject'].values))[0]
  845. metrics_per_model = []
  846. N = subject_mask.size
  847. for i in range(self.n_models):
  848. data = self.PD_data[subject_mask]
  849. cl_metrics = pd.DataFrame(cl_feature_mat, columns=fnames)
  850. cl_metrics['ClusterID'] = self.models[i].fit_predict(data)
  851. cl_metrics['Model'] = self.model_names[i]
  852. metrics_per_model.append(cl_metrics)
  854. # Adding healthy controls for reference
  855. HC_subject_mask = self.cl_data['Subject'].str.startswith('HC')
  856. cl_features = self.cl_data.loc[HC_subject_mask, fnames]
  857. cl_feature_mat = cl_features.values
  858. cl_feature_mat[np.isnan(cl_feature_mat)] = 0
  859. for i in range(self.n_models):
  860. cl_metrics = pd.DataFrame(cl_feature_mat, columns=fnames)
  861. cl_metrics['ClusterID'] = 'HC'
  862. cl_metrics['Model'] = self.model_names[i]
  863. metrics_per_model.append(cl_metrics)
  865. all_clinical_metrics = pd.concat(metrics_per_model)
  866. return all_clinical_metrics.groupby(['Model', 'ClusterID']).mean().reset_index(drop=False)
  868. def _visualize_metrics(self, metrics: pd.DataFrame):
  869. """
  870. Visualize metrics in table in heatmap-figure forms, grouped by model.
  871. :param metrics: a dataframe with a Model CLusterID columns, and other numerical columns to be visualized.
  872. """
  873. select_cols = [col for col in metrics.columns if col != 'Model']
  875. if self.n_models == 1: # Only one model to evaluate
  876. model_cl_data = metrics[select_cols].set_index('ClusterID')
  877. norm_cl_data = model_cl_data.apply(PreprocessUtilz.normalize_col, axis=0)
  878. plt.figure(figsize=(8, 3))
  879. sns.heatmap(norm_cl_data.T, cmap='coolwarm', cbar_kws ={'label': 'Normalized Value'}, annot=model_cl_data.round(3).T)
  880. plt.title(self.model_names[0].title())
  882. else: # Several models to evaluate
  883. fig, axes = plt.subplots(1, self.n_models, figsize=(self.n_models * 8, 3))
  884. for i in range(self.n_models):
  885. name = self.model_names[i]
  886. model_cl_data = metrics.loc[metrics['Model'] == name, select_cols].set_index('ClusterID')
  887. norm_cl_data = model_cl_data.apply(PreprocessUtilz.normalize_col, axis=0)
  888. sns.heatmap(norm_cl_data.T, ax=axes[i], cmap='coolwarm', cbar_kws ={'label': 'Normalized Value'}, annot=model_cl_data.round(3).T)
  889. axes[i].set_title(name.title())
  891. plt.suptitle('Metrics Between Clusters')
  893. plt.show()
  894. plt.close()
  896. @staticmethod
  897. def _UMAP_project(data) -> np.ndarray:
  898. """
  899. Project the given data into a lower dimension using the UMAP algorithm.
  900. :param data: a dataset to evaluate models on.
  901. """
  902. reducer = umap.UMAP(n_components=2, n_neighbors=5, random_state=SEED)
  903. return reducer.fit_transform(data)
  905. def _UMAP_visualize(self):
  906. """
  907. Creates a scatter plot for each model where UMAP-projected data points are colored by hard-clustering labels.
  908. Assumes models are already fitted to the data.
  909. :param data: a dataset to evaluate models on.
  910. """
  911. plt.figure()
  912. if self.n_models == 1:
  913. prepped = PreprocessUtilz._preprocess_pipeline(self.models[0], self.PD_data)
  914. UMAP_proj = ClusterReport._UMAP_project(prepped)
  915. x = UMAP_proj[:, 0]
  916. y = UMAP_proj[:, 1]
  917. labels = self.models[0].fit_predict(self.PD_data)
  918. scatter = plt.scatter(x=x, y=y, c=labels)
  919. plt.title('UMAP Projected Clusters (PD only)')
  920. plt.xlabel('UMAP-comp1')
  921. plt.ylabel('UMAP-comp2')
  922. legend = plt.legend(*scatter.legend_elements(), title="ClusterID")
  923. plt.gca().add_artist(legend)
  924. else:
  925. fig, axes = plt.subplots(1, self.n_models, figsize=(self.n_models * 6, 5))
  926. for i, model in enumerate(self.models):
  927. prepped = PreprocessUtilz._preprocess_pipeline(model, self.PD_data)
  928. UMAP_proj = ClusterReport._UMAP_project(prepped)
  929. x = UMAP_proj[:, 0]
  930. y = UMAP_proj[:, 1]
  931. labels = model.fit_predict(self.PD_data)
  932. scatter = axes[i].scatter(x=x, y=y, c=labels)
  933. axes[i].set_title(self.model_names[i].title(), fontsize=14)
  934. axes[i].set_xlabel('UMAP-comp1')
  935. axes[i].set_ylabel('UMAP-comp2')
  936. legend = axes[i].legend(*scatter.legend_elements(), title="ClusterID")
  937. axes[i].add_artist(legend)
  939. plt.suptitle('UMAP Projected Clusters (PD only)', fontsize=18)
  941. plt.show()
  942. plt.close()
  944. def report(self, include_clinical: bool = True):
  945. print("======================================== Model ClusterEval Report ========================================\n")
  947. print("--------------------------------------------- General Metrics --------------------------------------------\n")
  948. intrinsic_metrics = self._get_intrinsic_metrics(self.PD_data)
  949. print(intrinsic_metrics, '\n')
  951. print("-------------------------------------- Low Dimensional Visualization -------------------------------------\n")
  952. self._UMAP_visualize()
  954. if include_clinical:
  955. print("-------------------------------------- Clinical Metrics Differences -------------------------------------\n")
  956. clinical_metrics = self._get_clinical_means(self.cl_fnames)
  957. self._visualize_metrics(clinical_metrics)
  960. class FeatureImportance:
  961. def __init__(self, features_metadata: pd.DataFrame, n_trials: int, model_path: str, data: np.ndarray):
  962. """
  963. :param features_metadata:
  964. :param n_trials:
  965. :param model_path:
  966. :param data:
  967. """
  968. self.features_metadata = features_metadata
  969. self.n_trials = n_trials
  970. self.data = data
  972. self.base_model = joblib.load(model_path)
  973. self.dr_name = self.base_model.steps[-2][0]
  974. self.cluster_name = self.base_model.steps[-1][0]
  976. def create_splits(self, split_conditions: bool = True, split_regions: bool = False):
  977. """
  978. Creates data splits for feature importance analysis (via calc_feature_importance).
  979. :param split_conditions: If True, splits data between behavioral conditions. Deafult is True.
  980. :param split_conditions: If True, splits data between anatomical regions. Default is False.
  981. :return splits: a tuple (ex_ids, keep_ids) with the indices of the corresponding columns in the data to exlucde and keep.
  982. :return labels: a corresponding label for each split (feature_name, condition, region).
  983. If data was not split by conditions or regions, then condition | region = 'ALL' accordingly.
  984. """
  985. splits, labels = [], []
  987. features = self.features_metadata['FeatureName'].unique()
  988. conditions = self.features_metadata['Condition'].unique() if split_conditions else np.array(['ALL'])
  989. regions = self.features_metadata['Region'].unique() if split_regions else np.array(['ALL'])
  991. for fname in features:
  992. fmask = (self.features_metadata['FeatureName'] == fname)
  993. for cond in conditions:
  994. cmask = (self.features_metadata['Condition'] == cond)
  995. for reg in regions:
  996. rmask = (self.features_metadata['Region'] == reg)
  998. if split_conditions and split_regions: mask = fmask & cmask & rmask
  999. elif split_conditions: mask = fmask & cmask
  1000. elif split_regions: mask = fmask & rmask
  1001. else: mask = fmask
  1003. ex_idx = np.where(mask)[0]
  1004. keep_idx= np.where(~mask)[0]
  1005. splits.append((ex_idx, keep_idx))
  1006. labels.append((fname, cond, reg)) # ALWAYS 3-tuple
  1008. return splits, labels
  1010. @staticmethod
  1011. def fit_predict(model, X: np.ndarray, fit: bool = True):
  1012. """
  1013. Fits and predicts a clustering model on the input data (X), and returns the cluster labels and silhouette score.
  1014. :param model: a clustering model Pipeline object with DR (-2) and clustering (-1) steps.
  1015. :param X: the data matrix to predict.
  1016. :param fit: If True, fits the model on X, then predicts the labels on X. If False, only performs label prediction.
  1017. """
  1018. if fit: labels = model.fit_predict(X)
  1019. else: labels = model.predict(X)
  1021. embed = model[:-1].transform(X)
  1022. non_noise = labels != -1
  1023. n_clusters = len(np.unique(labels[non_noise]))
  1024. if n_clusters > 1 and embed.shape[0] > n_clusters:
  1025. sil = silhouette_score(embed[non_noise], labels[non_noise])
  1026. else:
  1027. sil = np.nan
  1029. return labels, sil
  1031. def calc_feature_importance(self, split_conditions: bool = True, split_regions: bool = False):
  1032. """
  1033. Calculates feature importance on the input data using a permutation importance approach.
  1034. :param split_conditions: If True, splits data between behavioral conditions. Deafult is True.
  1035. :param split_conditions: If True, splits data between anatomical regions. Default is False.
  1036. :return: a dataframe with the silhouette drop and ARI for each permuted data split, within each trial.
  1037. """
  1038. splits, labels = self.create_splits(split_conditions, split_regions)
  1040. shape = (len(splits), self.n_trials)
  1041. ARI = np.zeros(shape, dtype=float)
  1042. silhouette = np.zeros(shape, dtype=float)
  1044. rng = np.random.default_rng(42)
  1046. for j in tqdm(range(self.n_trials)):
  1047. seed = int(rng.integers(0, 1_000_000))
  1049. # ---- FULL baseline once per trial ----
  1050. if self.cluster_name == 'DBSCAN': seed_params = {f"{dr_name}__random_state": seed}
  1051. else: seed_params = {
  1052. f"{self.dr_name}__random_state": seed,
  1053. f"{self.cluster_name}__random_state": seed
  1054. }
  1055. raw_model = clone(self.base_model).set_params(**seed_params)
  1057. raw_labels, raw_sil = FeatureImportance.fit_predict(model=raw_model, X=self.data)
  1059. # ---- Permutation ----
  1060. for i, (ex_idx, keep_idx) in enumerate(splits):
  1061. Xp = self.data.copy()
  1062. for c in np.atleast_1d(ex_idx):
  1063. Xp[:, c] = rng.permutation(Xp[:, c])
  1065. perm_labels, perm_sil = FeatureImportance.fit_predict(model=raw_model,
  1066. X=Xp,
  1067. fit=self.cluster_name == 'DBSCAN')
  1069. ARI[i, j] = adjusted_rand_score(raw_labels, perm_labels)
  1070. silhouette[i, j] = (raw_sil - perm_sil) if (not np.isnan(raw_sil) and not np.isnan(perm_sil)) else 0.0
  1072. labels_arr = np.array(labels, dtype=object) # shape (M, 2)
  1073. feat_series = np.repeat(labels_arr[:, 0], self.n_trials) # FeatureName repeated
  1074. cond_series = np.repeat(labels_arr[:, 1], self.n_trials) # Condition repeated
  1075. reg_series = np.repeat(labels_arr[:, 2], self.n_trials) # region repeated
  1076. trials = np.tile(np.arange(1, self.n_trials + 1), len(splits))
  1078. return pd.DataFrame({
  1079. "FeatureName": feat_series,
  1080. "Condition": cond_series,
  1081. "Region": reg_series,
  1082. "Trial": trials,
  1083. "ARI": ARI.reshape(-1),
  1084. "SilhouetteDrop": silhouette.reshape(-1)
  1085. })
  1087. @staticmethod
  1088. def plot_feature_importance(importance: pd.DataFrame, order = None, xlabels = None):
  1089. """
  1090. Plots feature importance data that were pre-calculated.
  1091. :param importance: a dataframe with the feature importance results.
  1092. :param order: a list representing the ordr of features for display.
  1093. :param xlables: an alternative list of xlabels, corresponding to each feature in the given order.
  1094. """
  1095. fig, axes = plt.subplots(1, 2, figsize=(15, 5))
  1097. sns.barplot(importance, hue='Condition', x='FeatureName', y='ARI',
  1098. ax=axes[0], capsize=.2, errorbar='se', order=order)
  1099. axes[0].set_title('Adjusted Rand Index (ARI) between Raw and Permuted Data')
  1100. axes[0].set_ylim(0.7, 0.95)
  1101. axes[0].set_xlabel('Feature Name')
  1102. if xlabels is not None:
  1103. axes[0].set_xticklabels(xlabels)
  1105. axes[0].legend(loc='upper right')
  1107. sns.barplot(importance, hue='Condition', x='FeatureName', y='SilhouetteDrop',
  1108. ax=axes[1], capsize=.2, errorbar='se', order=order)
  1109. axes[1].set_title('Drop in Silhouette Score After Permutation')
  1110. axes[1].set_xlabel('Feature Name')
  1111. axes[1].set_ylabel('$\Delta$ Silhouette')
  1112. if xlabels is not None:
  1113. axes[1].set_xticklabels(xlabels)
  1115. axes[1].legend(loc='upper right')
  1117. plt.suptitle('Permutation Feature Importance between Conditions', fontsize=16, y=1.02)
  1118. plt.show()

ClusteringUtilz.py at commit e8c6780, under MIT · at the source

Overview

Authors: Daniel Vered1,2, Zoya Katzir1,3, Idan Daniel Grosbard2, Avner Thaler1,2,3, Inbal Maidan1,2,3
ORCID iDs: Avner Thaler
  1. Laboratory of Early Markers of Neurodegeneration, Neurological Institute, Tel Aviv Sourasky Medical Center, Tel Aviv, Israel
  2. Sagol School of Neuroscience, Tel Aviv University, Tel Aviv, Israel
  3. Department of Neurology, Gray Faculty of Medical and Health Sciences, Tel Aviv University, Tel Aviv, Israel
Institutions: Tel Aviv University (Israel); Tel Aviv Sourasky Medical Center (Israel)
Journal: Movement disorders : official journal of the Movement Disorder Society, volume 41, issue 8, pages 2107-2116
Dates: received 27 January 2026; accepted 16 April 2026; published online 6 May 2026; in print August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1002/mds.70348 · PMID 42089402 · PMCID PMC13518290 · OpenAlex W7160381131
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), human (organism), Parkinson's (population)
Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Complexity
Keywords: Parkinson's disease, electroencephalography (EEG), neurophysiological subtypes, clustering analysis, neural heterogeneity
MeSH: Brain*, Electroencephalography*, Parkinson Disease*, Aged, Cluster Analysis, Clustering Algorithms, Female, Humans, Leucine-Rich Repeat Serine-Threonine Protein Kinase-2, Male, Middle Aged (* major topic)
Topic: Neurological disorders and treatments (Neurology, Medicine), according to OpenAlex
Funding: Israel Science Foundation (1157/20)
Citations: cited by 1 paper (Europe PMC); 48 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 14 matches between paragraphs and lines of code.

ronmon1994/eeg-preprocessing-toolkit

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: ce581702def5bb04ca8b92b9622dcd1f44989690, 28 January 2025
Languages: MATLAB (11)
Size: 12 files, 11 scripts
Software Heritage: not archived
Found in: the text, “EEG Acquisition and Preprocessing”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: EEGLAB (9 files), ICLabel (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
12 files
At the source: github.com/RonMon1994/EEG-Preprocessing-Toolkit/tree/main

danielvered/eeg-based-pd-subtyping

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: e8c6780629c56360123ce0f0dc44e7af41780667, 3 May 2026
Languages: Python (3), Jupyter (3)
Size: 18 files, 6 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, environment (requirements.txt), 3 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (6 files), pandas (5 files), scikit-learn (5 files), Matplotlib (4 files), SciPy (4 files), seaborn (4 files), UMAP (4 files), MNE-Python (3 files), specparam (formerly FOOOF) (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
8 files
At the source: github.com/DanielVered/EEG-based-PD-Subtyping/tree/main

The paper's code and data availability statement is in the Data section.

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

Code and data availability statement

The paper has a code and 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.1002/mds.70348.

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 2, 28 September 2026

  • Publisher: n/a → Wiley

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 5 keywords, 11 MeSH terms, 1 funder, 40 references.

Cite

This paper

Vered, D., Katzir, Z., Grosbard, I. D., Thaler, A., & Maidan, I. (2026). Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease. Movement disorders : official journal of the Movement Disorder Society, 41(8), 2107-2116. https://doi.org/10.1002/mds.70348

BibTeX

@article{vered2026electroencephalograp,
author = {Vered, Daniel and Katzir, Zoya and Grosbard, Idan Daniel and Thaler, Avner and Maidan, Inbal},
title = {{Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease}},
journal = {Movement disorders : official journal of the Movement Disorder Society},
year = {2026},
month = may,
volume = {41},
number = {8},
pages = {2107--2116},
publisher = {Wiley},
issn = {0885-3185},
doi = {10.1002/mds.70348},
url = {https://doi.org/10.1002/mds.70348},
pmid = {42089402},
pmcid = {PMC13518290}
}

RIS

TY - JOUR
AU - Vered, Daniel
AU - Katzir, Zoya
AU - Grosbard, Idan Daniel
AU - Thaler, Avner
AU - Maidan, Inbal
TI - Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease
T2 - Movement disorders : official journal of the Movement Disorder Society
J2 - Mov Disord
PY - 2026
DA - 2026/05/06
VL - 41
IS - 8
SP - 2107
EP - 2116
SN - 0885-3185
PB - Wiley
DO - 10.1002/mds.70348
UR - https://doi.org/10.1002/mds.70348
LA - en
ER -

CSL-JSON

{
"id": "10.1002/mds.70348",
"type": "article-journal",
"title": "Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease",
"container-title": "Movement disorders : official journal of the Movement Disorder Society",
"author": [
{
"family": "Vered",
"given": "Daniel"
},
{
"family": "Katzir",
"given": "Zoya"
},
{
"family": "Grosbard",
"given": "Idan Daniel"
},
{
"family": "Thaler",
"given": "Avner"
},
{
"family": "Maidan",
"given": "Inbal"
}
],
"container-title-short": "Mov Disord",
"volume": "41",
"issue": "8",
"page": "2107-2116",
"DOI": "10.1002/mds.70348",
"PMID": "42089402",
"PMCID": "PMC13518290",
"ISSN": "0885-3185",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/mds.70348",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
6
]
]
}
}

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.3389/fncom.2026.1786996 [code]
Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: specparam (formerly FOOOF), UMAP, MNE-Python, 7 other tools, EEG, 1 reference
[2] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 7 other tools, EEG, 1 reference
[3] doi:10.1093/cercor/bhag113 [code]
Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 6 other tools, EEG, 1 reference
[4] doi:10.1162/imag.a.1169 [code]
Diazepam alters the shape of alpha oscillations recorded from human cortex using EEG.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: specparam (formerly FOOOF), EEGLAB, MNE-Python, 5 other tools, EEG, 2 references
[5] doi:10.1371/journal.pcbi.1014043 [code]
EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data.
Journal: PLoS computational biology
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 4 other tools, EEG, 2 references
[6] doi:10.1097/j.pain.0000000000004044 [code]
No effect of rhythmic visual stimulation on experimental pain perception.
Journal: Pain
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 6 other tools, EEG
[7] doi:10.1093/cercor/bhag077 [code]
The longitudinal development of intrinsic timescales in infancy and their relation to alpha brain rhythm.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: EEGLAB, MNE-Python, statsmodels, 6 other tools, EEG, 2 references
[8] doi:10.7554/elife.107088 [code]
Development of auditory and spontaneous movement responses to music over the first postnatal year.
Journal: eLife
In common: ICLabel, EEGLAB, MNE-Python, 6 other tools, EEG, 1 reference
[9] doi:10.1016/j.ebiom.2026.106375 [code]
Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs.
Journal: EBioMedicine
In common: ICLabel, EEGLAB, MNE-Python, 6 other tools, EEG
[10] doi:10.1371/journal.pcbi.1014154 [code]
Complexity of resting cortical activity predicts neurophysiological responses to theta-burst stimulation but fails to generalize: A rigorous machine-learning approach.
Journal: PLoS computational biology
In common: specparam (formerly FOOOF), MNE-Python, statsmodels, 5 other tools, EEG, 1 reference

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.