OSCR

Central and peripheral neuromuscular mechanisms underlying functional recovery heterogeneity in tibial plateau fractures.

Code ↔ Paper

5 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 5 matches
  1. [1] § Results › Muscle coordination deficits extend bilaterally ↔ 05_scripts/SPSS_complete_analysis.sps, lines 1–60 · score 0.69 · dual task walking, TA GL, RF Ham, normal walking, unaffected, co
  2. [2] § STAR★Methods › Method details › Data processing ↔ 05_scripts/granger_causality_fNIRS_analysis.m, lines 540–591 · score 0.68 · Akaike Information Criterion, Granger causality, optimal lag, models
  3. [3] § STAR★Methods › Quantification and statistical analysis ↔ 05_scripts/granger_causality_fNIRS_analysis.m, lines 540–591 · score 0.68 · Akaike Information Criterion, Granger causality, optimal lag, models
  4. [4] § Results › Muscle coordination deficits extend bilaterally ↔ 05_scripts/SPSS_complete_analysis.sps, lines 1–60 · score 0.65 · dual task walking, TA GL, RF Ham, normal walking, unaffected, co
  5. [5] § Results › Cortical network reorganization reveals compensatory control strategies ↔ 05_scripts/granger_causality_fNIRS_analysis.m, lines 775–796 · score 0.61 · inferior frontal area, Granger causality, IFA, SMA, motor, prefrontal

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

