diff --git a/demos/efficiency/demo_efficiency.m b/demos/efficiency/demo_efficiency.m new file mode 100644 index 000000000..d506c3830 --- /dev/null +++ b/demos/efficiency/demo_efficiency.m @@ -0,0 +1,106 @@ +% (C) Copyright 2023 bidspm developers + +% Example to compute efficiency from +% https://www.nature.com/articles/s41598-024-52967-8 +% +% Related to this thread: +% https://twitter.com/m_wall/status/1757311725777387739 + +clear; +clc; +close all; + +addpath(fullfile(pwd, '..', '..')); +bidspm(); + +conditions = {'alcohol', 'neutral', 'negative', 'positive'}; +for i = 1:numel(conditions) + design_matrix{i} = ['trial_type.' conditions{i}]; +end + +%% create stats model JSON +bm = BidsModel(); +NodeIdx = 1; + +% high pass filter +bm.Nodes{NodeIdx}.Model.Options.HighPassFilterCutoffHz = 1 / 128; + +for i = 1:numel(conditions) + bm.Nodes{NodeIdx}.Model.X = design_matrix; +end + +% contrast for each condition against baseline +bm.Nodes{NodeIdx}.DummyContrasts = struct('type', 't', ... + 'Contrasts', {design_matrix}); + +% contrast of interests +% alcohol > neural +% negative > neutral +% positive > neutral +conditions_to_check = {'alcohol', 'negative', 'positive'}; +for i = 1:numel(conditions_to_check) + contrast = struct('Test', 't', ... + 'Name', [conditions_to_check{i} '_gt_neutral'], ... + 'Weights', [1, -1], ... + 'ConditionList', {{['trial_type.' conditions_to_check{i}], 'trial_type.neutral'}}); + bm.Nodes{NodeIdx}.Contrasts{i} = contrast; +end + +bm.write('smdl.json'); + +%% create events TSV file +IBI = 25; +ISI = 0; +stimDuration = 4; +stimPerBlock = 5; +nbBlocks = 8; + +TR = 1; + +randomize = true; + +trial_type = {}; +onset = []; +duration = []; + +time = 0; + +for iBlock = 1:nbBlocks + + condition_order = 1:4; + if randomize + condition_order = randperm(4); + end + + for i = 1:numel(conditions) + + condition_idx = condition_order(i); + + for iTrial = 1:stimPerBlock + trial_type{end + 1} = conditions{condition_idx}; %#ok<*SAGROW> + onset(end + 1) = time; + duration(end + 1) = stimDuration; + time = time + stimDuration + ISI; + end + + time = time + IBI; + + end +end + +tsv = struct('trial_type', {trial_type}, 'onset', onset, 'duration', duration'); + +bids.util.tsvwrite('events.tsv', tsv); + +opt.TR = TR; + +opt.model.file = fullfile(pwd, 'smdl.json'); + +[e, X] = computeDesignEfficiency(fullfile(pwd, 'events.tsv'), opt); + +%% +% saveas(gcf(), 'design_matrix.png'); +% close(gcf()); +% +% saveas(gcf(), 'events.png'); +% close(gcf()); diff --git a/demos/efficiency/events.tsv b/demos/efficiency/events.tsv new file mode 100644 index 000000000..4098319b3 --- /dev/null +++ b/demos/efficiency/events.tsv @@ -0,0 +1,161 @@ +trial_type onset duration +neutral 0 4 +neutral 4 4 +neutral 8 4 +neutral 12 4 +neutral 16 4 +positive 45 4 +positive 49 4 +positive 53 4 +positive 57 4 +positive 61 4 +negative 90 4 +negative 94 4 +negative 98 4 +negative 102 4 +negative 106 4 +alcohol 135 4 +alcohol 139 4 +alcohol 143 4 +alcohol 147 4 +alcohol 151 4 +positive 180 4 +positive 184 4 +positive 188 4 +positive 192 4 +positive 196 4 +negative 225 4 +negative 229 4 +negative 233 4 +negative 237 4 +negative 241 4 +neutral 270 4 +neutral 274 4 +neutral 278 4 +neutral 282 4 +neutral 286 4 +alcohol 315 4 +alcohol 319 4 +alcohol 323 4 +alcohol 327 4 +alcohol 331 4 +neutral 360 4 +neutral 364 4 +neutral 368 4 +neutral 372 4 +neutral 376 4 +negative 405 4 +negative 409 4 +negative 413 4 +negative 417 4 +negative 421 4 +positive 450 4 +positive 454 4 +positive 458 4 +positive 462 4 +positive 466 4 +alcohol 495 4 +alcohol 499 4 +alcohol 503 4 +alcohol 507 4 +alcohol 511 4 +negative 540 4 +negative 544 4 +negative 548 4 +negative 552 4 +negative 556 4 +positive 585 4 +positive 589 4 +positive 593 4 +positive 597 4 +positive 601 4 +neutral 630 4 +neutral 634 4 +neutral 638 4 +neutral 642 4 +neutral 646 4 +alcohol 675 4 +alcohol 679 4 +alcohol 683 4 +alcohol 687 4 +alcohol 691 4 +negative 720 4 +negative 724 4 +negative 728 4 +negative 732 4 +negative 736 4 +positive 765 4 +positive 769 4 +positive 773 4 +positive 777 4 +positive 781 4 +neutral 810 4 +neutral 814 4 +neutral 818 4 +neutral 822 4 +neutral 826 4 +alcohol 855 4 +alcohol 859 4 +alcohol 863 4 +alcohol 867 4 +alcohol 871 4 +positive 900 4 +positive 904 4 +positive 908 4 +positive 912 4 +positive 916 4 +alcohol 945 4 +alcohol 949 4 +alcohol 953 4 +alcohol 957 4 +alcohol 961 4 +neutral 990 4 +neutral 994 4 +neutral 998 4 +neutral 1002 4 +neutral 1006 4 +negative 1035 4 +negative 1039 4 +negative 1043 4 +negative 1047 4 +negative 1051 4 +alcohol 1080 4 +alcohol 1084 4 +alcohol 1088 4 +alcohol 1092 4 +alcohol 1096 4 +negative 1125 4 +negative 1129 4 +negative 1133 4 +negative 1137 4 +negative 1141 4 +neutral 1170 4 +neutral 1174 4 +neutral 1178 4 +neutral 1182 4 +neutral 1186 4 +positive 1215 4 +positive 1219 4 +positive 1223 4 +positive 1227 4 +positive 1231 4 +negative 1260 4 +negative 1264 4 +negative 1268 4 +negative 1272 4 +negative 1276 4 +neutral 1305 4 +neutral 1309 4 +neutral 1313 4 +neutral 1317 4 +neutral 1321 4 +alcohol 1350 4 +alcohol 1354 4 +alcohol 1358 4 +alcohol 1362 4 +alcohol 1366 4 +positive 1395 4 +positive 1399 4 +positive 1403 4 +positive 1407 4 +positive 1411 4 diff --git a/demos/efficiency/smdl.json b/demos/efficiency/smdl.json new file mode 100644 index 000000000..5b4346423 --- /dev/null +++ b/demos/efficiency/smdl.json @@ -0,0 +1,105 @@ +{ + "Name": "empty_model", + "BIDSModelVersion": "1.0.0", + "Description": "This is an empty BIDS stats model.", + "Input": { + "task": [ + "" + ] + }, + "Nodes": [ + { + "Level": "Run", + "Name": "run", + "GroupBy": [ + "" + ], + "Transformations": { + "Transformer": "", + "Instructions": [ + { + "Name": "", + "Inputs": [ + "" + ] + } + ] + }, + "Model": { + "X": [ + "trial_type.alcohol", + "trial_type.neutral", + "trial_type.negative", + "trial_type.positive" + ], + "Type": "glm", + "HRF": { + "Variables": [ + "" + ], + "Model": "DoubleGamma" + }, + "Options": { + "HighPassFilterCutoffHz": 0.0078125, + "Mask": { + "desc": [ + "brain" + ], + "suffix": [ + "mask" + ] + } + }, + "Software": [] + }, + "Contrasts": [ + { + "Test": "t", + "Name": "alcohol_gt_neutral", + "Weights": [ + 1, + -1 + ], + "ConditionList": [ + "trial_type.alcohol", + "trial_type.neutral" + ] + }, + { + "Test": "t", + "Name": "negative_gt_neutral", + "Weights": [ + 1, + -1 + ], + "ConditionList": [ + "trial_type.negative", + "trial_type.neutral" + ] + }, + { + "Test": "t", + "Name": "positive_gt_neutral", + "Weights": [ + 1, + -1 + ], + "ConditionList": [ + "trial_type.positive", + "trial_type.neutral" + ] + } + ], + "DummyContrasts": { + "type": "t", + "Contrasts": [ + "trial_type.alcohol", + "trial_type.neutral", + "trial_type.negative", + "trial_type.positive" + ] + } + } + ], + "Edges": [] +} diff --git a/src/QA/computeDesignEfficiency.m b/src/QA/computeDesignEfficiency.m index 034120cf4..c58cc5593 100644 --- a/src/QA/computeDesignEfficiency.m +++ b/src/QA/computeDesignEfficiency.m @@ -1,4 +1,4 @@ -function e = computeDesignEfficiency(tsvFile, opt) +function [e, X] = computeDesignEfficiency(tsvFile, opt) % % Calculate efficiency for fMRI GLMs. Relies on Rik Henson's fMRI_GLM_efficiency function. % @@ -46,59 +46,59 @@ % % EXAMPLE: % - % .. code-block:: guess + % .. code-block:: matlab % - % %% create stats model JSON - % json = createEmptyStatsModel(); - % runStepIdx = 1; - % json.Steps{runStepIdx}.Model.X = {'trial_type.cdt_A', 'trial_type.cdt_B'}; - % json.Steps{runStepIdx}.DummyContrasts = {'trial_type.cdt_A', 'trial_type.cdt_B'}; + % %% create stats model JSON + % bm = BidsModel(); + % NodeIdx = 1; + % bm.Nodes{NodeIdx}.Model.X = {'trial_type.cdt_A', 'trial_type.cdt_B'}; + % bm.Nodes{NodeIdx}.DummyContrasts = struct('type', 't', ... + % 'Contrasts', {{'trial_type.cdt_A', 'trial_type.cdt_B'}}); % - % contrast = struct('type', 't', ... - % 'Name', 'A_gt_B', ... - % 'weights', [1, -1], ... - % 'ConditionList', {{'trial_type.cdt_A', 'trial_type.cdt_B'}}); + % contrast = struct('type', 't', ... + % 'Name', 'A_gt_B', ... + % 'Weights', [1, -1], ... + % 'ConditionList', {{'trial_type.cdt_A', 'trial_type.cdt_B'}}); % - % json.Steps{runStepIdx}.Contrasts = contrast; + % bm.Nodes{NodeIdx}.Contrasts{1} = contrast; % - % bids.util.jsonwrite('smdl.json', json); + % bm.write('smdl.json'); % - % %% create events TSV file - % conditions = {'cdt_A', 'cdt_B'}; - % IBI = 5; - % ISI = 0.1; - % stimDuration = 1.5; - % stimPerBlock = 12; - % nbBlocks = 10; + % %% create events TSV file + % conditions = {'cdt_A', 'cdt_B'}; + % IBI = 5; + % ISI = 0.1; + % stimDuration = 1.5; + % stimPerBlock = 12; + % nbBlocks = 10; % - % trial_type = {}; - % onset = []; - % duration = []; + % trial_type = {}; + % onset = []; + % duration = []; % - % time = 0; + % time = 0; % - % for iBlock = 1:nbBlocks - % for cdt = 1:numel(conditions) - % for iTrial = 1:stimPerBlock - % trial_type{end + 1} = conditions{cdt}; - % onset(end + 1) = time; - % duration(end + 1) = stimDuration; - % time = time + stimDuration + ISI; - % end - % time = time + IBI; - % end - % end + % for iBlock = 1:nbBlocks + % for cdt = 1:numel(conditions) + % for iTrial = 1:stimPerBlock + % trial_type{end + 1} = conditions{cdt}; %#ok<*SAGROW> + % onset(end + 1) = time; + % duration(end + 1) = stimDuration; + % time = time + stimDuration + ISI; + % end + % time = time + IBI; + % end + % end % - % tsv = struct('trial_type', {trial_type}, 'onset', onset, 'duration', duration'); + % tsv = struct('trial_type', {trial_type}, 'onset', onset, 'duration', duration'); % - % bids.util.tsvwrite('events.tsv', tsv); + % bids.util.tsvwrite('events.tsv', tsv); % - % opt.TR = 2; + % opt.TR = 2; % - % opt.model.file = fullfile(pwd, 'smdl.json'); - % - % e = computeDesignEfficiency(fullfile(pwd, 'events.tsv'), opt); + % opt.model.file = fullfile(pwd, 'smdl.json'); % + % e = computeDesignEfficiency(fullfile(pwd, 'events.tsv'), opt); % % (C) Copyright 2021 Remi Gau @@ -147,7 +147,9 @@ end %% - [e, ~, ~, X] = fMRI_GLM_efficiency(opt); + [e, ~, ~, X_detrended] = fMRI_GLM_efficiency(opt); + + X = createDesignMatrix(opt); if opt.verbosity for i = 1:numel(opt.CM) @@ -157,9 +159,10 @@ if opt.verbosity && ~bids.internal.is_github_ci() - plotEvents(tsvFile); + bids.util.plot_events(tsvFile); - figure('name', 'design matrix', 'position', [50 50 1000 1000]); + figure('name', 'detrended filtered design matrix', ... + 'position', [50 50 1000 1000]); % TODO add a second axis with scale in seconds @@ -227,7 +230,9 @@ function plotFft(signal, rt, HPF) plot(Hz(q), gX(q, :)); - patch([0 1 1 0] / HPF, [0 0 1 1] * max(max(gX)), [1 1 1] * .9); + p1 = patch([0 1 1 0] / HPF, [0 0 1 1] * max(max(gX)), [1 1 1] * 0.5); + p1.FaceVertexAlphaData = 0.5; + p1.FaceAlpha = 'flat'; xlabel('Frequency (Hz)'); ylabel('relative spectral density'); diff --git a/src/QA/createDesignMatrix.m b/src/QA/createDesignMatrix.m new file mode 100644 index 000000000..571955646 --- /dev/null +++ b/src/QA/createDesignMatrix.m @@ -0,0 +1,81 @@ +function X = createDesignMatrix(S) + % + % Return the design matrix for the structure S. + % + % Adapted from fMRI_GLM_efficiency + % that only returns a filtered and detrended design matrix. + % + % USAGE:: + % + % X = createDesignMatrix(S) + % + % :param S: input structure with the following fields + % :type S: structure + % + % - ``S.Ns`` : number of scans + % - ``S.TR`` : inter-scan interval (s) + % - ``S.sots`` : cell array of onset times for each condition in units of scans + % - ``S.durs`` : cell array of durations for each condition in units of scans + % - ``S.bf`` : type of HRF basis function (based on spm_get_bf.m), + % default = "hrf" + + % (C) Copyright 2021 Remi Gau + + Ns = S.Ns; + + TR = S.TR; + + % Minimum time for convolution space (seconds) + dt = 0.2; + st = round(TR / dt); + + sots = S.sots; + + durs = S.durs; + for j = 1:length(sots) + if length(durs{j}) == 1 + durs{j} = repmat(durs{j}, 1, length(sots{j})); + elseif length(durs{j}) ~= length(sots{j}) + error('Number of durations (%d) does not match number of onsets (%d) for condition %d', ... + length(durs{j}), length(sots{j}), j); + end + end + + try + bf = S.bf; + catch + bf = 'hrf'; % SPM's canonical HRF + % bf = 'Finite Impulse Response' % FIR + end + + if ischar(bf) + xBF.dt = dt; + xBF.name = bf; + xBF.length = 30; + xBF.order = 30 / TR; + bf = spm_get_bf(xBF); + bf = bf.bf; + else + % Assume bf is a TxNk matrix of basis functions in units of dt + end + + Nj = length(sots); + Nk = size(bf, 2); + Nt = Ns * st; + + %% Create neural timeseries and concolve + s = 1:st:Nt; + X = zeros(Ns, Nj * Nk); + for j = 1:Nj + + u = zeros(Nt, 1); + for i = 1:length(sots{j}) + o = [sots{j}(i) * st:st:(sots{j}(i) + durs{j}(i)) * st]; + u(round(o) + 1) = 1; + end + + for k = 1:Nk + b = conv(u, bf(:, k)); + X(:, (j - 1) * Nk + k) = b(s); + end + end