OSCR

Demographics-robust spontaneous eye blinking slowing in patients with severe acquired brain injury.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods and materials › EOG blink waveform processing ↔ src/eogtools/eog.py, lines 51–145 · score 0.74 · filtered signal, EOG signal, MNE, rejected, log, threshold
  2. [2] § Methods and materials › EOG blink waveform processing ↔ scripts/feature_stats.py, lines 97–147 · score 0.72 · standard deviation, blink interval variability, overlap, logarithm, onset, LIBIV
  3. [3] § Methods and materials › Statistical analysis ↔ scripts/feature_stats.py, lines 291–365 · score 0.69 · Wilcoxon signed rank, Kruskal Wallis, post hoc, eMCS, Holm, pDoC
  4. [4] § Results › Blink features ↔ scripts/feature_stats.py, lines 291–365 · score 0.66 · Wilcoxon signed rank, Kruskal Wallis, Post hoc, eMCS, KW, Holm
  5. [5] § Methods and materials › Statistical analysis ↔ src/utils/features.py, lines 915–1064 · score 0.64 · multivariate logistic regression, logistic regression models, predictor, binary, fitted, HC
  6. [6] § Methods and materials › Statistical analysis ↔ scripts/demographics_data.py, lines 225–242 · score 0.63 · Mann Whitney, Kruskal Wallis, eMCS, pDoC, Demographic, age

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 · 651 lines · 29 KB · MIT · 3 matches

  1. import os
  2. import re
  3. import traceback
  4. from collections import namedtuple
  6. import matplotlib.pyplot as plt
  7. import numpy as np
  8. import pandas as pd
  9. import yaml
  10. from dotenv import load_dotenv
  12. from utils.features import (
  13. aggregate_features,
  14. annotate_panels,
  15. export_significant_html_table,
  16. libiv,
  17. logreg_covariate,
  18. plot_logreg_effects,
  19. merge_and_check_features,
  20. plot_features_summary,
  21. plot_recording_durations,
  22. run_stat_tests_ncheck,
  23. run_stat_tests_wilcoxon,
  24. waveforms_plot,
  25. )
  27. #%% Load opt and environment variables
  29. load_dotenv() # Load environment variables from a .env file
  30. results_dir = os.getenv('RESULTS_DIR', './results')
  31. data_dir = os.getenv('DATA_DIR', './data')
  32. config_path = os.getenv('CONFIG_PATH', './configs/opt.yml')
  34. # Load options from YAML configuration file
  35. with open(os.path.abspath(config_path)) as file:
  36. opt = yaml.safe_load(file)
  37. opt = namedtuple('opt', opt.keys())(*opt.values()) # dict2namedtuple
  39. print(f"\nLoaded options from {config_path}:")
  40. _ = [print(f" - {k}: {v}") for k, v in opt._asdict().items()]
  43. def nonagg_EBR(grouped, start_time, end_time):
  44. """
  45. Compute the Eye Blink Rate (EBR) for each group in the grouped DataFrame.
  46. EBR is defined as the number of blinks per minute.
  47. """
  48. all_dfs = []
  49. for (grp,cnd, sbj), gdf in grouped:
  51. #print(grp, cnd,sbj,gdf.shape)
  52. base = f"sub-{sbj}_task-{cnd.upper()}_eog"
  53. to_load = f"{os.path.dirname(__file__)}/../data/sub-{sbj}/eog/{base}.edf"
  54. from mne.io import read_raw_edf
  56. raw = read_raw_edf(to_load,verbose='ERROR')
  57. remove_movimento_s = 0
  58. for ann in raw.annotations:
  59. if ann['description'] == 'Movimento':
  60. # Find the preceding annotation that is either 'Movimento' or 'OHMETER MEASURE'
  61. prev_ann = None
  62. for candidate in reversed([a for a in raw.annotations if a['onset'] < ann['onset']]):
  63. if candidate['description'] in ['Movimento', 'OHMETER MEASURE']:
  64. prev_ann = candidate
  65. break
  67. if prev_ann:
  68. # Calculate the delta with the preceding annotation's onset + duration
  69. delta_prev_ann = ann['onset'] - (prev_ann['onset'] + prev_ann.get('duration', 0))
  70. else:
  71. # If no preceding annotation, calculate delta with the start
  72. delta_prev_ann = ann['onset']
  74. # Add the minimum of 10 seconds or the calculated delta
  75. remove_movimento_s += min(10, max(0, delta_prev_ann))
  76. elif ann['description'] == 'OHMETER MEASURE':
  77. # Add the duration of the OHMETER MEASURE annotation
  78. remove_movimento_s += ann['duration']
  79. else:
  80. continue
  82. T_min_nocorr = ( (end_time
  83. -start_time)
  84. / 60)
  85. T_min_w_corr = ( (min(end_time, max(raw.times))
  86. -max(start_time, min(raw.times))
  87. -remove_movimento_s)
  88. / 60 )
  89. ebr_nocorr = gdf.shape[0] / T_min_nocorr
  90. ebr_w_corr = gdf.shape[0] / T_min_w_corr
  91. print(f"[{os.path.basename(to_load)}] - ({min(raw.times):.0f},{max(raw.times):.0f}) out of ({start_time},{end_time}) - {remove_movimento_s}; corr: {ebr_nocorr:.1f} >> {ebr_w_corr}")
  92. all_dfs.append(pd.DataFrame({'Group':grp,'Condition':cnd,'Subject':sbj,'EBR':ebr_w_corr},index=[0]))
  94. return pd.concat(all_dfs, ignore_index=True)
  97. def nonagg_libiv(grouped, start_time, end_time):
  98. """
  99. Compute the Inter-Blink Interval Variability (IBIV) for each group in the grouped DataFrame.
  100. IBIV is defined as the standard deviation of the logarithm of inter-blink intervals.
  101. """
  102. # libiv use bb interval from grouped; we need to remove bb intervals that
  103. # overlap with artifacts such as:
  104. # Movimento.onset-10:Movimento.onset
  105. # OHMETER.onset:OHMETER.onset+OHMETER.duration
  106. def _disjoint(a,b,strict=False):
  107. a, b = sorted(a), sorted(b)
  108. if strict:
  109. return a[1] < b[0] or b[1] < a[0]
  110. else:
  111. return a[1] <= b[0] or b[1] <= a[0]
  114. all_dfs = []
  115. for (grp,cnd, sbj), gdf in grouped:
  116. libiv_no_corr = libiv(gdf.copy())
  117. #print(grp, cnd,sbj,gdf.shape)
  118. base = f"sub-{sbj}_task-{cnd.upper()}_eog"
  119. to_load = f"{os.path.dirname(__file__)}/../data/sub-{sbj}/eog/{base}.edf"
  120. from mne.io import read_raw_edf
  122. raw = read_raw_edf(to_load,verbose='ERROR')
  123. for ann in raw.annotations:
  124. if ann['description'] == 'Movimento':
  125. a = (ann['onset'] - 10, ann['onset'])
  126. def bb_int(row):
  127. if np.isnan(row['BB']):
  128. return (row['Time'], row['Time'])
  129. else:
  130. return (row['Time'] - row['BB'], row['Time'])
  131. # filter rows of gdf
  132. gdf = gdf[gdf.apply(lambda row: _disjoint(a, bb_int(row)),
  133. axis=1)]
  134. elif ann['description'] == 'OHMETER MEASURE':
  135. a = (ann['onset'], ann['onset'] + ann['duration'])
  136. # filter rows of gdf
  137. gdf = gdf[gdf.apply(lambda row: _disjoint(a, bb_int(row)),
  138. axis=1)]
  139. else:
  140. continue
  142. libiv_w_corr = libiv(gdf)
  144. print(f"[{os.path.basename(to_load)}] - libiv corr from {libiv_no_corr} to {libiv_w_corr}")
  145. all_dfs.append(pd.DataFrame({'Group':grp,'Condition':cnd,'Subject':sbj,'LIBIV':libiv_w_corr},index=[0]))
  147. return pd.concat(all_dfs, ignore_index=True)
  151. #%% Main
  152. if __name__ == '__main__':
  153. demographics_df = pd.read_csv(os.path.join(data_dir, "demographics.csv"),dtype={"Subject": str})
  154. subject2group = demographics_df.copy().set_index("Subject")["Group"].to_dict()
  156. print(f"Working on output folder: {results_dir}")
  157. os.makedirs(results_dir, exist_ok=True)
  159. print("Checking record durations")
  160. fig = plot_recording_durations(data_dir, subject2group=subject2group,
  161. palette=opt.palette)
  162. #plt.show(block=False)
  163. fig.savefig(os.path.join(results_dir,'recording_duration.svg'))
  165. print("Merging and checking features")
  166. data = merge_and_check_features(os.path.join(results_dir,
  167. 'eog\\*_features.csv'),
  168. subject2group=subject2group,
  169. filter_similarity=opt.filter_similarity)
  171. print("Visualizing waveforms")
  172. fig = waveforms_plot(os.path.join(results_dir,
  173. 'eog\\*_waveforms.csv'),
  174. subject2group=subject2group,
  175. palette=opt.palette,
  176. filter_similarity=opt.filter_similarity
  177. )
  178. #plt.show(block=False)
  179. fig.savefig(os.path.join(results_dir,'waveforms.svg'))
  181. print('Extracting aggregated features for each subject')
  182. window_definitions = {"6mins": {'Resting': (30,30+6*60), 'Oddball': (30,30+6*60)}}
  184. features_agg_funcs = {
  185. "Mean Amplitude (µV)": ('Amplitude', 'mean'),
  186. "Mean Duration (ms)": ('Duration', 'mean'),
  187. "Mean Rise Time (ms)": ('Rise', 'mean'),
  188. "Mean Fall Time (ms)": ('Fall', 'mean'),
  189. }
  191. features_nonagg_funcs = {
  192. "EBR (blink/min)": nonagg_EBR,
  193. "LIBIV": nonagg_libiv,
  194. }
  196. # first time column 0 for resting, first stimulation for oddball
  197. records = []
  198. for sub in data['Subject'].unique():
  199. for cond in data['Condition'].unique():
  200. if cond.upper() == 'RESTING':
  201. first_time = 0
  202. elif cond.upper() == 'ODDBALL':
  203. stim_path = os.path.join(results_dir,"STIM", f"sub-{sub}_task-{cond}_eog_stims.csv")
  204. try:
  205. first_time = float(pd.read_csv(stim_path)['onset'].iloc[0])
  206. except Exception as e:
  207. print(f"Failed reading stim file for {sub} at "
  208. f"{stim_path}: {e}.")
  209. first_time = np.nan
  210. records.append({'Subject': sub, 'Condition': cond, 'FirstTime': first_time})
  211. first_time_df = pd.DataFrame.from_records(records)
  212. data = data.merge(first_time_df, on=['Subject','Condition'], how='left')
  214. dagg = [aggregate_features(data.query(f"Condition == '{cond}'"),
  215. window_name=name,
  216. time_bounds=bounds[cond],
  217. first_time_col='FirstTime',
  218. groupby_columns=['Group','Condition','Subject'],
  219. bound_window_on=['Time'],
  220. aggregate_funcs=features_agg_funcs,
  221. nonaggregate_funcs=features_nonagg_funcs
  222. )
  223. for name, bounds in window_definitions.items()
  224. for cond in bounds]
  226. dagg = pd.concat(dagg, ignore_index=True)
  228. # Keep only the 6mins window and remove Window column
  229. dagg = dagg[dagg['Window'] == '6mins'].drop(columns=['Window']).reset_index(drop=True)
  231. print("Imputation of missing data")
  232. # Imputation of missing values by group, subject, condition, window
  233. how_many_nans = dagg.isna().sum().sum()
  234. # Report data to be imputed
  235. for col in dagg.columns:
  236. if dagg[col].isna().any():
  237. nan_rows = dagg[dagg[col].isna()]
  238. for ridx, _ in nan_rows.iterrows():
  239. print(f"{pd.DataFrame(dagg.loc[ridx,['Group','Subject','Condition','Window',col]]).T}")
  240. # Impute with median stratifying by group, condition and window
  241. for col2impute in [col for col in dagg.columns if col.startswith('std_')]:
  242. # Impute missing values in the column col2impute
  243. dagg[col2impute] = dagg.groupby(['Group', 'Condition', 'Window'])[col2impute].transform(
  244. lambda x: x.fillna(x.median()))
  245. print(f"Imputed {how_many_nans} NaN values, remaining NaNs: {dagg.isna().sum().sum()}")
  247. dagg['Group'] = pd.Categorical(dagg['Group'], categories=opt.palette.keys(), ordered=True)
  248. dagg = dagg.sort_values(by=['Group', 'Subject', 'Condition'],
  249. ascending=[True, True, False])
  251. # sort columns
  252. first_cols = ['Group', 'Subject', 'Condition']
  253. dagg = dagg[ first_cols +
  254. # with EBR and LIBIV first if present
  255. [col for col in dagg.columns
  256. if col not in first_cols and
  257. col.startswith(('EBR (blink/min)', 'LIBIV'))] +
  258. # then the rest
  259. [col for col in dagg.columns
  260. if col not in first_cols and
  261. not col.startswith(('EBR (blink/min)', 'LIBIV'))]
  262. ]
  264. dagg.to_csv(os.path.join(results_dir, 'aggregated_features.csv'),
  265. index=False)
  266. print(f'Aggregated features, saved in {os.path.join(results_dir, "aggregated_features.csv")}:')
  267. print(dagg)
  270. # Join aggregated results with demographics
  271. dagg_and_demo = dagg.merge(demographics_df, on=['Group', 'Subject'], how='left')
  272. if 'Window' in dagg_and_demo.columns:
  273. dagg_and_demo = dagg_and_demo.drop(columns=['Window'])
  274. print(f"Results with demographics: {dagg_and_demo.shape[0]} rows, {dagg_and_demo.shape[1]} columns")
  276. # Save joined results
  277. out_path = os.path.join(results_dir, 'aggregated_features_with_demographics.csv')
  278. dagg_and_demo.to_csv(out_path, index=False)
  279. print(f"Saved: {out_path}")
  281. # Fold conditions into columns
  282. demo_cols = [c for c in demographics_df.columns if c in dagg_and_demo.columns]
  283. value_cols = [c for c in dagg_and_demo.columns if c not in set(['Group', 'Subject', 'Condition']) | set(demo_cols)]
  284. dagg_and_demo_folded = (
  285. dagg_and_demo.copy().fillna(-1, inplace = False) # needed to pivot
  286. .pivot_table(index=demo_cols, columns='Condition', values=value_cols, aggfunc='first')
  287. .reset_index()
  288. )
  289. dagg_and_demo_folded.replace(-1, np.nan, inplace=True)
  291. # Flatten MultiIndex columns like "<Condition> <Feature>"
  292. new_cols = []
  293. for col in dagg_and_demo_folded.columns:
  294. if isinstance(col, tuple):
  295. feat, cond = col
  296. new_cols.append(f"{cond} {feat}".strip())
  297. else:
  298. new_cols.append(col)
  299. dagg_and_demo_folded.columns = new_cols
  301. # Revert placeholder demographics back to NaN
  302. dagg_and_demo_folded[demo_cols] = dagg_and_demo_folded[demo_cols].replace(-1, pd.NA)
  304. print(f"Results with demographics and folded conditions: {dagg_and_demo_folded.shape[0]} rows, {dagg_and_demo_folded.shape[1]} columns")
  306. # Save folded results
  307. out_path_folded = os.path.join(results_dir, 'aggregated_features_folded_with_demographics.csv')
  308. dagg_and_demo_folded.to_csv(out_path_folded, index=False)
  309. print(f"Saved: {out_path_folded}")
  312. group_stats_df = run_stat_tests_ncheck(dagg,
  313. condition_col='Condition',
  314. group_col='Group',
  315. alpha=0.05,
  316. p_adjust='holm',
  317. group_comparisons=[('HC', 'eMCS'),
  318. ('eMCS', 'pDoC'),
  319. ('HC', 'pDoC')])
  321. print('\nStatistical test results:\n')
  322. group_stats_df = group_stats_df.loc[:, [col for col in group_stats_df.columns
  323. if col[1] != 'Statistic']]
  324. group_stats_df = group_stats_df.applymap(lambda x: round(x, 4)
  325. if isinstance(x, float) else x)
  327. cols_mdn_iqr_nrm = [c for c in group_stats_df.columns if c[0].startswith('Median')]
  328. group_stats_df_mdn_iqr_nrm = group_stats_df[cols_mdn_iqr_nrm]
  329. print("Median, IQR, and Normality test results:\n")
  330. print(group_stats_df_mdn_iqr_nrm.to_csv(sep=';'))
  332. cols_kruskal_ph = [c for c in group_stats_df.columns if not c[0].startswith('Median') and not c[0].startswith('Normality')]
  333. group_stats_df_kruskal_ph = group_stats_df[cols_kruskal_ph]
  334. print("Kruskal-Wallis and Conover post-hoc test results:\n")
  335. group_stats_df_kruskal_ph = group_stats_df_kruskal_ph.applymap(lambda x: round(x, 2) if isinstance(x, float) else x)
  336. print(group_stats_df_kruskal_ph.to_csv(sep=';'))
  338. export_significant_html_table(group_stats_df,
  339. os.path.join(results_dir, 'KW-Dunn_test.html'),
  340. hlines=[4,6,8,10,12],
  341. vlines=[])
  343. try:
  345. wilcoxon_stats_df = run_stat_tests_wilcoxon(dagg,
  346. split_by='Group',
  347. pair_by='Condition')
  348. print("\nWilcoxon Signed-Rank Test Results:\n")
  350. wilcoxon_stats_df.index.name = None
  351. wilcoxon_stats_df = wilcoxon_stats_df.applymap(lambda x: round(x, 2)
  352. if isinstance(x, float) else x)
  354. print(wilcoxon_stats_df.to_csv(sep=';'))
  356. except Exception as e:
  357. print(f"Error occurred while processing Wilcoxon test: {e}")
  358. print(traceback.format_exc())
  360. plot_layout = np.array([
  361. [('6mins', 'EBR (blink/min)'), ('6mins', 'LIBIV'), ('6mins', 'Mean Amplitude (µV)'),
  362. ],
  363. [ ('6mins', 'Mean Duration (ms)'), ('6mins', 'Mean Rise Time (ms)'), ('6mins', 'Mean Fall Time (ms)')
  364. ],
  365. ])
  367. # Plot the features summary for each group and condition
  368. plot_features_summary(dagg.query("Condition == 'Resting' or Condition == 'Oddball'"),
  369. plot_layout,
  370. palette=opt.palette)
  372. # hack for duration fall time and rise time axes
  373. plt.gcf().get_axes()[-1].sharey(plt.gcf().get_axes()[-2]) # rise and fall time share y with duration
  374. plt.gcf().set_size_inches(7.7, 5.5)
  375. [ax.set_title('') for ax in plt.gcf().get_axes()]
  377. # Add statistical annotations to the boxplots
  378. from statannotations.Annotator import Annotator
  379. # Define pairs for comparisons
  380. pairs = [
  381. (("HC", "Resting"), ("eMCS", "Resting")),
  382. (("HC", "Resting"), ("pDoC", "Resting")),
  383. (("eMCS", "Resting"), ("pDoC", "Resting")),
  384. (("HC", "Oddball"), ("eMCS", "Oddball")),
  385. (("HC", "Oddball"), ("pDoC", "Oddball")),
  386. (("eMCS", "Oddball"), ("pDoC", "Oddball")),
  387. ]
  389. # Add annotations to the features summary plot
  390. fig = plt.gcf()
  391. axes = fig.get_axes()
  392. for ax, (window, feature) in zip(axes, plot_layout.reshape(-1, 2), strict=True):
  393. if window is None or feature is None:
  394. continue
  395. # Filter data for the current feature and window
  396. data_to_annotate = dagg[['Group', 'Condition', feature]]
  397. data_to_annotate = data_to_annotate.rename(columns={feature: "Value"})
  399. # Initialize Annotator
  400. annotator = Annotator(
  401. ax=ax,
  402. pairs=pairs,
  403. data=data_to_annotate,
  404. x="Group",
  405. hue="Condition",
  406. y="Value",
  407. order=["HC", "eMCS", "pDoC"],
  408. hue_order=["Resting", "Oddball"],
  409. perform_stat_test=False, # We will set p-values manually
  410. verbose=False,
  411. )
  413. # Get post-hoc test results
  414. custom_pvalues = []
  415. for pair in pairs:
  416. group1, condition1 = pair[0]
  417. group2, condition2 = pair[1]
  418. feature_name = feature.split('[')[0].strip()
  420. # Extract the p-value
  421. try:
  422. p_value = group_stats_df.loc[(feature_name, condition1), (f"Conover", f"{group1} vs. {group2}")]
  423. if isinstance(p_value, str):
  424. match = re.findall(r'(0\.\d+)', p_value)
  425. p_value = float(match[1]) if match else np.nan# take the adjusted one
  427. else:
  428. p_value = np.nan # If the p-value is not found, set it to NaN
  429. except KeyError:
  430. p_value = np.nan # If the p-value is not found, set it to NaN
  432. if pd.isna(p_value) or p_value == '':
  433. p_value = 1
  434. # Add p-value to the list
  435. custom_pvalues.append(p_value)
  437. # Ensure the number of custom_pvalues matches the number of pairs
  438. if len(custom_pvalues) != len(pairs):
  439. raise ValueError("Mismatch between the number of pairs and custom p-values.")
  442. # Customize annotation format
  443. # Ensure the annotator object is properly initialized
  444. if not hasattr(annotator, "configure"):
  445. raise AttributeError("The annotator object is not properly initialized. Ensure it is created with a valid plotter.")
  447. # Configure the annotator with the correct parameters
  448. annotator.configure(
  449. text_format="star",
  450. loc="inside",
  451. hide_non_significant=True,
  452. )
  453. print(f"Adding annotations: { {pair: p
  454. for pair, p in zip(pairs, custom_pvalues, strict=True) } }")
  456. # Apply the p-values and annotate
  457. if custom_pvalues:
  458. annotator.set_pvalues_and_annotate(custom_pvalues)
  459. else:
  460. raise ValueError("Custom p-values are missing. Ensure they are calculated and passed correctly.")
  462. annotate_panels(plt.gcf().get_axes(),xy =(-0.35, 1.1))
  463. plt.tight_layout()
  464. plt.savefig(os.path.abspath(os.path.join(results_dir,'features_boxplot.svg')))
  465. plt.show(block=False)
  466. #'''
  469. # Logistic regression with age/sex covariates
  470. try:
  471. # Work on merged features + demographics
  472. df_lr = dagg_and_demo.copy()
  474. # Find and standardize age/sex columns
  475. def _find_col(df, names):
  476. for n in names:
  477. if n in df.columns:
  478. return n
  479. for c in df.columns:
  480. if c.lower() in {n.lower() for n in names}:
  481. return c
  482. return None
  484. age_src = _find_col(df_lr, ['age', 'Age', 'AGE'])
  485. sex_src = _find_col(df_lr, ['sex', 'Sex', 'SEX', 'Gender', 'gender'])
  487. if age_src is None or sex_src is None:
  488. print("Age/sex columns not found in demographics. Skipping logistic regression.")
  489. else:
  490. # Normalize age to numeric
  491. df_lr['age'] = pd.to_numeric(df_lr[age_src], errors='coerce')
  493. # Normalize sex to numeric (0/1), try common string encodings first
  494. if df_lr[sex_src].dtype == object:
  495. sex_map = {
  496. 'm': 1, 'male': 1, 'man': 1, 'maschio': 1,
  497. 'f': 0, 'female': 0, 'woman': 0, 'femmina': 0
  498. }
  499. sex_norm = (
  500. df_lr[sex_src].astype(str).str.strip().str.lower().map(sex_map)
  501. )
  502. # Fallback: attempt numeric coercion
  503. sex_norm = sex_norm.where(~sex_norm.isna(),
  504. pd.to_numeric(df_lr[sex_src], errors='coerce'))
  505. df_lr['sex'] = sex_norm
  506. else:
  507. df_lr['sex'] = pd.to_numeric(df_lr[sex_src], errors='coerce')
  509. # Ensure required columns
  510. features2screen = [
  511. 'EBR (blink/min)',
  512. 'LIBIV',
  513. 'Mean Amplitude (µV)',
  514. 'Mean Duration (ms)',
  515. 'Mean Rise Time (ms)',
  516. 'Mean Fall Time (ms)'
  517. ]
  518. features2screen = [f for f in features2screen if f in df_lr.columns]
  519. if not features2screen:
  520. print("No matching features found for logistic regression. Skipping.")
  521. else:
  522. # Run per condition
  523. def _run_condition(cond):
  524. sub = df_lr.query("Condition == @cond").copy()
  525. # Drop rows missing essential covariates or target
  526. sub = sub.dropna(subset=['age', 'sex', 'Group'])
  527. if sub.empty:
  528. return None
  529. dfs = []
  530. for feat in features2screen:
  531. sub_feat = sub.dropna(subset=[feat])
  532. if sub_feat.empty:
  533. continue
  534. try:
  535. res = logreg_covariate(sub_feat, feat, ['age', 'sex'], 'Group',
  536. verbose=True, plot=False)
  537. figro = plot_logreg_effects(res, feat, ['age','sex'],'Group', sub,
  538. group2col = {g:d['Resting']
  539. for g,d in opt.palette.items()},
  540. show_prob=False
  541. )
  542. figmp = plot_logreg_effects(res, feat, ['age','sex'],'Group', sub,
  543. group2col = {g:d['Resting']
  544. for g,d in opt.palette.items()},
  545. show_prob=True
  546. )
  547. os.makedirs(results_dir+'/logreg_plots/', exist_ok=True)
  548. # save relative odds fig
  549. featpath = feat.replace(" (blink/min)", "").replace(" ","_")
  550. figro.savefig(os.path.join(results_dir+'/logreg_plots/',
  551. f'logreg_{cond}_{featpath}_odds_ratio.svg'))
  552. # save mapped probabilities fig
  553. figmp.savefig(os.path.join(results_dir+'/logreg_plots/',
  554. f'logreg_{cond}_{featpath}_mapped_probabilities.svg'))
  555. #plt.show(block=True)
  556. # Attach feature name to index if needed
  557. if isinstance(res, pd.Series):
  558. res = res.to_frame().T
  559. res.index = pd.Index([feat], name='Feature')
  560. dfs.append(res)
  561. except Exception as _e:
  562. print(f"logreg failed for {cond} / {feat}: {_e}")
  563. if not dfs:
  564. return None
  565. out = pd.concat(dfs, axis=0, ignore_index=False)
  566. out['Condition'] = cond
  567. return out
  569. lr_age_df_resting = _run_condition('Resting')
  570. lr_age_df_oddball = _run_condition('Oddball')
  572. frames = [d for d in [lr_age_df_resting, lr_age_df_oddball] if d is not None]
  573. if not frames:
  574. print("No logistic regression output produced.")
  575. else:
  576. lr_age_df = pd.concat(frames, axis=0, ignore_index=False)
  578. # Formatting table
  579. req_cols = {
  580. 'feat_OR', 'feat_CI_lower', 'feat_CI_upper',
  581. 'feat_p-value', 'feat_R2_Nagelkerke', 'feat_AUC',
  582. 'feat_corrected_OR', 'feat_corrected_CI_lower', 'feat_corrected_CI_upper',
  583. 'feat_corrected_p-value', 'feat_corrected_R2_Nagelkerke', 'feat_corrected_AUC'
  584. }
  585. if req_cols.issubset(set(lr_age_df.columns)):
  586. lr_age_df_fmt = lr_age_df.round(4).assign(
  587. Uncorrected=lr_age_df['feat_OR'].round(4).astype(str) + " (" +
  588. lr_age_df['feat_CI_lower'].round(4).astype(str) + ", " +
  589. lr_age_df['feat_CI_upper'].round(4).astype(str) + ")",
  590. Corrected=lr_age_df['feat_corrected_OR'].round(4).astype(str) + " (" +
  591. lr_age_df['feat_corrected_CI_lower'].round(4).astype(str) + ", " +
  592. lr_age_df['feat_corrected_CI_upper'].round(4).astype(str) + ")"
  593. )
  594. lr_age_df_fmt = lr_age_df_fmt[['Condition',
  595. 'Uncorrected', 'feat_p-value',
  596. 'feat_R2_Nagelkerke', 'feat_AUC',
  597. 'Corrected', 'feat_corrected_p-value',
  598. 'feat_corrected_R2_Nagelkerke', 'feat_corrected_AUC']]
  600. # Prepare MultiIndex table
  601. lr_age_df_midx = lr_age_df_fmt.copy().reset_index(names=['Feature'])
  602. featorder = [feat.split(' [')[0] for feat in features2screen]
  603. # Enforce sorting by Feature (using features2screen order) and then by Condition
  604. lr_age_df_midx['Feature'] = pd.Categorical(lr_age_df_midx['Feature'],
  605. categories=featorder,
  606. ordered=True)
  607. lr_age_df_midx.sort_values(by=['Feature', 'Condition'],
  608. ascending=[True, False],
  609. inplace=True)
  610. lr_age_df_midx = lr_age_df_midx.set_index(
  611. pd.MultiIndex.from_arrays(
  612. [lr_age_df_midx['Feature'].values,
  613. lr_age_df_midx['Condition'].values],
  614. names=['Feature', 'Condition']
  615. )
  616. ).drop(columns=['Feature', 'Condition'])
  618. lr_age_df_midx.columns = pd.MultiIndex.from_tuples([
  619. ('Uncorrected', 'OR, median (CI 95%)'),
  620. ('Uncorrected', 'p'),
  621. ('Uncorrected', 'Nagelkerke R2'),
  622. ('Uncorrected', 'AUC'),
  623. ('Corrected', 'OR, median (CI 95%)'),
  624. ('Corrected', 'p'),
  625. ('Corrected', 'Nagelkerke R2'),
  626. ('Corrected', 'AUC')
  627. ])
  629. print("Logistic regression results:")
  630. print(lr_age_df_midx.to_csv(sep=';'))
  632. # Save
  633. out_csv = os.path.join(results_dir, 'logreg_covariate_results.csv')
  634. lr_age_df_midx.to_csv(out_csv)
  635. print(f"Saved: {out_csv}")
  636. else:
  637. # Fallback: save raw output
  638. print("Unexpected columns from logreg_covariate; saving raw output.")
  639. out_csv = os.path.join(results_dir, 'logreg_covariate_results_raw.csv')
  640. lr_age_df.to_csv(out_csv)
  641. print(f"Saved: {out_csv}")
  643. except Exception as e:
  644. print(f"Error during logistic regression analysis: {e}")
  645. print(traceback.format_exc())