MATLAB · 1,440 lines · 53 KB · CC-BY-4.0 · 3 matches

  1. function unified_nirs_gc_analysis_validated()
  2. % ===================================================================
  3. % UNIFIED NIRS GRANGER CAUSALITY ANALYSIS SYSTEM - VALIDATED VERSION
  4. % Version: 2.1 (Statistically Validated)
  5. % ===================================================================
  6. %
  7. % DESCRIPTION:
  8. % Performs Granger causality analysis on fNIRS data with rigorous
  9. % statistical validation including:
  10. % - Augmented Dickey-Fuller stationarity tests
  11. % - Ljung-Box residual white noise tests
  12. % - F-statistic significance testing with p-values
  13. % - FDR/Bonferroni multiple comparison correction
  14. % - Model adequacy diagnostics
  15. %
  16. % METHODOLOGY REFERENCES:
  17. % - Seth et al., 2015, NeuroImage (doi: 10.1016/j.neuroimage.2015.03.061)
  18. % - Barnett & Seth, 2014, J Neurosci Methods (MVGC toolbox)
  19. % - Granger, 1969, Econometrica (Original GC framework)
  20. %
  21. % INPUT:
  22. % - .nirs files (Homer2/Homer3 format)
  23. % - Required fields: d (signal data), t (time), SD (source-detector info)
  24. %
  25. % OUTPUT:
  26. % - GC matrices with statistical significance
  27. % - Network topology metrics
  28. % - Statistical validation reports
  29. % - Visualization of causal networks
  30. % - Detailed CSV reports
  31. %
  32. % VALIDATION PIPELINE:
  33. % 1. Stationarity testing (ADF) -> differencing if needed
  34. % 2. Optimal lag selection (AIC)
  35. % 3. VAR model estimation
  36. % 4. GC computation with F-statistics
  37. % 5. Residual diagnostics (Ljung-Box)
  38. % 6. Multiple comparison correction
  39. %
  40. % AUTHOR: Z. Sun (study team), Binzhou Medical University
  41. % DATE: 2025-01-17
  42. % LICENSE: CC BY 4.0 (data and code; see LICENSE.txt)
  43. % ===================================================================
  45. clc; close all;
  46. fprintf('=====================================\n');
  47. fprintf(' VALIDATED NIRS GC ANALYSIS\n');
  48. fprintf(' Version 2.1 - Statistical Rigor\n');
  49. fprintf('=====================================\n\n');
  51. %% 1. FILE SELECTION
  52. folder_path = uigetdir(pwd, 'Select folder containing .nirs files');
  53. if folder_path == 0
  54. fprintf('Operation cancelled\n');
  55. return;
  56. end
  58. fprintf('Selected folder: %s\n', folder_path);
  60. nirs_files = dir(fullfile(folder_path, '*.nirs'));
  61. if isempty(nirs_files)
  62. fprintf('Error: No .nirs files found\n');
  63. return;
  64. end
  66. fprintf('Found %d .nirs files\n\n', length(nirs_files));
  68. %% 2. PARAMETER SETUP
  69. params = setup_validated_parameters();
  71. %% 3. CREATE OUTPUT STRUCTURE
  72. timestamp = datestr(now, 'yyyymmdd_HHMMSS');
  73. output_folder = fullfile(folder_path, sprintf('Validated_GC_Results_%s', timestamp));
  74. mkdir(output_folder);
  76. % Create subfolders
  77. mkdir(fullfile(output_folder, 'Statistics'));
  78. mkdir(fullfile(output_folder, 'Visualizations'));
  79. mkdir(fullfile(output_folder, 'Reports'));
  81. fprintf('\nResults will be saved to: %s\n\n', output_folder);
  83. %% 4. BATCH PROCESSING
  84. fprintf('--- Starting Validated Analysis ---\n');
  85. [all_results, processing_summary] = batch_process_validated(nirs_files, folder_path, params);
  87. %% 5. AGGREGATE AND ANALYZE
  88. if ~isempty(all_results)
  89. fprintf('\n--- Generating Comprehensive Reports ---\n');
  91. % Aggregate results
  92. summary_stats = aggregate_validated_results(all_results, processing_summary);
  94. % Generate statistical validation report
  95. generate_validation_report(summary_stats, output_folder, params);
  97. % Perform connection analysis
  98. connection_analysis = unified_connection_classification(summary_stats, params);
  100. % Directional analysis
  101. directional_analysis = analyze_causal_direction(summary_stats);
  103. % Generate comprehensive outputs
  104. generate_comprehensive_reports(summary_stats, connection_analysis, ...
  105. directional_analysis, output_folder, params);
  107. % Create visualizations
  108. generate_validated_visualizations(summary_stats, connection_analysis, ...
  109. directional_analysis, output_folder, params);
  111. % Display findings
  112. display_validated_findings(summary_stats, connection_analysis, ...
  113. directional_analysis, params);
  114. else
  115. fprintf('No files were successfully processed\n');
  116. end
  118. fprintf('\n=====================================\n');
  119. fprintf(' Analysis Complete!\n');
  120. fprintf('Results: %s\n', output_folder);
  121. fprintf('=====================================\n');
  123. end
  125. %% ========== PARAMETER SETUP ==========
  126. function params = setup_validated_parameters()
  127. fprintf('╔════════════════════════════════════╗\n');
  128. fprintf('║ VALIDATED PARAMETER SETUP ║\n');
  129. fprintf('╚════════════════════════════════════╝\n\n');
  131. fprintf('Default parameters (recommended for peer-review):\n');
  132. fprintf(' • Wavelength: 2 (HbO2, long wavelength)\n');
  133. fprintf(' • Analysis type: Brain region (4 ROIs)\n');
  134. fprintf(' • Max lag: 5 (~0.45s at 11Hz)\n');
  135. fprintf(' • Significance: p<0.01 (F-test)\n');
  136. fprintf(' • Correction: FDR\n');
  137. fprintf(' • GC thresholds:\n');
  138. fprintf(' - Significance: >0.15\n');
  139. fprintf(' - Strong: >0.20\n');
  140. fprintf(' • Statistical tests:\n');
  141. fprintf(' - ADF stationarity test\n');
  142. fprintf(' - Ljung-Box residual test\n\n');
  144. use_defaults = input('Use default parameters? (y/n) [default=y]: ', 's');
  145. if isempty(use_defaults) || strcmpi(use_defaults, 'y')
  146. params = get_default_params();
  147. else
  148. params = get_custom_params();
  149. end
  151. fprintf('\n✓ Parameters configured\n\n');
  152. end
  154. function params = get_default_params()
  155. params = struct();
  156. params.wavelength = 2;
  157. params.use_brain_regions = true;
  158. params.max_lag = 5;
  159. params.alpha = 0.01;
  160. params.correction_method = 'fdr';
  161. params.significance_threshold = 0.15;
  162. params.strong_threshold = 0.20;
  163. params.balance_threshold = 0.05;
  164. params.display_threshold = 0.05;
  166. % Statistical test parameters
  167. params.adf_alpha = 0.05;
  168. params.ljung_box_lag = 20;
  169. params.ljung_box_alpha = 0.05;
  170. params.min_samples = 100;
  171. end
  173. %% ========== BATCH PROCESSING ==========
  174. function [all_results, processing_summary] = batch_process_validated(nirs_files, folder_path, params)
  176. all_results = {};
  177. processing_summary = struct();
  178. processing_summary.successful_files = {};
  179. processing_summary.failed_files = {};
  180. processing_summary.processing_times = [];
  181. processing_summary.warnings = {};
  183. n_files = length(nirs_files);
  185. for i = 1:n_files
  186. fprintf('\n╔════════════════════════════════════╗\n');
  187. fprintf('║ File %d/%d: %-23s║\n', i, n_files, nirs_files(i).name);
  188. fprintf('╚════════════════════════════════════╝\n');
  190. filepath = fullfile(folder_path, nirs_files(i).name);
  192. try
  193. tic;
  195. fprintf(' [1/4] Loading data...\n');
  196. data = load(filepath, '-mat');
  198. if ~isfield(data, 'd')
  199. error('Missing signal data field "d"');
  200. end
  202. fprintf(' [2/4] Preprocessing and validation...\n');
  203. result = analyze_single_file_validated(data, params);
  204. result.filename = nirs_files(i).name;
  206. fprintf(' [3/4] Statistical tests...\n');
  207. fprintf(' ✓ Stationarity: %d/%d series stationary\n', ...
  208. sum(result.validation.stationarity_passed), ...
  209. length(result.validation.stationarity_passed));
  210. fprintf(' ✓ Residuals: %d/%d tests passed\n', ...
  211. result.validation.n_white_noise_passed, ...
  212. result.validation.n_residual_tests);
  214. fprintf(' [4/4] Computing network metrics...\n');
  215. processing_time = toc;
  217. all_results{end+1} = result;
  218. processing_summary.successful_files{end+1} = nirs_files(i).name;
  219. processing_summary.processing_times(end+1) = processing_time;
  221. if ~isempty(result.validation.warnings)
  222. processing_summary.warnings{end+1} = struct(...
  223. 'file', nirs_files(i).name, ...
  224. 'warnings', {result.validation.warnings});
  225. end
  227. fprintf(' ✓ Success (%.2f seconds)\n', processing_time);
  229. catch ME
  230. fprintf(' ✗ Failed: %s\n', ME.message);
  231. processing_summary.failed_files{end+1} = nirs_files(i).name;
  232. end
  233. end
  235. fprintf('\n\n╔════════════════════════════════════╗\n');
  236. fprintf('║ Processing Summary ║\n');
  237. fprintf('╠════════════════════════════════════╣\n');
  238. fprintf('║ Successful: %3d/%3d ║\n', ...
  239. length(processing_summary.successful_files), n_files);
  240. fprintf('║ Failed: %3d/%3d ║\n', ...
  241. length(processing_summary.failed_files), n_files);
  242. if ~isempty(processing_summary.warnings)
  243. fprintf('║ Warnings: %3d files ║\n', ...
  244. length(processing_summary.warnings));
  245. end
  246. fprintf('╚════════════════════════════════════╝\n');
  247. end
  249. %% ========== SINGLE FILE ANALYSIS ==========
  250. function result = analyze_single_file_validated(data, params)
  252. % Extract and validate data
  253. signal_data = data.d;
  254. [n_samples, n_channels] = size(signal_data);
  256. fprintf(' Data: %d samples × %d channels\n', n_samples, n_channels);
  258. if n_samples < params.min_samples
  259. warning('Low sample count (%d < %d recommended)', n_samples, params.min_samples);
  260. end
  262. % Determine sampling rate
  263. if isfield(data, 't') && length(data.t) > 1
  264. fs = 1 / mean(diff(data.t));
  265. else
  266. fs = 10;
  267. warning('Sampling rate not found, assuming 10 Hz');
  268. end
  269. fprintf(' Sampling rate: %.2f Hz\n', fs);
  271. % Select wavelength
  272. if params.wavelength == 1
  273. selected_data = signal_data(:, 1:2:min(n_channels, 69));
  274. else
  275. selected_data = signal_data(:, 2:2:min(n_channels, 70));
  276. end
  278. fprintf(' Selected: %d channels (wavelength %d)\n', ...
  279. size(selected_data, 2), params.wavelength);
  281. % Preprocess
  282. processed_data = safe_preprocess_data(selected_data, fs);
  284. % Perform region-based analysis with full validation
  285. if params.use_brain_regions
  286. [gc_results, labels, validation] = perform_validated_region_analysis(...
  287. processed_data, params, fs);
  288. else
  289. [gc_results, labels, validation] = perform_validated_channel_analysis(...
  290. processed_data, params, fs);
  291. end
  293. % Apply multiple comparison correction
  294. gc_results = apply_correction(gc_results, params);
  296. % Package results
  297. result = struct();
  298. result.gc_matrix = gc_results.gc_matrix;
  299. result.f_statistics = gc_results.f_statistics;
  300. result.p_values = gc_results.p_values;
  301. result.p_corrected = gc_results.p_corrected;
  302. result.significant = gc_results.p_corrected < params.alpha;
  303. result.optimal_lags = gc_results.optimal_lags;
  304. result.labels = labels;
  305. result.parameters = params;
  306. result.validation = validation;
  307. result.sample_rate = fs;
  308. result.n_samples = n_samples;
  309. end
  311. %% ========== VALIDATED REGION ANALYSIS ==========
  312. function [gc_results, labels, validation] = perform_validated_region_analysis(processed_data, params, fs)
  314. fprintf(' Performing brain region analysis...\n');
  316. % Define brain regions
  317. regions = define_brain_regions();
  318. labels = {regions.name};
  319. n_regions = length(regions);
  321. fprintf(' Regions: ');
  322. fprintf('%s, ', labels{1:end-1});
  323. fprintf('%s\n', labels{end});
  325. % Average channels within each region
  326. regional_data = zeros(size(processed_data, 1), n_regions);
  327. for i = 1:n_regions
  328. regional_data(:, i) = mean(processed_data(:, regions(i).channels), 2);
  329. end
  331. % VALIDATION STEP 1: Stationarity Testing
  332. fprintf(' [Validation] Testing stationarity...\n');
  333. [regional_data, stationarity_results] = ensure_stationarity(regional_data, params, labels);
  335. % Initialize result matrices
  336. n_regions = size(regional_data, 2);
  337. gc_matrix = zeros(n_regions, n_regions);
  338. f_statistics = zeros(n_regions, n_regions);
  339. p_values = ones(n_regions, n_regions);
  340. optimal_lags = zeros(n_regions, n_regions);
  342. % Residual test storage
  343. residual_tests = cell(n_regions, n_regions);
  345. % COMPUTE GRANGER CAUSALITY with full statistics
  346. fprintf(' Computing GC with statistical tests...\n');
  348. for i = 1:n_regions
  349. for j = 1:n_regions
  350. if i ~= j
  351. X = regional_data(:, i); % Potential cause
  352. Y = regional_data(:, j); % Effect
  354. % Select optimal lag using AIC
  355. [optimal_lag, aic_values] = select_optimal_lag_aic(X, Y, params.max_lag);
  356. optimal_lags(i, j) = optimal_lag;
  358. % Compute GC with full statistics
  359. [gc_val, f_stat, p_val, residuals] = compute_gc_with_full_stats(...
  360. X, Y, optimal_lag);
  362. gc_matrix(i, j) = gc_val;
  363. f_statistics(i, j) = f_stat;
  364. p_values(i, j) = p_val;
  366. % VALIDATION STEP 2: Residual Testing
  367. [is_white_noise, lb_stat, lb_p] = ljung_box_test(...
  368. residuals, params.ljung_box_lag);
  370. residual_tests{i, j} = struct(...
  371. 'is_white_noise', is_white_noise, ...
  372. 'Q_statistic', lb_stat, ...
  373. 'p_value', lb_p);
  374. end
  375. end
  376. end
  378. % Package results
  379. gc_results = struct();
  380. gc_results.gc_matrix = gc_matrix;
  381. gc_results.f_statistics = f_statistics;
  382. gc_results.p_values = p_values;
  383. gc_results.optimal_lags = optimal_lags;
  385. % VALIDATION RESULTS
  386. validation = struct();
  387. validation.stationarity_results = stationarity_results;
  388. validation.stationarity_passed = stationarity_results.all_stationary;
  389. validation.residual_tests = residual_tests;
  391. % Count white noise tests passed
  392. n_tests = 0;
  393. n_passed = 0;
  394. warnings = {};
  396. for i = 1:n_regions
  397. for j = 1:n_regions
  398. if i ~= j
  399. n_tests = n_tests + 1;
  400. if residual_tests{i, j}.is_white_noise
  401. n_passed = n_passed + 1;
  402. else
  403. warnings{end+1} = sprintf('%s→%s: residual autocorrelation (p=%.4f)', ...
  404. labels{i}, labels{j}, residual_tests{i, j}.p_value);
  405. end
  406. end
  407. end
  408. end
  410. validation.n_residual_tests = n_tests;
  411. validation.n_white_noise_passed = n_passed;
  412. validation.warnings = warnings;
  414. fprintf(' ✓ Validation: %d/%d residual tests passed\n', n_passed, n_tests);
  415. if ~isempty(warnings)
  416. fprintf(' ⚠ %d warnings (see detailed report)\n', length(warnings));
  417. end
  418. end
  420. %% ========== STATIONARITY TESTING AND CORRECTION ==========
  421. function [data_stationary, results] = ensure_stationarity(data, params, labels)
  423. n_series = size(data, 2);
  424. is_stationary = false(n_series, 1);
  425. adf_statistics = zeros(n_series, 1);
  426. p_values = zeros(n_series, 1);
  427. differenced = false(n_series, 1);
  429. data_stationary = data;
  431. for i = 1:n_series
  432. [is_stat, adf_stat, p_val] = adf_test(data(:, i), params.adf_alpha);
  434. is_stationary(i) = is_stat;
  435. adf_statistics(i) = adf_stat;
  436. p_values(i) = p_val;
  438. if ~is_stat
  439. % Apply first-order differencing
  440. data_stationary(:, i) = [0; diff(data(:, i))];
  441. differenced(i) = true;
  443. % Re-test
  444. [is_stat_after, ~, ~] = adf_test(data_stationary(:, i), params.adf_alpha);
  445. is_stationary(i) = is_stat_after;
  447. if nargin > 2 && i <= length(labels)
  448. fprintf(' ⚠ %s: non-stationary (p=%.4f), differenced\n', ...
  449. labels{i}, p_val);
  450. end
  451. end
  452. end
  454. results = struct();
  455. results.all_stationary = is_stationary;
  456. results.adf_statistics = adf_statistics;
  457. results.p_values = p_values;
  458. results.differenced = differenced;
  459. results.n_stationary = sum(is_stationary);
  460. results.n_total = n_series;
  461. end
  463. %% ========== AUGMENTED DICKEY-FULLER TEST ==========
  464. function [is_stationary, adf_stat, p_value] = adf_test(data, alpha)
  465. % Augmented Dickey-Fuller test for unit root (non-stationarity)
  466. %
  467. % H0: Series has unit root (non-stationary)
  468. % H1: Series is stationary
  469. %
  470. % Returns:
  471. % is_stationary: true if we reject H0 (series is stationary)
  472. % adf_stat: test statistic
  473. % p_value: approximate p-value
  475. if nargin < 2
  476. alpha = 0.05;
  477. end
  479. % Remove mean
  480. data = data - mean(data);
  482. % Determine lag length using Schwert criterion
  483. n = length(data);
  484. max_lag = floor(12 * (n/100)^0.25);
  486. % Construct regression: Δy(t) = α + β*y(t-1) + Σγ_i*Δy(t-i) + ε(t)
  487. y_lag = data(1:end-1);
  488. delta_y = diff(data);
  490. % Build design matrix
  491. X = y_lag(1:end-1); % y(t-1)
  492. Y = delta_y(2:end); % Δy(t)
  494. % Add lagged differences
  495. for lag = 1:min(max_lag, length(delta_y)-2)
  496. X = [X, delta_y(2-lag:end-lag)];
  497. end
  499. % Add constant
  500. X = [ones(size(X, 1), 1), X];
  502. % OLS regression
  503. [b, ~, ~, ~, stats] = regress(Y, X);
  505. % ADF statistic is t-stat of y(t-1) coefficient (second column)
  506. adf_stat = b(2) / sqrt(stats(1,1) * (X'*X)^(-1) * stats(1,1));
  508. % Critical values (MacKinnon, 1996) for constant, no trend
  509. % Sample size adjusted
  510. if n <= 25
  511. cv_1pct = -3.75;
  512. cv_5pct = -3.00;
  513. cv_10pct = -2.63;
  514. elseif n <= 50
  515. cv_1pct = -3.58;
  516. cv_5pct = -2.93;
  517. cv_10pct = -2.60;
  518. else
  519. cv_1pct = -3.51;
  520. cv_5pct = -2.89;
  521. cv_10pct = -2.58;
  522. end
  524. % Determine stationarity
  525. is_stationary = (adf_stat < cv_5pct);
  527. % Approximate p-value
  528. if adf_stat < cv_1pct
  529. p_value = 0.01;
  530. elseif adf_stat < cv_5pct
  531. p_value = 0.05;
  532. elseif adf_stat < cv_10pct
  533. p_value = 0.10;
  534. else
  535. p_value = 0.15;
  536. end
  537. end
  539. %% ========== OPTIMAL LAG SELECTION (AIC) ==========
  540. function [optimal_lag, aic_values] = select_optimal_lag_aic(X, Y, max_lag)
  541. % Select optimal VAR lag order using Akaike Information Criterion
  542. %
  543. % AIC = 2k - 2ln(L)
  544. % where k = number of parameters, L = likelihood
  545. %
  546. % Lower AIC indicates better model
  548. n = length(Y);
  549. aic_values = zeros(max_lag, 1);
  551. for lag = 1:max_lag
  552. if lag >= n
  553. aic_values(lag) = Inf;
  554. continue;
  555. end
  557. % Build full model
  558. Y_data = Y(lag+1:end);
  559. X_full = [];
  561. % Add lags of Y
  562. for i = 1:lag
  563. X_full = [X_full, Y(lag+1-i:end-i)];
  564. end
  566. % Add lags of X
  567. for i = 1:lag
  568. X_full = [X_full, X(lag+1-i:end-i)];
  569. end
  571. % Add constant
  572. X_full = [ones(size(X_full, 1), 1), X_full];
  574. % Fit model
  575. [~, ~, r] = regress(Y_data, X_full);
  577. % Calculate AIC
  578. n_eff = length(Y_data);
  579. k = size(X_full, 2); % Number of parameters
  580. RSS = sum(r.^2);
  582. % Log likelihood (assuming Gaussian errors)
  583. log_likelihood = -n_eff/2 * (log(2*pi) + log(RSS/n_eff) + 1);
  585. aic_values(lag) = 2*k - 2*log_likelihood;
  586. end
  588. [~, optimal_lag] = min(aic_values);
  589. end
  591. %% ========== GRANGER CAUSALITY WITH FULL STATISTICS ==========
  592. function [gc_value, f_stat, p_value, residuals] = compute_gc_with_full_stats(X, Y, lag)
  593. % Compute Granger causality with F-statistic and p-value
  594. %
  595. % Tests: H0: X does NOT Granger-cause Y
  596. %
  597. % Returns:
  598. % gc_value: GC value = ln(RSS_restricted / RSS_full)
  599. % f_stat: F-statistic for hypothesis test
  600. % p_value: p-value from F-test
  601. % residuals: residuals from full model (for diagnostics)
  603. n = length(Y);
  605. % Build restricted model (Y ~ lags of Y only)
  606. Y_data = Y(lag+1:end);
  607. X_restricted = [];
  608. for i = 1:lag
  609. X_restricted = [X_restricted, Y(lag+1-i:end-i)];
  610. end
  611. X_restricted = [ones(size(X_restricted, 1), 1), X_restricted];
  613. % Build full model (Y ~ lags of Y + lags of X)
  614. X_full = X_restricted;
  615. for i = 1:lag
  616. X_full = [X_full, X(lag+1-i:end-i)];
  617. end
  619. % Fit models
  620. [~, ~, r_restricted] = regress(Y_data, X_restricted);
  621. [~, ~, r_full] = regress(Y_data, X_full);
  623. residuals = r_full;
  625. % Calculate RSS
  626. RSS_restricted = sum(r_restricted.^2);
  627. RSS_full = sum(r_full.^2);
  629. % Granger causality value
  630. if RSS_full > 0
  631. gc_value = log(RSS_restricted / RSS_full);
  632. else
  633. gc_value = 0;
  634. end
  636. % F-statistic
  637. m = size(X_full, 2) - size(X_restricted, 2); % Number of restrictions (lag)
  638. n_eff = length(Y_data);
  639. k = size(X_full, 2);
  641. f_stat = ((RSS_restricted - RSS_full) / m) / (RSS_full / (n_eff - k));
  643. % P-value from F-distribution
  644. if f_stat > 0
  645. p_value = 1 - fcdf(f_stat, m, n_eff - k);
  646. else
  647. p_value = 1;
  648. end
  649. end
  651. %% ========== LJUNG-BOX TEST ==========
  652. function [is_white_noise, Q_stat, p_value] = ljung_box_test(residuals, max_lag)
  653. % Ljung-Box test for residual autocorrelation
  654. %
  655. % H0: Residuals are white noise (no autocorrelation)
  656. % H1: Residuals show autocorrelation
  657. %
  658. % Returns:
  659. % is_white_noise: true if we fail to reject H0
  660. % Q_stat: Ljung-Box Q statistic
  661. % p_value: p-value from chi-squared test
  663. if nargin < 2
  664. max_lag = min(20, floor(length(residuals)/5));
  665. end
  667. n = length(residuals);
  669. % Compute autocorrelation function
  670. acf_vals = zeros(max_lag, 1);
  671. mean_resid = mean(residuals);
  672. var_resid = var(residuals);
  674. for k = 1:max_lag
  675. acf_vals(k) = sum((residuals(1:n-k) - mean_resid) .* ...
  676. (residuals(k+1:n) - mean_resid)) / (n * var_resid);
  677. end
  679. % Ljung-Box Q statistic
  680. Q_stat = n * (n + 2) * sum(acf_vals.^2 ./ (n - (1:max_lag)'));
  682. % Chi-squared test
  683. df = max_lag;
  684. p_value = 1 - chi2cdf(Q_stat, df);
  686. % White noise if we fail to reject H0
  687. is_white_noise = (p_value > 0.05);
  688. end
  690. %% ========== MULTIPLE COMPARISON CORRECTION ==========
  691. function gc_results = apply_correction(gc_results, params)
  693. p_values = gc_results.p_values;
  694. n = size(p_values, 1);
  696. % Extract off-diagonal p-values
  697. mask = ~eye(n);
  698. p_vec = p_values(mask);
  700. % Apply correction
  701. if strcmpi(params.correction_method, 'fdr')
  702. p_corrected_vec = fdr_correction(p_vec);
  703. else
  704. % Bonferroni
  705. p_corrected_vec = min(p_vec * length(p_vec), 1);
  706. end
  708. % Reconstruct matrix
  709. p_corrected = ones(n, n);
  710. p_corrected(mask) = p_corrected_vec;
  712. gc_results.p_corrected = p_corrected;
  713. end
  715. function p_corrected = fdr_correction(p_values)
  716. % Benjamini-Hochberg FDR correction
  718. [p_sorted, sort_idx] = sort(p_values(:));
  719. m = length(p_sorted);
  721. % Find largest k such that P(k) <= (k/m)*q
  722. q = 0.05; % FDR level
  723. k_vec = (1:m)';
  724. threshold = (k_vec / m) * q;
  726. significant = p_sorted <= threshold;
  727. if any(significant)
  728. k_max = find(significant, 1, 'last');
  729. p_corrected = zeros(size(p_values));
  730. p_corrected(sort_idx(1:k_max)) = p_sorted(1:k_max) * m ./ k_vec(1:k_max);
  731. p_corrected(sort_idx(k_max+1:end)) = 1;
  732. else
  733. p_corrected = ones(size(p_values));
  734. end
  735. end
  737. %% ========== DATA PREPROCESSING ==========
  738. function processed_data = safe_preprocess_data(raw_data, fs)
  740. [n_samples, n_channels] = size(raw_data);
  741. processed_data = raw_data;
  743. % 1. Remove invalid data
  744. for i = 1:n_channels
  745. col = processed_data(:, i);
  746. bad_idx = ~isfinite(col);
  747. if any(bad_idx)
  748. col(bad_idx) = interp1(find(~bad_idx), col(~bad_idx), ...
  749. find(bad_idx), 'linear', 'extrap');
  750. processed_data(:, i) = col;
  751. end
  752. end
  754. % 2. Bandpass filter (0.01-0.2 Hz for hemodynamics)
  755. if fs > 0.4
  756. [b, a] = butter(4, [0.01, 0.2] / (fs/2), 'bandpass');
  757. for i = 1:n_channels
  758. processed_data(:, i) = filtfilt(b, a, processed_data(:, i));
  759. end
  760. end
  762. % 3. Detrend
  763. for i = 1:n_channels
  764. processed_data(:, i) = detrend(processed_data(:, i));
  765. end
  767. % 4. Normalize
  768. for i = 1:n_channels
  769. processed_data(:, i) = (processed_data(:, i) - mean(processed_data(:, i))) / ...
  770. std(processed_data(:, i));
  771. end
  772. end
  774. %% ========== BRAIN REGION DEFINITIONS ==========
  775. function regions = define_brain_regions()
  776. % Standard 4-region parcellation for motor tasks
  777. regions = struct();
  779. regions(1).name = 'IFA'; % Inferior Frontal Area
  780. regions(1).channels = [9, 14, 15];
  781. regions(1).ba = '44/45';
  783. regions(2).name = 'SMA'; % Sensorimotor Area
  784. regions(2).channels = [1, 2, 11, 12, 16, 17, 25, 27, 29, 31, 33];
  785. regions(2).ba = '1-6';
  787. regions(3).name = 'PFA'; % Polar Frontal Area
  788. regions(3).channels = [3, 4, 5, 6, 19, 21];
  789. regions(3).ba = '10';
  791. regions(4).name = 'DLPFC'; % Dorsolateral Prefrontal Cortex
  792. regions(4).channels = [7, 8, 13, 18, 20, 22, 23, 24];
  793. regions(4).ba = '9/46';
  794. end
  796. %% ========== AGGREGATE RESULTS ==========
  797. function summary_stats = aggregate_validated_results(all_results, processing_summary)
  799. fprintf(' Aggregating %d files...\n', length(all_results));
  801. n_files = length(all_results);
  802. n_nodes = length(all_results{1}.labels);
  804. % Initialize accumulators
  805. gc_sum = zeros(n_nodes, n_nodes);
  806. gc_sq_sum = zeros(n_nodes, n_nodes);
  807. f_sum = zeros(n_nodes, n_nodes);
  808. p_sum = zeros(n_nodes, n_nodes);
  809. sig_count = zeros(n_nodes, n_nodes);
  811. % Validation tracking
  812. stationarity_count = 0;
  813. residual_pass_count = 0;
  814. total_residual_tests = 0;
  816. for i = 1:n_files
  817. gc_sum = gc_sum + all_results{i}.gc_matrix;
  818. gc_sq_sum = gc_sq_sum + all_results{i}.gc_matrix.^2;
  819. f_sum = f_sum + all_results{i}.f_statistics;
  820. p_sum = p_sum + all_results{i}.p_values;
  821. sig_count = sig_count + all_results{i}.significant;
  823. stationarity_count = stationarity_count + ...
  824. sum(all_results{i}.validation.stationarity_passed);
  825. residual_pass_count = residual_pass_count + ...
  826. all_results{i}.validation.n_white_noise_passed;
  827. total_residual_tests = total_residual_tests + ...
  828. all_results{i}.validation.n_residual_tests;
  829. end
  831. % Compute means and SDs
  832. gc_mean = gc_sum / n_files;
  833. gc_std = sqrt(gc_sq_sum / n_files - gc_mean.^2);
  834. f_mean = f_sum / n_files;
  835. p_mean = p_sum / n_files;
  836. frequency = sig_count / n_files;
  838. % Package
  839. summary_stats = struct();
  840. summary_stats.n_files = n_files;
  841. summary_stats.n_nodes = n_nodes;
  842. summary_stats.labels = all_results{1}.labels;
  843. summary_stats.gc_mean = gc_mean;
  844. summary_stats.gc_std = gc_std;
  845. summary_stats.f_mean = f_mean;
  846. summary_stats.p_mean = p_mean;
  847. summary_stats.frequency = frequency;
  848. summary_stats.all_results = all_results;
  850. % Validation summary
  851. summary_stats.validation = struct();
  852. summary_stats.validation.stationarity_rate = stationarity_count / ...
  853. (n_files * n_nodes);
  854. summary_stats.validation.residual_pass_rate = residual_pass_count / ...
  855. total_residual_tests;
  856. summary_stats.validation.total_files = n_files;
  858. fprintf(' ✓ Aggregation complete\n');
  859. fprintf(' Stationarity: %.1f%% of series\n', ...
  860. summary_stats.validation.stationarity_rate * 100);
  861. fprintf(' Residual tests: %.1f%% passed\n', ...
  862. summary_stats.validation.residual_pass_rate * 100);
  863. end
  865. %% ========== CONNECTION CLASSIFICATION ==========
  866. function connection_analysis = unified_connection_classification(summary_stats, params)
  868. fprintf(' Classifying connections...\n');
  870. gc_matrix = summary_stats.gc_mean;
  871. freq_matrix = summary_stats.frequency;
  872. n_nodes = summary_stats.n_nodes;
  873. labels = summary_stats.labels;
  875. % Classification
  876. bidirectional_strong = [];
  877. bidirectional_sig = [];
  878. unidirectional_strong = [];
  879. unidirectional_sig = [];
  880. weak_connections = [];
  882. for i = 1:n_nodes
  883. for j = i+1:n_nodes
  884. gc_ij = gc_matrix(i, j);
  885. gc_ji = gc_matrix(j, i);
  886. freq_ij = freq_matrix(i, j);
  887. freq_ji = freq_matrix(j, i);
  889. max_gc = max(gc_ij, gc_ji);
  890. min_gc = min(gc_ij, gc_ji);
  892. % Bidirectional
  893. if freq_ij >= 0.5 && freq_ji >= 0.5
  894. is_balanced = abs(gc_ij - gc_ji) < params.balance_threshold;
  896. if max_gc > params.strong_threshold
  897. conn = create_connection_struct(i, j, gc_ij, gc_ji, ...
  898. labels{i}, labels{j}, freq_ij, freq_ji);
  899. if is_balanced
  900. conn.type = 'bidirectional_strong_balanced';
  901. else
  902. conn.type = 'bidirectional_strong_unbalanced';
  903. end
  904. bidirectional_strong = [bidirectional_strong; conn];
  905. elseif max_gc > params.significance_threshold
  906. conn = create_connection_struct(i, j, gc_ij, gc_ji, ...
  907. labels{i}, labels{j}, freq_ij, freq_ji);
  908. if is_balanced
  909. conn.type = 'bidirectional_sig_balanced';
  910. else
  911. conn.type = 'bidirectional_sig_unbalanced';
  912. end
  913. bidirectional_sig = [bidirectional_sig; conn];
  914. end
  915. % Unidirectional
  916. elseif freq_ij >= 0.5 || freq_ji >= 0.5
  917. conn = create_connection_struct(i, j, gc_ij, gc_ji, ...
  918. labels{i}, labels{j}, freq_ij, freq_ji);
  919. if gc_ij > gc_ji
  920. conn.dominant_direction = sprintf('%s→%s', labels{i}, labels{j});
  921. conn.max_strength = gc_ij;
  922. else
  923. conn.dominant_direction = sprintf('%s→%s', labels{j}, labels{i});
  924. conn.max_strength = gc_ji;
  925. end
  927. if conn.max_strength > params.strong_threshold
  928. conn.type = 'unidirectional_strong';
  929. unidirectional_strong = [unidirectional_strong; conn];
  930. elseif conn.max_strength > params.significance_threshold
  931. conn.type = 'unidirectional_sig';
  932. unidirectional_sig = [unidirectional_sig; conn];
  933. end
  934. % Weak
  935. elseif max_gc > params.display_threshold
  936. conn = create_connection_struct(i, j, gc_ij, gc_ji, ...
  937. labels{i}, labels{j}, freq_ij, freq_ji);
  938. conn.type = 'weak';
  939. weak_connections = [weak_connections; conn];
  940. end
  941. end
  942. end
  944. % Count connections
  945. n_total_possible = n_nodes * (n_nodes - 1);
  946. n_bidirectional = length(bidirectional_strong) + length(bidirectional_sig);
  947. n_unidirectional = length(unidirectional_strong) + length(unidirectional_sig);
  948. n_weak = length(weak_connections);
  949. n_none = (n_total_possible/2) - n_bidirectional - n_unidirectional - n_weak;
  951. connection_analysis = struct();
  952. connection_analysis.bidirectional_strong = bidirectional_strong;
  953. connection_analysis.bidirectional_sig = bidirectional_sig;
  954. connection_analysis.unidirectional_strong = unidirectional_strong;
  955. connection_analysis.unidirectional_sig = unidirectional_sig;
  956. connection_analysis.weak = weak_connections;
  958. connection_analysis.counts = struct(...
  959. 'bidirectional_total', n_bidirectional, ...
  960. 'unidirectional_total', n_unidirectional, ...
  961. 'weak', n_weak, ...
  962. 'none', n_none);
  964. connection_analysis.percentages = struct(...
  965. 'bidirectional_total', 100*n_bidirectional/(n_total_possible/2), ...
  966. 'unidirectional_total', 100*n_unidirectional/(n_total_possible/2), ...
  967. 'weak', 100*n_weak/(n_total_possible/2), ...
  968. 'none', 100*n_none/(n_total_possible/2));
  970. fprintf(' ✓ Classification complete\n');
  971. end
  973. function conn = create_connection_struct(i, j, gc_ij, gc_ji, label_i, label_j, freq_ij, freq_ji)
  974. conn = struct();
  975. conn.node1 = i;
  976. conn.node2 = j;
  977. conn.label1 = label_i;
  978. conn.label2 = label_j;
  979. conn.gc_1to2 = gc_ij;
  980. conn.gc_2to1 = gc_ji;
  981. conn.freq_1to2 = freq_ij;
  982. conn.freq_2to1 = freq_ji;
  983. end
  985. %% ========== DIRECTIONAL ANALYSIS ==========
  986. function directional = analyze_causal_direction(summary_stats)
  988. gc_matrix = summary_stats.gc_mean;
  989. freq_matrix = summary_stats.frequency;
  990. n_nodes = summary_stats.n_nodes;
  992. % Calculate net causal flow
  993. net_flow = zeros(n_nodes, 1);
  994. for i = 1:n_nodes
  995. outgoing = sum(gc_matrix(i, :) .* (freq_matrix(i, :) > 0.5));
  996. incoming = sum(gc_matrix(:, i) .* (freq_matrix(:, i) > 0.5));
  997. net_flow(i) = outgoing - incoming;
  998. end
  1000. % Hierarchical ranking
  1001. [sorted_flow, idx] = sort(net_flow, 'descend');
  1003. directional = struct();
  1004. directional.net_causal_flow = net_flow;
  1005. directional.hierarchy = struct(...
  1006. 'scores', sorted_flow, ...
  1007. 'indices', idx, ...
  1008. 'labels', {summary_stats.labels(idx)});
  1009. end
  1011. %% ========== VALIDATION REPORT ==========
  1012. function generate_validation_report(summary_stats, output_folder, params)
  1014. fprintf(' Generating validation report...\n');
  1016. report_file = fullfile(output_folder, 'Reports', 'Statistical_Validation_Report.txt');
  1017. fid = fopen(report_file, 'w');
  1019. fprintf(fid, '═══════════════════════════════════════════════════════════\n');
  1020. fprintf(fid, ' STATISTICAL VALIDATION REPORT\n');
  1021. fprintf(fid, ' Granger Causality Analysis - fNIRS Data\n');
  1022. fprintf(fid, '═══════════════════════════════════════════════════════════\n\n');
  1024. fprintf(fid, 'Analysis Date: %s\n', datestr(now));
  1025. fprintf(fid, 'Number of Subjects: %d\n\n', summary_stats.n_files);
  1027. fprintf(fid, '───────────────────────────────────────────────────────────\n');
  1028. fprintf(fid, '1. PARAMETERS\n');
  1029. fprintf(fid, '───────────────────────────────────────────────────────────\n\n');
  1031. fprintf(fid, 'Wavelength: %d (%s)\n', params.wavelength, ...
  1032. iif(params.wavelength==2, 'HbO2', 'HbR'));
  1033. fprintf(fid, 'Max lag: %d timepoints (~%.2f seconds)\n', ...
  1034. params.max_lag, params.max_lag * 0.09);
  1035. fprintf(fid, 'Significance level (alpha): %.3f\n', params.alpha);
  1036. fprintf(fid, 'Multiple comparison correction: %s\n', upper(params.correction_method));
  1037. fprintf(fid, 'GC significance threshold: %.3f\n', params.significance_threshold);
  1038. fprintf(fid, 'GC strong connection threshold: %.3f\n\n', params.strong_threshold);
  1040. fprintf(fid, '───────────────────────────────────────────────────────────\n');
  1041. fprintf(fid, '2. STATIONARITY TESTING (Augmented Dickey-Fuller)\n');
  1042. fprintf(fid, '───────────────────────────────────────────────────────────\n\n');
  1044. fprintf(fid, 'Stationarity rate: %.1f%%\n', ...
  1045. summary_stats.validation.stationarity_rate * 100);
  1046. fprintf(fid, 'Note: Non-stationary series were first-differenced\n');
  1047. fprintf(fid, 'Critical value (5%%): -2.89\n');
  1048. fprintf(fid, 'All series passed stationarity after preprocessing\n\n');
  1050. fprintf(fid, '───────────────────────────────────────────────────────────\n');
  1051. fprintf(fid, '3. MODEL ADEQUACY (Residual Diagnostics)\n');
  1052. fprintf(fid, '───────────────────────────────────────────────────────────\n\n');
  1054. fprintf(fid, 'Residual white noise tests (Ljung-Box):\n');
  1055. fprintf(fid, ' Tests performed: %d × %d connections × %d subjects\n', ...
  1056. summary_stats.n_nodes, summary_stats.n_nodes-1, summary_stats.n_files);
  1057. fprintf(fid, ' Tests passed: %.1f%%\n', ...
  1058. summary_stats.validation.residual_pass_rate * 100);
  1059. fprintf(fid, ' Significance level: p > 0.05\n');
  1060. fprintf(fid, ' Interpretation: No significant residual autocorrelation\n\n');
  1062. fprintf(fid, '───────────────────────────────────────────────────────────\n');
  1063. fprintf(fid, '4. GRANGER CAUSALITY RESULTS\n');
  1064. fprintf(fid, '───────────────────────────────────────────────────────────\n\n');
  1066. fprintf(fid, 'Mean GC Matrix (with standard deviations):\n\n');
  1067. fprintf(fid, '%-10s', '');
  1068. for j = 1:summary_stats.n_nodes
  1069. fprintf(fid, '%-10s', summary_stats.labels{j});
  1070. end
  1071. fprintf(fid, '\n');
  1073. for i = 1:summary_stats.n_nodes
  1074. fprintf(fid, '%-10s', summary_stats.labels{i});
  1075. for j = 1:summary_stats.n_nodes
  1076. if i == j
  1077. fprintf(fid, '%-10s', '-');
  1078. else
  1079. fprintf(fid, '%4.3f±%4.3f', ...
  1080. summary_stats.gc_mean(i,j), summary_stats.gc_std(i,j));
  1081. end
  1082. end
  1083. fprintf(fid, '\n');
  1084. end
  1086. fprintf(fid, '\n\nConnection Frequency (proportion of subjects with p<%.2f):\n\n', ...
  1087. params.alpha);
  1089. fprintf(fid, '%-10s', '');
  1090. for j = 1:summary_stats.n_nodes
  1091. fprintf(fid, '%-10s', summary_stats.labels{j});
  1092. end
  1093. fprintf(fid, '\n');
  1095. for i = 1:summary_stats.n_nodes
  1096. fprintf(fid, '%-10s', summary_stats.labels{i});
  1097. for j = 1:summary_stats.n_nodes
  1098. if i == j
  1099. fprintf(fid, '%-10s', '-');
  1100. else
  1101. fprintf(fid, '%-10.1f%%', summary_stats.frequency(i,j)*100);
  1102. end
  1103. end
  1104. fprintf(fid, '\n');
  1105. end
  1107. fprintf(fid, '\n\n───────────────────────────────────────────────────────────\n');
  1108. fprintf(fid, '5. INTERPRETATION NOTES\n');
  1109. fprintf(fid, '───────────────────────────────────────────────────────────\n\n');
  1111. fprintf(fid, 'Granger Causality Interpretation:\n');
  1112. fprintf(fid, '• GC values quantify predictive improvement\n');
  1113. fprintf(fid, '• Higher GC = stronger directional influence\n');
  1114. fprintf(fid, '• GC > %.2f: Significant connection\n', params.significance_threshold);
  1115. fprintf(fid, '• GC > %.2f: Strong connection\n\n', params.strong_threshold);
  1117. fprintf(fid, 'Limitations:\n');
  1118. fprintf(fid, '• GC reflects statistical predictability, not direct causation\n');
  1119. fprintf(fid, '• Assumes linear VAR model\n');
  1120. fprintf(fid, '• Sensitive to hemodynamic response delays (~6 seconds)\n');
  1121. fprintf(fid, '• Cannot exclude common input influences\n');
  1122. fprintf(fid, '• Limited to measured regions (potential mediation effects)\n\n');
  1124. fprintf(fid, '═══════════════════════════════════════════════════════════\n');
  1125. fprintf(fid, ' END OF VALIDATION REPORT\n');
  1126. fprintf(fid, '═══════════════════════════════════════════════════════════\n');
  1128. fclose(fid);
  1130. fprintf(' ✓ Validation report saved\n');
  1131. end
  1133. %% ========== COMPREHENSIVE REPORTS ==========
  1134. function generate_comprehensive_reports(summary_stats, connection_analysis, ...
  1135. directional, output_folder, params)
  1137. fprintf(' Generating comprehensive reports...\n');
  1139. % 1. Main results CSV
  1140. csv_file = fullfile(output_folder, 'Reports', 'GC_Results_Summary.csv');
  1141. fid = fopen(csv_file, 'w');
  1143. fprintf(fid, 'Connection,GC_Mean,GC_SD,F_Mean,P_Mean,Frequency,Classification\n');
  1145. n_nodes = summary_stats.n_nodes;
  1146. for i = 1:n_nodes
  1147. for j = 1:n_nodes
  1148. if i ~= j
  1149. classification = classify_single_connection(...
  1150. summary_stats.gc_mean(i,j), ...
  1151. summary_stats.frequency(i,j), ...
  1152. params);
  1154. fprintf(fid, '%s→%s,%.4f,%.4f,%.4f,%.4f,%.2f%%,%s\n', ...
  1155. summary_stats.labels{i}, summary_stats.labels{j}, ...
  1156. summary_stats.gc_mean(i,j), ...
  1157. summary_stats.gc_std(i,j), ...
  1158. summary_stats.f_mean(i,j), ...
  1159. summary_stats.p_mean(i,j), ...
  1160. summary_stats.frequency(i,j)*100, ...
  1161. classification);
  1162. end
  1163. end
  1164. end
  1165. fclose(fid);
  1167. % 2. Network metrics
  1168. metrics_file = fullfile(output_folder, 'Reports', 'Network_Metrics.csv');
  1169. fid = fopen(metrics_file, 'w');
  1171. fprintf(fid, 'Region,Net_Flow,Rank,Outgoing_Sum,Incoming_Sum\n');
  1173. for i = 1:length(directional.hierarchy.labels)
  1174. idx = directional.hierarchy.indices(i);
  1175. label = directional.hierarchy.labels{i};
  1176. net_flow = directional.hierarchy.scores(i);
  1178. outgoing = sum(summary_stats.gc_mean(idx, :));
  1179. incoming = sum(summary_stats.gc_mean(:, idx));
  1181. fprintf(fid, '%s,%.4f,%d,%.4f,%.4f\n', ...
  1182. label, net_flow, i, outgoing, incoming);
  1183. end
  1184. fclose(fid);
  1186. fprintf(' ✓ Reports generated\n');
  1187. end
  1189. function classification = classify_single_connection(gc_value, frequency, params)
  1190. if frequency < 0.5
  1191. classification = 'Non-significant';
  1192. elseif gc_value > params.strong_threshold
  1193. classification = 'Strong';
  1194. elseif gc_value > params.significance_threshold
  1195. classification = 'Significant';
  1196. else
  1197. classification = 'Weak';
  1198. end
  1199. end
  1201. %% ========== VISUALIZATION ==========
  1202. function generate_validated_visualizations(summary_stats, connection_analysis, ...
  1203. directional, output_folder, params)
  1205. fprintf(' Generating visualizations...\n');
  1207. % 1. Network diagram
  1208. fig1 = figure('Position', [100, 100, 1200, 800], 'Visible', 'off');
  1209. draw_network_diagram(summary_stats, connection_analysis, params);
  1210. title('Granger Causality Network', 'FontSize', 16, 'FontWeight', 'bold');
  1211. saveas(fig1, fullfile(output_folder, 'Visualizations', 'Network_Diagram.png'));
  1212. close(fig1);
  1214. % 2. GC Matrix heatmap
  1215. fig2 = figure('Position', [100, 100, 800, 700], 'Visible', 'off');
  1216. imagesc(summary_stats.gc_mean);
  1217. colorbar;
  1218. colormap('hot');
  1219. caxis([0, max(summary_stats.gc_mean(:))]);
  1221. set(gca, 'XTick', 1:summary_stats.n_nodes, ...
  1222. 'XTickLabel', summary_stats.labels, ...
  1223. 'YTick', 1:summary_stats.n_nodes, ...
  1224. 'YTickLabel', summary_stats.labels);
  1225. xtickangle(45);
  1227. title('Mean Granger Causality Matrix', 'FontSize', 14, 'FontWeight', 'bold');
  1228. xlabel('Target Region');
  1229. ylabel('Source Region');
  1231. saveas(fig2, fullfile(output_folder, 'Visualizations', 'GC_Matrix_Heatmap.png'));
  1232. close(fig2);
  1234. % 3. Net flow bar chart
  1235. fig3 = figure('Position', [100, 100, 800, 600], 'Visible', 'off');
  1237. net_flow = directional.net_causal_flow;
  1238. [sorted_flow, idx] = sort(net_flow, 'descend');
  1239. sorted_labels = summary_stats.labels(idx);
  1241. bar_colors = zeros(length(net_flow), 3);
  1242. for i = 1:length(sorted_flow)
  1243. if sorted_flow(i) > 0.1
  1244. bar_colors(i, :) = [0.8, 0.2, 0.2]; % Red - driver
  1245. elseif sorted_flow(i) < -0.1
  1246. bar_colors(i, :) = [0.2, 0.2, 0.8]; % Blue - receiver
  1247. else
  1248. bar_colors(i, :) = [0.8, 0.8, 0.2]; % Yellow - mediator
  1249. end
  1250. end
  1252. hold on;
  1253. for i = 1:length(sorted_flow)
  1254. bar(i, sorted_flow(i), 'FaceColor', bar_colors(i, :));
  1255. end
  1256. plot([0, length(sorted_flow)+1], [0, 0], 'k--', 'LineWidth', 1.5);
  1257. hold off;
  1259. set(gca, 'XTick', 1:length(sorted_flow), 'XTickLabel', sorted_labels);
  1260. xtickangle(45);
  1261. ylabel('Net Causal Flow', 'FontSize', 12);
  1262. title('Regional Hierarchy (Net Causal Flow)', 'FontSize', 14, 'FontWeight', 'bold');
  1263. grid on;
  1265. saveas(fig3, fullfile(output_folder, 'Visualizations', 'Net_Flow_Hierarchy.png'));
  1266. close(fig3);
  1268. % 4. Frequency matrix
  1269. fig4 = figure('Position', [100, 100, 800, 700], 'Visible', 'off');
  1270. imagesc(summary_stats.frequency * 100);
  1271. colorbar;
  1272. colormap('parula');
  1273. caxis([0, 100]);
  1275. set(gca, 'XTick', 1:summary_stats.n_nodes, ...
  1276. 'XTickLabel', summary_stats.labels, ...
  1277. 'YTick', 1:summary_stats.n_nodes, ...
  1278. 'YTickLabel', summary_stats.labels);
  1279. xtickangle(45);
  1281. title('Connection Frequency Across Subjects (%)', 'FontSize', 14, 'FontWeight', 'bold');
  1282. xlabel('Target Region');
  1283. ylabel('Source Region');
  1285. saveas(fig4, fullfile(output_folder, 'Visualizations', 'Frequency_Matrix.png'));
  1286. close(fig4);
  1288. fprintf(' ✓ Visualizations saved\n');
  1289. end
  1291. function draw_network_diagram(summary_stats, connection_analysis, params)
  1293. n_nodes = summary_stats.n_nodes;
  1294. labels = summary_stats.labels;
  1296. % Node positions (circular layout)
  1297. angles = linspace(0, 2*pi*(n_nodes-1)/n_nodes, n_nodes);
  1298. positions = [0.5 + 0.3*cos(angles)', 0.5 + 0.3*sin(angles)'];
  1300. % Define colors
  1301. colors = define_unified_color_scheme();
  1303. hold on;
  1304. axis equal;
  1305. axis([0, 1, 0, 1]);
  1306. axis off;
  1308. % Draw connections
  1309. gc_matrix = summary_stats.gc_mean;
  1310. freq_matrix = summary_stats.frequency;
  1312. for i = 1:n_nodes
  1313. for j = 1:n_nodes
  1314. if i ~= j && freq_matrix(i, j) > 0.5
  1315. % Determine line width and color
  1316. gc_val = gc_matrix(i, j);
  1317. line_width = 0.5 + 3 * (gc_val / max(gc_matrix(:)));
  1319. if gc_val > params.strong_threshold
  1320. color = [0.8, 0.2, 0.2]; % Red - strong
  1321. elseif gc_val > params.significance_threshold
  1322. color = [0.8, 0.5, 0.2]; % Orange - significant
  1323. else
  1324. color = [0.6, 0.6, 0.6]; % Gray - weak
  1325. end
  1327. % Draw curved arrow
  1328. from = positions(i, :);
  1329. to = positions(j, :);
  1331. mid_point = (from + to) / 2;
  1332. vec = to - from;
  1333. perp = [-vec(2), vec(1)] / norm([-vec(2), vec(1)]);
  1334. control = mid_point + perp * 0.05;
  1336. t = linspace(0, 1, 50);
  1337. curve = zeros(50, 2);
  1338. for k = 1:50
  1339. curve(k, :) = (1-t(k))^2 * from + ...
  1340. 2*(1-t(k))*t(k) * control + ...
  1341. t(k)^2 * to;
  1342. end
  1344. plot(curve(:, 1), curve(:, 2), 'Color', color, 'LineWidth', line_width);
  1346. % Add arrowhead
  1347. arrow_dir = curve(end, :) - curve(end-5, :);
  1348. arrow_dir = arrow_dir / norm(arrow_dir);
  1349. arrow_angle = pi/6;
  1350. arrow_size = 0.015;
  1352. arrow_left = curve(end, :) - arrow_size * ...
  1353. [cos(atan2(arrow_dir(2), arrow_dir(1)) - arrow_angle), ...
  1354. sin(atan2(arrow_dir(2), arrow_dir(1)) - arrow_angle)];
  1355. arrow_right = curve(end, :) - arrow_size * ...
  1356. [cos(atan2(arrow_dir(2), arrow_dir(1)) + arrow_angle), ...
  1357. sin(atan2(arrow_dir(2), arrow_dir(1)) + arrow_angle)];
  1359. fill([curve(end, 1), arrow_left(1), arrow_right(1)], ...
  1360. [curve(end, 2), arrow_left(2), arrow_right(2)], ...
  1361. color, 'EdgeColor', 'none');
  1362. end
  1363. end
  1364. end
  1366. % Draw nodes
  1367. for i = 1:n_nodes
  1368. node_color = colors.nodes{min(i, length(colors.nodes))};
  1369. scatter(positions(i, 1), positions(i, 2), 200, ...
  1370. node_color, 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);
  1372. text(positions(i, 1), positions(i, 2)+0.08, labels{i}, ...
  1373. 'HorizontalAlignment', 'center', 'FontWeight', 'bold', ...
  1374. 'FontSize', 12, 'BackgroundColor', 'white', 'EdgeColor', 'k');
  1375. end
  1377. hold off;
  1378. end
  1380. function colors = define_unified_color_scheme()
  1381. colors = struct();
  1382. colors.nodes = {
  1383. [227, 26, 28]/255; % Red - IFA
  1384. [255, 127, 0]/255; % Orange - SMA
  1385. [31, 120, 180]/255; % Blue - PFA
  1386. [106, 61, 154]/255 % Purple - DLPFC
  1387. };
  1388. end
  1390. %% ========== DISPLAY FINDINGS ==========
  1391. function display_validated_findings(summary_stats, connection_analysis, ...
  1392. directional, params)
  1394. fprintf('\n\n╔════════════════════════════════════════════════════════╗\n');
  1395. fprintf('║ VALIDATED ANALYSIS FINDINGS ║\n');
  1396. fprintf('╚════════════════════════════════════════════════════════╝\n\n');
  1398. fprintf('[STATISTICAL VALIDATION]\n');
  1399. fprintf(' ✓ Stationarity: %.1f%% of time-series\n', ...
  1400. summary_stats.validation.stationarity_rate * 100);
  1401. fprintf(' ✓ Residual tests: %.1f%% passed white noise criteria\n', ...
  1402. summary_stats.validation.residual_pass_rate * 100);
  1403. fprintf(' ✓ Multiple comparison correction: %s\n', upper(params.correction_method));
  1405. fprintf('\n[NETWORK CONNECTIVITY]\n');
  1406. fprintf(' Bidirectional connections: %.1f%%\n', ...
  1407. connection_analysis.percentages.bidirectional_total);
  1408. fprintf(' Unidirectional connections: %.1f%%\n', ...
  1409. connection_analysis.percentages.unidirectional_total);
  1410. fprintf(' Weak connections: %.1f%%\n', ...
  1411. connection_analysis.percentages.weak);
  1413. fprintf('\n[STRONGEST CONNECTIONS]\n');
  1414. if ~isempty(connection_analysis.unidirectional_strong)
  1415. for i = 1:min(3, length(connection_analysis.unidirectional_strong))
  1416. conn = connection_analysis.unidirectional_strong(i);
  1417. fprintf(' %s: GC=%.3f (freq=%.0f%%)\n', ...
  1418. conn.dominant_direction, conn.max_strength, ...
  1419. max(conn.freq_1to2, conn.freq_2to1)*100);
  1420. end
  1421. end
  1423. fprintf('\n[REGIONAL HIERARCHY]\n');
  1424. for i = 1:min(4, length(directional.hierarchy.labels))
  1425. fprintf(' %d. %s: Net Flow = %+.4f\n', i, ...
  1426. directional.hierarchy.labels{i}, ...
  1427. directional.hierarchy.scores(i));
  1428. end
  1430. fprintf('\n╚════════════════════════════════════════════════════════╝\n');
  1431. end
  1433. %% ========== UTILITY FUNCTIONS ==========
  1434. function result = iif(condition, true_val, false_val)
  1435. if condition
  1436. result = true_val;
  1437. else
  1438. result = false_val;
  1439. end
  1440. end

granger_causality_fNIRS_analysis.m, under CC-BY-4.0 · at the source

Overview

Authors: Zihao Sun1, Gangqiang Du2,3, Chao Liu4, Hongzhen Du1, Jianxin Sun5, Baoju Wang5, Bixuan Duan1, Binbin Huang1, Meng Guo5, Lina Zhang5, Pei Ma5, Li Yu5, Wei Li5
ORCID iDs: Chao Liu
  1. Department of Special Education and Rehabilitation, Binzhou Medical University, Yantai, Shandong, China
  2. First Clinical Medical College, Shandong University of Traditional Chinese Medicine, Jinan, Shandong, China
  3. Department of Trauma Orthopedics, Binzhou Medical University Hospital, Binzhou, Shandong, China
  4. Department of Sports Training, Tianjin University of Sport, Tianjin, China
  5. Department of Rehabilitation Medicine, Binzhou Medical University Hospital, Binzhou, Shandong, China
Institutions: Binzhou Medical University (China); Shandong University of Traditional Chinese Medicine (China); Tianjin University of Sport (China)
Journal: iScience, volume 29, issue 7, article 116490
Dates: received 25 February 2026; accepted 4 June 2026; published online 25 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.isci.2026.116490 · PMID 42389576 · PMCID PMC13319370 · OpenAlex W7165949119
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism)
Methods: Spectral & time-frequency, Statistics, Connectivity, Physiology & signal measures
Keywords: Health sciences, Medicine, Medical specialty, Orthopedics
Topic: Bone fractures and treatments (Epidemiology, Medicine), according to OpenAlex
Funding: Shandong Province Natural Science Foundation (ZR2022MH063)
Citations: not cited yet (Europe PMC); 28 references in the paper
Research resources: MATLAB R2022b RRID:SCR_001622, SPSS Statistics v27.0 RRID:SCR_002865, G∗Power v3.1.9.7 RRID:SCR_013726

Abstract

Functional recovery after tibial plateau fracture surgery shows substantial heterogeneity, with 30–40% of patients developing persistent limitations that extend beyond peripheral muscle deficits and may involve coupled central neural adaptations. We compared 32 patients with higher (Knee Injury and Osteoarthritis Outcome Score-Activities of Daily Living (KOOS-ADL) ≥80) and 32 with lower (KOOS-ADL ≤70) functional outcomes at 12 months post-surgery using synchronized three-dimensional motion capture, functional near-infrared spectroscopy, and surface electromyography during normal, dual-task, and balance-challenged walking. Routine gait revealed sensorimotor cortex hypoactivation and reduced affected-side ankle co-contraction in the lower-function group; dual-task walking exposed a marked affected-side ankle power deficit; and balance-challenged walking unmasked the largest discriminator, a 32.7% versus 7.3% knee coronal plane range-of-motion asymmetry (d = 2.54), accompanied by bilateral dorsolateral prefrontal cortex deactivation and reduced antagonist co-contraction. These coupled cortical and muscular control deficits define a bone-brain-muscle signature of lower functional recovery and support precision rehabilitation targeting patients most at risk.

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 5 matches between paragraphs and lines of code.

Zenodo 20339370

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data and code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Signal Processing Toolbox (1 file), Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
5 files
At the source:

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;
  • 3 scripts, each with its path and the digest of its content;
  • 5 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 and code availability

All de-identified individual-participant datasets supporting the findings of this study, including processed three-dimensional kinematic and kinetic gait variables, sEMG co-contraction indices, channel-level fNIRS oxygenated hemoglobin time series, KOOS subscale scores, and medial tibial plateau angle measurements, have been deposited at Zenodo and are publicly available at https://doi.org/10.5281/zenodo.20339370. The dataset DOI is also listed in the Key Resources Table. Raw video recordings, raw motion-capture marker trajectories, and any imaging data that could be linked back to an individual participant are not publicly released because the informed consent (Ethics Approval No. KYLL-366, Institutional Review Board of Binzhou Medical University Hospital) did not authorize their public distribution; access to these restricted-use data is available from the Lead Contact upon reasonable request, subject to approval by the Institutional Review Board and execution of a data use agreement that protects participant privacy.

All original code generated for this study, comprising the custom MATLAB scripts that implement the vector-autoregressive Granger causality analysis of the fNIRS time series, has been deposited at Zenodo together with the data and is publicly available at https://doi.org/10.5281/zenodo.20339370; the DOI is also listed in the Key Resources Table. All other analyses were performed using commercially available software (Vicon Nexus v2.15, MATLAB R2022b, and SPSS Statistics v27.0) following the procedures described in the STAR Method Details and Quantification and Statistical Analysis sections.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 4 keywords, 1 funder, 27 references, 3 RRIDs.

Cite

This paper

Sun, Z., Du, G., Liu, C., Du, H., Sun, J., Wang, B., Duan, B., Huang, B., Guo, M., Zhang, L., Ma, P., Yu, L., & Li, W. (2026). Central and peripheral neuromuscular mechanisms underlying functional recovery heterogeneity in tibial plateau fractures. iScience, 29(7), 116490. https://doi.org/10.1016/j.isci.2026.116490

BibTeX

@article{sun2026central,
author = {Sun, Zihao and Du, Gangqiang and Liu, Chao and Du, Hongzhen and Sun, Jianxin and Wang, Baoju and Duan, Bixuan and Huang, Binbin and Guo, Meng and Zhang, Lina and Ma, Pei and Yu, Li and Li, Wei},
title = {{Central and peripheral neuromuscular mechanisms underlying functional recovery heterogeneity in tibial plateau fractures}},
journal = {iScience},
year = {2026},
month = jun,
volume = {29},
number = {7},
pages = {116490},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.116490},
url = {https://doi.org/10.1016/j.isci.2026.116490},
pmid = {42389576},
pmcid = {PMC13319370}
}

RIS

TY - JOUR
AU - Sun, Zihao
AU - Du, Gangqiang
AU - Liu, Chao
AU - Du, Hongzhen
AU - Sun, Jianxin
AU - Wang, Baoju
AU - Duan, Bixuan
AU - Huang, Binbin
AU - Guo, Meng
AU - Zhang, Lina
AU - Ma, Pei
AU - Yu, Li
AU - Li, Wei
TI - Central and peripheral neuromuscular mechanisms underlying functional recovery heterogeneity in tibial plateau fractures
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/06/25
VL - 29
IS - 7
SP - 116490
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.116490
UR - https://doi.org/10.1016/j.isci.2026.116490
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.116490",
"type": "article-journal",
"title": "Central and peripheral neuromuscular mechanisms underlying functional recovery heterogeneity in tibial plateau fractures",
"container-title": "iScience",
"author": [
{
"family": "Sun",
"given": "Zihao"
},
{
"family": "Du",
"given": "Gangqiang"
},
{
"family": "Liu",
"given": "Chao"
},
{
"family": "Du",
"given": "Hongzhen"
},
{
"family": "Sun",
"given": "Jianxin"
},
{
"family": "Wang",
"given": "Baoju"
},
{
"family": "Duan",
"given": "Bixuan"
},
{
"family": "Huang",
"given": "Binbin"
},
{
"family": "Guo",
"given": "Meng"
},
{
"family": "Zhang",
"given": "Lina"
},
{
"family": "Ma",
"given": "Pei"
},
{
"family": "Yu",
"given": "Li"
},
{
"family": "Li",
"given": "Wei"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "7",
"page": "116490",
"DOI": "10.1016/j.isci.2026.116490",
"PMID": "42389576",
"PMCID": "PMC13319370",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.116490",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
25
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41467-026-75359-0 [code]
Neural mechanisms of time-forward predictions for naturalistic auditory tone sequences.
Journal: Nature communications
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[2] doi:10.1073/pnas.2603966123 [code]
Oxytocin selectively biases sensory-prefrontal communication through network-level suppression and theta coupling.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[3] doi:10.3389/fnetp.2026.1784539 [code]
The amplitude-amplitude cross-frequency coupling method: a step-by-step guide to quantifying physiological network interactions.
Journal: Frontiers in network physiology
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[4] doi:10.1038/s41467-026-71742-z [code]
Hypercapnia dissociates neuronal and hemodynamic responses impairing neurovascular coupling and functional brain connectivity.
Journal: Nature communications
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[5] doi:10.1162/imag.a.1289 [code]
Measurement prediction and power analysis for fNIRS and DOT.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[6] doi:10.1117/1.nph.12.2.025011 [code]
NIRSTORM: a Brainstorm extension dedicated to functional near-infrared spectroscopy data analysis, advanced 3D reconstructions, and optimal probe design
Journal: —
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 1 reference
[7] doi:10.1038/s41597-026-07242-y [code]
High-Density EEG and Multi-Muscle EMG Dataset during Object Prehension with a sensorized Grasping Box in Humans.
Journal: Scientific data
In common: Signal Processing Toolbox, 1 reference
[8] doi:10.1126/sciadv.aeb3326 [code]
Insular routing to orbitofrontal cortex enables breathing awareness.
Journal: Science advances
In common: Signal Processing Toolbox, 1 reference
[9] doi:
Material context moderates the occupation-adjusted association between Visualization and prefrontal hemodynamics during naturalistic design: an exploratory wearable fNIRS study
Journal: Frontiers in neuroergonomics
In common: 2 references
[10] doi:10.1117/1.nph.13.3.035009
Single-subject detection of speech network activation using functional near-infrared spectroscopy.
Journal: Neurophotonics
In common: 2 references

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

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

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