feature_stats.py at commit e3894ca, under MIT · at the source

Overview

Authors: Alfonso Magliacano1, Leonardo Corsi1,2, Piergiuseppe Liuzzi1, Calogero Maria Oddo2, Andrea Mannini1, Anna Estraneo1
ORCID iDs: Alfonso Magliacano, Leonardo Corsi, Piergiuseppe Liuzzi, Calogero Maria Oddo, Andrea Mannini
  1. IRCCS Fondazione Don Carlo Gnocchi ETS, Via di Scandicci 269, 50143 Firenze, FI, Italy
  2. The BioRobotics Institute, Scuola Superiore Sant’Anna, Viale Rinaldo Piaggio, 34, 56025 Pontedera PI, Italy
Institutions: Scuola Superiore Sant'Anna (Italy)
Journal: Neuroscience of consciousness, volume 2026, issue 1, article niag043
Dates: received 25 November 2025; accepted 2 July 2026; published online 12 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1093/nc/niag043 · PMID 42592187 · PMCID PMC13464740 · OpenAlex W7202297577
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), other condition (population), clinical / translational (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning, Preprocessing, Smoothing, state filtering, decompositions, Physiology & signal measures
Keywords: severe acquired brain injury, disorders of consciousness, eye blink, electro-oculography, rehabilitation
Topic: Vestibular and auditory disorders (Neurology, Neuroscience), according to OpenAlex
Funding: Ministero della Salute (H43C21000180001)
Citations: not cited yet (Europe PMC); 41 references in the paper

Abstract

This study aimed to explore whether spontaneous eye blinking features from a brief vertical electro-oculogram recording can distinguish patients with severe acquired brain injury from healthy volunteers. Given evidence for shared neurobiological substrates in blinking control and attention, we expected altered blinking features in patients. In this observational cross-sectional study we extracted eye blink rate, amplitude, rise time, fall time, duration, and blink-blink interval variability from electro-oculography in 19 patients with severe acquired brain injury and 21 healthy controls during resting state and an auditory oddball task. We tested group differences with Kruskal-Wallis and Conover-Iman post-hoc comparisons. Differences covarying age and sex were assessed with Wald tests for odds ratios from logistic regression. Univariate Kruskal-Wallis testing showed longer rise time and duration and smaller amplitude in patients relative to healthy controls. When considering demographics as covariate, we show that amplitude is strongly confounded by demographics whereas rise time, duration and, to a lesser extent, fall time, are associated to the brain injury condition regardless of age or sex. Eye blink rate and blink-blink interval variability did not show statistically significant differences. Condition moderated effects: group contrasts were strongest at rest and attenuated during oddball. We found preliminary evidence that rise time and total duration are prolonged in patients independently of demographic factors, although with no difference between different levels of consciousness. Our results suggest a general slowing of blink dynamics consistent with altered basal ganglia dopaminergic activity, contributing physiological information relevant for monitoring and rehabilitation planning in brain injury.

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

Repository

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

Leonardo-Corsi/blink-waveforms-sabi

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e3894ca9960220cead95803de90fbb816992a7aa, 23 April 2026
Languages: Python (13)
Size: 114 files, 13 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, CITATION.cff, environment (pyproject.toml)
Not found: tests, continuous integration, documentation
Tools: pandas (9 files), NumPy (8 files), Matplotlib (6 files), MNE-Python (6 files), SciPy (4 files), seaborn (2 files), statsmodels (2 files), scikit-learn (1 file), scikit-posthocs (1 file), statannotations (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
15 files
At the source: github.com/Leonardo-Corsi/blink-waveforms-sabi

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:

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

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

Data

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

Data availability

The full data supporting this study are available with the code reproducing the article results at https://github.com/Leonardo-Corsi/blink-waveforms-sabi.

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

Versions

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

Version 3, 28 September 2026

  • Funding: added Ministero della Salute: H43C21000180001

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 5 keywords, 40 references.

Cite

This paper

Magliacano, A., Corsi, L., Liuzzi, P., Oddo, C. M., Mannini, A., & Estraneo, A. (2026). Demographics-robust spontaneous eye blinking slowing in patients with severe acquired brain injury. Neuroscience of consciousness, 2026(1), niag043. https://doi.org/10.1093/nc/niag043

BibTeX

@article{magliacano2026demographics,
author = {Magliacano, Alfonso and Corsi, Leonardo and Liuzzi, Piergiuseppe and Oddo, Calogero Maria and Mannini, Andrea and Estraneo, Anna},
title = {{Demographics-robust spontaneous eye blinking slowing in patients with severe acquired brain injury}},
journal = {Neuroscience of consciousness},
year = {2026},
month = aug,
volume = {2026},
number = {1},
pages = {niag043},
publisher = {Oxford University Press},
issn = {2057-2107},
doi = {10.1093/nc/niag043},
url = {https://doi.org/10.1093/nc/niag043},
pmid = {42592187},
pmcid = {PMC13464740}
}

RIS

TY - JOUR
AU - Magliacano, Alfonso
AU - Corsi, Leonardo
AU - Liuzzi, Piergiuseppe
AU - Oddo, Calogero Maria
AU - Mannini, Andrea
AU - Estraneo, Anna
TI - Demographics-robust spontaneous eye blinking slowing in patients with severe acquired brain injury
T2 - Neuroscience of consciousness
J2 - Neurosci Conscious
PY - 2026
DA - 2026/08/12
VL - 2026
IS - 1
SP - niag043
SN - 2057-2107
PB - Oxford University Press
DO - 10.1093/nc/niag043
UR - https://doi.org/10.1093/nc/niag043
LA - en
ER -

CSL-JSON

{
"id": "10.1093/nc/niag043",
"type": "article-journal",
"title": "Demographics-robust spontaneous eye blinking slowing in patients with severe acquired brain injury",
"container-title": "Neuroscience of consciousness",
"author": [
{
"family": "Magliacano",
"given": "Alfonso"
},
{
"family": "Corsi",
"given": "Leonardo"
},
{
"family": "Liuzzi",
"given": "Piergiuseppe"
},
{
"family": "Oddo",
"given": "Calogero Maria"
},
{
"family": "Mannini",
"given": "Andrea"
},
{
"family": "Estraneo",
"given": "Anna"
}
],
"container-title-short": "Neurosci Conscious",
"volume": "2026",
"issue": "1",
"page": "niag043",
"DOI": "10.1093/nc/niag043",
"PMID": "42592187",
"PMCID": "PMC13464740",
"ISSN": "2057-2107",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/nc/niag043",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
12
]
]
}
}

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.1093/brain/awaf412 [code]
Multimodal multicentre investigation of diagnostic and prognostic markers in disorders of consciousness.
Journal: Brain : a journal of neurology
In common: statannotations, MNE-Python, seaborn, 5 other tools, other condition, 3 references
[2] doi:10.1016/j.isci.2026.115329 [code]
Brain metastases converge on shared geometric architecture and transcriptomic landscape yet remain distinct from gliomas.
Journal: iScience
In common: scikit-posthocs, statannotations, statsmodels, 6 other tools, clinical / translational, other condition
[3] doi:10.1038/s41531-026-01380-1 [code]
Identifying maximal beta power from directional subthalamic local field potentials in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: statannotations, MNE-Python, statsmodels, 6 other tools, clinical / translational
[4] doi:10.1038/s41467-026-75662-w [code]
Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices.
Journal: Nature communications
In common: scikit-posthocs, MNE-Python, statsmodels, 6 other tools
[5] doi:10.1117/1.nph.13.3.035006 [code]
Characterizing developmental changes in infant habituation using functional change point detection.
Journal: Neurophotonics
In common: scikit-posthocs, statannotations, statsmodels, 5 other tools
[6] doi:10.1038/s43856-026-01518-5 [code]
Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states.
Journal: Communications medicine
In common: statannotations, statsmodels, seaborn, 5 other tools, clinical / translational, other condition
[7] doi:10.1371/journal.pone.0345651 [code]
Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition.
Journal: PloS one
In common: scikit-posthocs, statsmodels, seaborn, 5 other tools, clinical / translational
[8] doi:10.1162/imag.a.105 [code]
Right posterior theta reflects human parahippocampal phase resetting by salient cues during goal-directed navigation
Journal: n/a
In common: statannotations, MNE-Python, statsmodels, 5 other tools
[9] doi:10.1002/epi4.70336 [code]
Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II.
Journal: Epilepsia open
In common: statannotations, statsmodels, seaborn, 5 other tools, clinical / translational
[10] doi:10.1002/alz.71649 [code]
Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: statannotations, statsmodels, seaborn, 5 other tools, other condition

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.