From 9911f08d6dfccfef8400c8d161f11814207cb6a6 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Wed, 17 Jun 2026 19:23:49 -0400 Subject: [PATCH 1/8] Add glm_map object class scaffold New scikit-learn-style estimator class that bundles GLM design specification, fitted result maps (statistic_image), and design diagnostics in one container. Composition over fmri_glm_design_matrix (wrapped in .design) and fmri_data.regress (the compute engine). - classdef with stored properties (design/level/is_timeseries, betas/t/contrast maps, vif/leverage/collinearity diagnostics, provenance) and true Dependent accessors (X, TR, onsets, durations, regressor_names, num_*, is_fitted) that read through to .design. - Implemented + MATLAB-verified: diagnostics (VIF/cVIF/leverage/ condition number/rank/redundant-column report), add_contrasts, threshold, montage, table, plot_design, summary, check_properties, private select_map. - Documented stubs with planned field mappings: fit (wraps regress), build_design (wraps fmri_glm_design_matrix.build), import_SPM (SPM12/SPM25 -> .design). Co-Authored-By: Claude Opus 4.8 (1M context) --- CanlabCore/@glm_map/add_contrasts.m | 77 ++++ CanlabCore/@glm_map/build_design.m | 49 +++ CanlabCore/@glm_map/check_properties.m | 73 ++++ CanlabCore/@glm_map/diagnostics.m | 179 +++++++++ CanlabCore/@glm_map/fit.m | 73 ++++ CanlabCore/@glm_map/glm_map.m | 446 +++++++++++++++++++++++ CanlabCore/@glm_map/import_SPM.m | 54 +++ CanlabCore/@glm_map/montage.m | 57 +++ CanlabCore/@glm_map/plot_design.m | 74 ++++ CanlabCore/@glm_map/private/select_map.m | 35 ++ CanlabCore/@glm_map/summary.m | 43 +++ CanlabCore/@glm_map/table.m | 50 +++ CanlabCore/@glm_map/threshold.m | 71 ++++ 13 files changed, 1281 insertions(+) create mode 100644 CanlabCore/@glm_map/add_contrasts.m create mode 100644 CanlabCore/@glm_map/build_design.m create mode 100644 CanlabCore/@glm_map/check_properties.m create mode 100644 CanlabCore/@glm_map/diagnostics.m create mode 100644 CanlabCore/@glm_map/fit.m create mode 100644 CanlabCore/@glm_map/glm_map.m create mode 100644 CanlabCore/@glm_map/import_SPM.m create mode 100644 CanlabCore/@glm_map/montage.m create mode 100644 CanlabCore/@glm_map/plot_design.m create mode 100644 CanlabCore/@glm_map/private/select_map.m create mode 100644 CanlabCore/@glm_map/summary.m create mode 100644 CanlabCore/@glm_map/table.m create mode 100644 CanlabCore/@glm_map/threshold.m diff --git a/CanlabCore/@glm_map/add_contrasts.m b/CanlabCore/@glm_map/add_contrasts.m new file mode 100644 index 00000000..34c47f0c --- /dev/null +++ b/CanlabCore/@glm_map/add_contrasts.m @@ -0,0 +1,77 @@ +function obj = add_contrasts(obj, C, names, varargin) +% Add one or more linear contrasts to a glm_map object. +% +% Appends rows of contrast weights (over regressors) to obj.contrasts and +% names to obj.contrast_names, validating that the weight vectors match the +% number of regressors in the design. +% +% :Usage: +% :: +% +% obj = add_contrasts(obj, C, names) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design available (obj.X non-empty). +% +% **C:** +% A [num_contrasts x num_regressors] matrix; each ROW is one contrast +% over the regressors. (Stored internally as [regressors x contrasts].) +% +% **names:** +% Cell array of contrast names, one per row of C. Optional; defaults +% to 'Con1', 'Con2', ... +% +% :Outputs: +% +% **obj:** +% glm_map with contrasts and contrast_names appended. +% +% :Examples: +% :: +% +% g = glm_map('X', [ones(30,1) zscore((1:30)') randn(30,1)], 'level', 2); +% g = add_contrasts(g, [0 1 0; 0 0 1], {'slope' 'nuisance'}); +% +% :See also: +% - diagnostics, fmri_data.regress +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. +% .. + +if nargin < 3 || isempty(names) + names = {}; +end +if ~iscell(names), names = cellstr(names); end + +nreg = obj.num_regressors; +if nreg == 0 + error('glm_map:NoDesign', 'No design available; set obj.X or build a design before adding contrasts.'); +end + +% Each row of C is a contrast; validate width against number of regressors +if size(C, 2) ~= nreg + error('glm_map:ContrastSize', ... + 'Each contrast must have %d weights (one per regressor); got %d columns in C.', nreg, size(C, 2)); +end + +ncon_new = size(C, 1); + +% Default names +if isempty(names) + start = numel(obj.contrast_names); + names = arrayfun(@(k) sprintf('Con%d', start + k), 1:ncon_new, 'UniformOutput', false); +elseif numel(names) ~= ncon_new + error('glm_map:ContrastNames', 'Number of names (%d) must match number of contrasts (%d).', numel(names), ncon_new); +end + +% Append. Internally contrasts are stored [regressors x contrasts]. +obj.contrasts = [obj.contrasts, C']; +obj.contrast_names = [obj.contrast_names(:); names(:)]'; + +obj.history{end + 1} = sprintf('add_contrasts: added %d contrast(s)', ncon_new); + +end % add_contrasts diff --git a/CanlabCore/@glm_map/build_design.m b/CanlabCore/@glm_map/build_design.m new file mode 100644 index 00000000..716086c0 --- /dev/null +++ b/CanlabCore/@glm_map/build_design.m @@ -0,0 +1,49 @@ +function obj = build_design(obj, varargin) +% Build the design matrix X for an event/1st-level glm_map from onsets. +% +% Delegates to the wrapped fmri_glm_design_matrix object (and ultimately +% onsets2fmridesign) to convolve onsets/durations with the basis set and +% assemble X, regressor names, and any nuisance covariates. After this call, +% the Dependent obj.X and obj.regressor_names read through to the built design. +% +% :Usage: +% :: +% +% obj = build_design(obj, varargin) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a non-empty .design (fmri_glm_design_matrix). +% +% :Optional Inputs: +% +% **'doplot' / 'noplot':** +% Plot the resulting design matrix (default false). +% +% :Outputs: +% +% **obj:** +% The glm_map object with the wrapped design built (design.xX.X populated). +% +% :See also: +% - fmri_glm_design_matrix, fmri_glm_design_matrix.build, onsets2fmridesign +% +% .. +% Programmers' notes: +% SCAFFOLD - not yet implemented. Planned behavior: +% 1. Require obj.level == 1 and ~isempty(obj.design). +% 2. obj.design = build(obj.design); % existing method +% 3. Sync obj.contrasts/contrast_names defaults if requested. +% 4. Append to obj.history. +% .. + +if isempty(obj.design) + error('glm_map:NoDesign', ... + 'build_design requires an fmri_glm_design_matrix in obj.design (event/1st-level mode).'); +end + +error('glm_map:NotImplemented', ... + 'build_design() is scaffolded but not yet implemented. Planned for Phase 3 (wraps fmri_glm_design_matrix.build).'); + +end % build_design diff --git a/CanlabCore/@glm_map/check_properties.m b/CanlabCore/@glm_map/check_properties.m new file mode 100644 index 00000000..42dc60ad --- /dev/null +++ b/CanlabCore/@glm_map/check_properties.m @@ -0,0 +1,73 @@ +function obj = check_properties(obj, varargin) +% Validate and enforce types on a glm_map object's properties. +% +% Lightweight consistency checks: ensures cell-array metadata fields are +% cells, the level is 1 or 2, contrast bookkeeping is consistent, and any +% contained statistic_image maps pass their own check_properties. +% +% :Usage: +% :: +% +% obj = check_properties(obj) +% +% :Inputs: +% +% **obj:** +% A glm_map object. +% +% :Outputs: +% +% **obj:** +% The glm_map object with types/fields coerced where needed. Warns on +% inconsistencies it cannot silently fix. +% +% :See also: +% - statistic_image.check_properties +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. +% .. + +% Cell-array fields +if ~iscell(obj.contrast_names), obj.contrast_names = cellstr(obj.contrast_names); end +if ~iscell(obj.warnings), obj.warnings = {}; end +if ~iscell(obj.history), obj.history = {}; end +if ~iscell(obj.regressor_names_direct) && ~isempty(obj.regressor_names_direct) + obj.regressor_names_direct = cellstr(obj.regressor_names_direct); +end + +% Level +if ~ismember(obj.level, [1 2]) + warning('glm_map:BadLevel', 'level should be 1 or 2; got %s.', num2str(obj.level)); +end + +% is_timeseries should be logical +obj.is_timeseries = logical(obj.is_timeseries); + +% AR only meaningful for timeseries +if obj.level == 2 && obj.is_timeseries + warning('glm_map:LevelTimeseries', 'is_timeseries is true but level is 2 (group); AR models are not appropriate for group images.'); +end + +% Contrast bookkeeping +if ~isempty(obj.contrasts) + if size(obj.contrasts, 1) ~= obj.num_regressors && obj.num_regressors > 0 + warning('glm_map:ContrastSize', 'contrasts has %d rows but design has %d regressors.', ... + size(obj.contrasts, 1), obj.num_regressors); + end + if ~isempty(obj.contrast_names) && numel(obj.contrast_names) ~= size(obj.contrasts, 2) + warning('glm_map:ContrastNames', 'Number of contrast_names (%d) does not match number of contrasts (%d).', ... + numel(obj.contrast_names), size(obj.contrasts, 2)); + end +end + +% Delegate to contained statistic_image maps where possible +for f = {'betas', 't', 'contrast_estimates', 'contrast_t'} + m = obj.(f{1}); + if ~isempty(m) && isa(m, 'statistic_image') && ismethod(m, 'check_properties') + obj.(f{1}) = check_properties(m); + end +end + +end % check_properties diff --git a/CanlabCore/@glm_map/diagnostics.m b/CanlabCore/@glm_map/diagnostics.m new file mode 100644 index 00000000..a9d40225 --- /dev/null +++ b/CanlabCore/@glm_map/diagnostics.m @@ -0,0 +1,179 @@ +function obj = diagnostics(obj, varargin) +% Compute and report design diagnostics for a glm_map object. +% +% Evaluates the conditioning of the design matrix X (and contrasts C): +% variance inflation factors (VIF) per regressor, contrast VIFs (cVIF), +% per-observation leverage, condition number, rank deficiency, and a +% redundant/near-collinear row report. Results are stored back into the +% object and (optionally) printed. +% +% :Usage: +% :: +% +% obj = diagnostics(obj, varargin) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design matrix available (obj.X non-empty). +% +% :Optional Inputs: +% +% **'doverbose' / 'noverbose':** +% Print a summary table (default true). +% +% **'vif_thresh', [value]:** +% VIF warning threshold (default 4). +% +% :Outputs: +% +% **obj:** +% glm_map with vif, contrast_vif, leverages, condition_number, +% rank_deficient, collinearity_report, and warnings populated. +% +% :Examples: +% :: +% +% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2); +% g = diagnostics(g); +% +% :See also: +% - VIF, cVIF, fmri_data.regress +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. Operates on the design only (no fit +% required), so it can be used to screen a design before fitting. +% .. + +% ------------------------------------------------------------------------- +% Parse inputs +% ------------------------------------------------------------------------- +doverbose = true; +if any(strcmpi(varargin, 'noverbose')), doverbose = false; end + +vif_thresh = 4; +wh = find(strcmpi(varargin, 'vif_thresh')); +if ~isempty(wh), vif_thresh = varargin{wh(1) + 1}; end + +X = obj.X; +if isempty(X) + error('glm_map:NoDesign', 'No design matrix available (obj.X is empty). Build or supply a design first.'); +end + +mywarnings = {}; + +% ------------------------------------------------------------------------- +% Variance inflation factors (per regressor) +% ------------------------------------------------------------------------- +obj.vif = VIF(X); + +if any(obj.vif > vif_thresh) + mywarnings{end + 1} = sprintf(['Design multicollinearity: %d regressor(s) have VIF > %g. ' ... + 'Check obj.vif and obj.regressor_names.'], sum(obj.vif > vif_thresh), vif_thresh); +end + +% ------------------------------------------------------------------------- +% Contrast variance inflation factors (per contrast), if contrasts defined +% ------------------------------------------------------------------------- +if ~isempty(obj.contrasts) + if size(obj.contrasts, 1) ~= size(X, 2) + mywarnings{end + 1} = sprintf(['Contrast matrix has %d rows but design has %d regressors; ' ... + 'skipping contrast VIFs.'], size(obj.contrasts, 1), size(X, 2)); + else + % cVIF expects one contrast per row -> transpose [P x K] to [K x P] + obj.contrast_vif = cVIF(X, obj.contrasts'); + end +end + +% ------------------------------------------------------------------------- +% Leverage (per observation) +% ------------------------------------------------------------------------- +H = X * pinv(X); +obj.leverages = diag(H)'; + +if any(abs(zscore(obj.leverages)) >= 3) + mywarnings{end + 1} = ['Some observations have high leverage (abs(z(leverage)) >= 3); ' ... + 'the fit may be unstable. Check obj.leverages.']; +end + +% ------------------------------------------------------------------------- +% Conditioning / rank +% ------------------------------------------------------------------------- +obj.condition_number = cond(X); +obj.rank_deficient = rank(X) < size(X, 2); + +if obj.rank_deficient + mywarnings{end + 1} = 'Design matrix X is rank deficient (rank(X) < number of regressors).'; +end + +% ------------------------------------------------------------------------- +% Redundant / near-collinear column report +% ------------------------------------------------------------------------- +report = struct(); +report.vif_threshold = vif_thresh; +report.high_vif_columns = find(obj.vif > vif_thresh); + +% Duplicate (identical) columns +ncol = size(X, 2); +dup_pairs = []; +for a = 1:ncol - 1 + for b = a + 1:ncol + if isequal(X(:, a), X(:, b)) + dup_pairs(end + 1, :) = [a b]; %#ok + end + end +end +report.duplicate_column_pairs = dup_pairs; +if ~isempty(dup_pairs) + mywarnings{end + 1} = sprintf('%d pair(s) of identical design columns detected (see obj.collinearity_report).', size(dup_pairs, 1)); +end + +% Near-collinear pairs by pairwise correlation magnitude +R = corrcoef(X); +R(logical(eye(ncol))) = 0; +[ia, ib] = find(triu(abs(R) > 0.95, 1)); +report.high_correlation_pairs = [ia ib]; + +obj.collinearity_report = report; + +% ------------------------------------------------------------------------- +% Store warnings +% ------------------------------------------------------------------------- +obj.warnings = [obj.warnings(:); mywarnings(:)]'; +obj.history{end + 1} = 'diagnostics: computed VIF, cVIF, leverage, condition number, collinearity report'; + +% ------------------------------------------------------------------------- +% Report +% ------------------------------------------------------------------------- +if doverbose + rn = obj.regressor_names; + fprintf('\n glm_map diagnostics\n %s\n', repmat('-', 1, 50)); + fprintf(' %-28s %s\n', 'Regressor', 'VIF'); + for i = 1:numel(obj.vif) + if i <= numel(rn) && ~isempty(rn{i}), name = rn{i}; else, name = sprintf('R%d', i); end + flag = ''; if obj.vif(i) > vif_thresh, flag = ' <-- high'; end + fprintf(' %-28s %6.2f%s\n', name, obj.vif(i), flag); + end + if ~isempty(obj.contrast_vif) + fprintf(' %s\n Contrast VIFs:\n', repmat('-', 1, 50)); + for i = 1:numel(obj.contrast_vif) + if i <= numel(obj.contrast_names) && ~isempty(obj.contrast_names{i}) + name = obj.contrast_names{i}; + else + name = sprintf('Con%d', i); + end + fprintf(' %-28s %6.2f\n', name, obj.contrast_vif(i)); + end + end + fprintf(' %s\n', repmat('-', 1, 50)); + fprintf(' condition number : %.2f\n', obj.condition_number); + fprintf(' rank deficient : %d\n', obj.rank_deficient); + if ~isempty(mywarnings) + fprintf(' %d warning(s):\n', numel(mywarnings)); + for i = 1:numel(mywarnings), fprintf(' - %s\n', mywarnings{i}); end + end + fprintf('\n'); +end + +end % diagnostics diff --git a/CanlabCore/@glm_map/fit.m b/CanlabCore/@glm_map/fit.m new file mode 100644 index 00000000..79f20af8 --- /dev/null +++ b/CanlabCore/@glm_map/fit.m @@ -0,0 +1,73 @@ +function obj = fit(obj, data, varargin) +% Fit a mass-univariate GLM and populate result maps in a glm_map object. +% +% Scikit-learn-style fit: builds the design (if event/1st-level), runs +% fmri_data.regress on the data, and unpacks the outputs into the object's +% statistic_image result maps (betas, t, contrast_estimates, contrast_t), +% along with sigma, dfe, and the design diagnostics. +% +% :Usage: +% :: +% +% obj = fit(obj, data, varargin) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design specified (either a .design event +% model or a direct .X matrix). +% +% **data:** +% An fmri_data object whose images correspond to the rows of X +% (timeseries for level 1, contrast/subject images for level 2). +% +% :Optional Inputs: +% +% **'robust':** +% Use robust (bisquare) regression. Passed through to fmri_data.regress. +% +% **'AR', [order]:** +% Use an autoregressive error model of the given order. Only valid +% when obj.is_timeseries is true. +% +% **'residuals':** +% Also store residuals in obj.residuals. +% +% **'doverbose' / 'noverbose':** +% Toggle verbose output (default true). +% +% :Outputs: +% +% **obj:** +% The glm_map object with result maps and diagnostics populated, and +% is_fitted == true. +% +% :Examples: +% :: +% +% dat = load_image_set('emotionreg'); +% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2); +% g = fit(g, dat); +% +% :See also: +% - fmri_data.regress, build_design, diagnostics +% +% .. +% Programmers' notes: +% SCAFFOLD - not yet implemented. Planned behavior: +% 1. If event/1st-level and X not built, call build_design(obj). +% 2. Attach obj.X to data.X (and obj.contrasts/contrast_names). +% 3. If obj.is_timeseries and 'AR' requested, pass 'AR' to regress; +% otherwise error if 'AR' requested on non-timeseries data. +% 4. out = regress(data, ...): unpack out.b -> obj.betas, out.t -> obj.t, +% out.contrast_images -> obj.contrast_estimates, out.con_t -> +% obj.contrast_t, out.sigma -> obj.sigma, out.df -> obj.dfe. +% 5. Populate obj.vif/contrast_vif/leverages/warnings from out.diagnostics +% (or call diagnostics(obj)). +% 6. Record options in obj.fit_parameters and append to obj.history. +% .. + +error('glm_map:NotImplemented', ... + 'fit() is scaffolded but not yet implemented. Planned for Phase 2 (wraps fmri_data.regress).'); + +end % fit diff --git a/CanlabCore/@glm_map/glm_map.m b/CanlabCore/@glm_map/glm_map.m new file mode 100644 index 00000000..e63610af --- /dev/null +++ b/CanlabCore/@glm_map/glm_map.m @@ -0,0 +1,446 @@ +% glm_map Object class for mass-univariate GLM / multiple-regression fits. +% +% A scikit-learn-style estimator object that bundles, in a single container: +% (1) the *design specification* (onsets, durations, event names/types, +% parametric modulators, basis set, covariates, design matrix X, and +% contrasts) -- wrapping an fmri_glm_design_matrix object, +% (2) the *fitted result maps* (betas, t, contrast estimates, contrast t), +% stored as statistic_image objects, plus sigma and residuals, and +% (3) the *design diagnostics* (variance inflation factors, contrast VIFs, +% leverage, condition number, rank, collinearity/redundancy checks). +% +% glm_map holds the outputs of a mass-univariate GLM fit produced by +% fmri_data.regress. The workflow mirrors scikit-learn estimators: +% +% g = glm_map(...) % construct: set design + fit options (no data) +% g = fit(g, fmri_data_obj) % build design (if needed) and run regression +% table(g); montage(g); % inspect results +% +% - This is a standalone container class (composition, not an image_vector +% subclass): the voxel data live inside the contained statistic_image / +% fmri_data objects (betas, t, contrast_estimates, ...). +% - It supports two design modes: +% * level 1 (event / 1st-level): an fmri_glm_design_matrix in .design +% builds X by HRF convolution of onsets. Mark .is_timeseries = true +% for within-run BOLD data so autoregressive ('AR') models are valid. +% * level 2 (direct / group): a design matrix X is supplied directly. +% - Several convenience attributes (X, regressor_names, onsets, durations, +% condition_names, TR, num_*) are implemented as true MATLAB *Dependent* +% properties that read through to the wrapped design object, so there is a +% single source of truth and no risk of stale duplicates. +% +% :Usage: +% :: +% +% g = glm_map(varargin) +% g = glm_map(fmri_glm_design_matrix_obj) % wrap a 1st-level design +% g = glm_map('X', X, 'level', 2) % direct/group design +% g = glm_map('fieldname', value, ...) % set any stored property +% +% :Inputs: +% +% **(optional first arg) fmri_glm_design_matrix object:** +% If the first argument is an fmri_glm_design_matrix, it is stored in +% .design and the object is marked as a 1st-level (event) model. +% +% **'fieldname', value pairs:** +% Any settable property of glm_map followed by a value. See +% properties(glm_map) and the property_descriptions field. +% +% :Optional Inputs: +% +% **'X', [obs x regressors] matrix:** +% A pre-built design matrix for direct/group (level-2) analysis. +% +% **'level', 1 | 2:** +% Analysis level. 1 = first-level (within-run); 2 = second-level +% (group). Default 2 for direct designs, 1 when a design object is given. +% +% **'is_timeseries', [logical]:** +% True if the data to be fit are a within-run BOLD timeseries; enables +% autoregressive (AR) error models in fit. Default false. +% +% **'contrasts', [regressors x contrasts] matrix:** +% Contrast matrix C (rows must match number of regressors in X). +% +% **'contrast_names', {cellstr}:** +% Names for each contrast (columns of C). +% +% :Outputs: +% +% **g:** +% A glm_map object. +% +% :Examples: +% :: +% +% % ---- Direct / group (2nd-level) design ---- +% dat = load_image_set('emotionreg'); % 30 contrast images +% X = [ones(30,1) zscore((1:30)')]; % intercept + a covariate +% g = glm_map('X', X, 'level', 2, ... +% 'regressor_names', {'intercept' 'cov'}); +% g = fit(g, dat); % runs fmri_data.regress +% diagnostics(g); % VIFs, leverage, etc. +% table(g); montage(g, 't'); +% +% % ---- Event / 1st-level design wrapping fmri_glm_design_matrix ---- +% d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... +% 'onsets', ons, 'condition_names', names); +% g = glm_map(d); % level 1, event mode +% g.is_timeseries = true; % enable AR models in fit +% g = build_design(g); % onsets -> X via convolution +% g.onsets % read-through to .design +% +% % ---- Import from an SPM (SPM12/SPM25) first-level model ---- +% g = import_SPM(glm_map, '/path/to/SPM.mat'); +% +% :See also: +% - fmri_data, fmri_data.regress +% - statistic_image +% - fmri_glm_design_matrix +% - VIF, cVIF +% +% .. +% Author and copyright information: +% +% Copyright (C) 2026 Tor Wager +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% .. +% +% Programmers' notes: +% 2026 - Initial scaffold. Constructor, Dependent accessors, and disp are +% functional; fit/build_design/diagnostics/import_SPM/display/table are +% stubbed for later implementation phases. +% .. + +classdef glm_map + + % --------------------------------------------------------------------- + % Stored properties + % --------------------------------------------------------------------- + properties + + % --- Design specification (1st-level / event mode) --------------- + design % fmri_glm_design_matrix object (wrapped). Holds onsets, durations, event names/types, parametric modulators, basis set, and built X. Empty for pure direct/group designs. + + % --- Design specification (shared / direct mode) ---------------- + level = 2; % Analysis level: 1 = first-level (within-run), 2 = second-level (group) + is_timeseries = false; % Logical; true if data are a within-run BOLD timeseries (enables AR error models in fit) + contrasts = []; % [regressors x contrasts] contrast matrix C + contrast_names = {}; % Cell array of contrast names, one per column of C + + % --- Fitted result maps (populated by fit) ---------------------- + betas % statistic_image (type 'Beta'), [voxels x regressors] + t % statistic_image (type 'T') for regressors, [voxels x regressors] + contrast_estimates % statistic_image (type 'Contrast'), [voxels x contrasts] + contrast_t % statistic_image (type 'T') for contrasts, [voxels x contrasts] + sigma % fmri_data object: residual standard deviation per voxel + dfe % Error degrees of freedom for the fit + residuals % fmri_data object: residuals [observations x voxels] (optional; only if requested) + + % --- Design diagnostics (populated by fit/diagnostics) ---------- + vif % Row vector of variance inflation factors, one per regressor (from VIF.m) + contrast_vif % Vector of contrast variance inflation factors, one per contrast (from cVIF.m) + leverages % Per-observation leverage values, diag(X*pinv(X)) + condition_number % Condition number of the design matrix X + rank_deficient = false; % Logical; true if rank(X) < number of regressors + collinearity_report % Struct with redundant/duplicate row checks, near-collinear pairs, and centering flags + warnings = {}; % Cell array of warning messages accumulated during build/fit + + % --- Provenance / metadata -------------------------------------- + analysis_name = ''; % Short descriptive name for this analysis + fit_parameters = struct(); % Struct recording options used in fit (robust, AR order, threshold, grandmeanscale, ...) + notes = ''; % Free-text notes + history = {}; % Cell array of one-line provenance strings + + % --- Stored backing for direct-mode design ---------------------- + % (Use the Dependent .X / .regressor_names accessors instead of + % touching these directly.) + Xdirect = []; % [obs x regressors] design matrix for direct/group mode (no event model) + regressor_names_direct = {}; % Cell array of regressor names for direct mode + + property_descriptions = { ... + 'design: fmri_glm_design_matrix object holding the 1st-level event design (onsets, durations, names, basis set, built X)' ... + 'level: 1 = first-level (within-run), 2 = second-level (group)' ... + 'is_timeseries: logical, true if data are a within-run BOLD timeseries (enables AR error models)' ... + 'contrasts: [regressors x contrasts] contrast matrix C' ... + 'contrast_names: cell array of contrast names' ... + 'betas/t: statistic_image maps of regression coefficients and their t-statistics, [voxels x regressors]' ... + 'contrast_estimates/contrast_t: statistic_image maps for linear contrasts, [voxels x contrasts]' ... + 'sigma: fmri_data object with residual standard deviation per voxel' ... + 'dfe: error degrees of freedom' ... + 'residuals: fmri_data object with residuals (optional)' ... + 'vif/contrast_vif: variance inflation factors for regressors and contrasts' ... + 'leverages: per-observation leverage values' ... + 'condition_number/rank_deficient: design matrix conditioning diagnostics' ... + 'collinearity_report: struct of redundant-row and near-collinearity checks' ... + 'warnings: cell array of warnings from build/fit' ... + 'analysis_name/fit_parameters/notes/history: provenance and metadata' ... + }; + + end % stored properties + + + % --------------------------------------------------------------------- + % Dependent properties (true MATLAB Dependent; computed accessors) + % --------------------------------------------------------------------- + properties (Dependent) + + TR % Repetition time (s). Reads/writes design.TR for event models. + X % [obs x regressors] design matrix. Reads built X from .design (event) or .Xdirect (direct). + regressor_names % Cell array of regressor (design column) names + onsets % Cell array of event onsets (read-through to .design; 1st-level only) + durations % Cell array of event durations (read-through to .design; 1st-level only) + condition_names % Cell array of condition/event names (read-through to .design; 1st-level only) + num_images % Number of images/observations (rows of X) + num_regressors % Number of regressors (columns of X) + num_contrasts % Number of contrasts (columns of C) + is_fitted % Logical; true once result maps (.betas) are populated + + end % dependent properties + + + methods + + % ================================================================= + % Constructor + % ================================================================= + function obj = glm_map(varargin) + + % Empty object: return defaults + if nargin == 0 + return + end + + % If first argument is an fmri_glm_design_matrix, wrap it as a + % 1st-level (event) design and consume that argument. + if ~isempty(varargin) && isa(varargin{1}, 'fmri_glm_design_matrix') + obj.design = varargin{1}; + obj.level = 1; + varargin(1) = []; + end + + % Names of stored (settable) properties for generic assignment + stored_names = properties('glm_map'); % stored props only (Dependent excluded by properties()) + + % Names of Dependent properties that have setters + settable_dependent = {'TR', 'X', 'regressor_names'}; + + for i = 1:length(varargin) + + if ~ischar(varargin{i}), continue, end + + fieldname = varargin{i}; + + % Map a couple of friendly aliases to backing storage + switch fieldname + case {'names', 'variable_names'} + fieldname = 'regressor_names'; + end + + if any(strcmp(fieldname, stored_names)) || any(strcmp(fieldname, settable_dependent)) + + obj.(fieldname) = varargin{i + 1}; + + % Blank the consumed value so a trailing char value is + % not re-interpreted as a keyword on the next iteration + if ischar(varargin{i + 1}) + varargin{i + 1} = []; + end + + else + warning('glm_map:UnknownField', 'Unknown glm_map field: %s', fieldname); + end + + end % parse inputs + + end % constructor + + + % ================================================================= + % Dependent property GET accessors + % ================================================================= + function val = get.TR(obj) + if ~isempty(obj.design) + val = obj.design.TR; + else + val = []; + end + end + + function val = get.X(obj) + % Prefer a built design matrix from the wrapped design object; + % fall back to the direct-mode matrix. + val = []; + if ~isempty(obj.design) && isstruct(obj.design.xX) ... + && isfield(obj.design.xX, 'X') && ~isempty(obj.design.xX.X) + val = obj.design.xX.X; + elseif ~isempty(obj.Xdirect) + val = obj.Xdirect; + end + end + + function val = get.regressor_names(obj) + val = {}; + if ~isempty(obj.design) && isstruct(obj.design.xX) ... + && isfield(obj.design.xX, 'name') && ~isempty(obj.design.xX.name) + val = obj.design.xX.name; + elseif ~isempty(obj.regressor_names_direct) + val = obj.regressor_names_direct; + end + end + + function val = get.onsets(obj) + val = local_collect_U_field(obj.design, 'ons'); + end + + function val = get.durations(obj) + val = local_collect_U_field(obj.design, 'dur'); + end + + function val = get.condition_names(obj) + val = local_collect_U_field(obj.design, 'name'); + end + + function val = get.num_images(obj) + val = size(obj.X, 1); + end + + function val = get.num_regressors(obj) + val = size(obj.X, 2); + end + + function val = get.num_contrasts(obj) + val = size(obj.contrasts, 2); + end + + function val = get.is_fitted(obj) + val = ~isempty(obj.betas); + end + + + % ================================================================= + % Dependent property SET accessors + % ================================================================= + function obj = set.TR(obj, val) + if isempty(obj.design) + error('glm_map:NoDesign', ... + 'Cannot set TR: no fmri_glm_design_matrix in .design. Create one first (event/1st-level mode).'); + end + obj.design.TR = val; + end + + function obj = set.X(obj, val) + % Setting X targets the direct/group design backing store. + if ~isempty(obj.design) && isstruct(obj.design.xX) ... + && isfield(obj.design.xX, 'X') && ~isempty(obj.design.xX.X) + warning('glm_map:DesignPresent', ... + ['This glm_map has a built event design in .design; setting .X stores a direct ' ... + 'matrix in .Xdirect that the event design will shadow. Use build_design instead.']); + end + obj.Xdirect = val; + end + + function obj = set.regressor_names(obj, val) + if ~iscell(val), val = cellstr(val); end + obj.regressor_names_direct = val; + end + + + % ================================================================= + % disp: concise object summary + % ================================================================= + function disp(obj) + + fprintf(' glm_map object\n'); + fprintf(' %s\n', repmat('-', 1, 60)); + + if ~isempty(obj.analysis_name) + fprintf(' analysis_name : %s\n', obj.analysis_name); + end + + switch obj.level + case 1, levelstr = '1 (first-level / within-run)'; + case 2, levelstr = '2 (second-level / group)'; + otherwise, levelstr = num2str(obj.level); + end + fprintf(' level : %s\n', levelstr); + fprintf(' is_timeseries : %d\n', obj.is_timeseries); + + if ~isempty(obj.design) + fprintf(' design : fmri_glm_design_matrix (event mode), TR = %g\n', obj.design.TR); + else + fprintf(' design : (direct/group mode, no event design)\n'); + end + + fprintf(' X : %d images x %d regressors\n', obj.num_images, obj.num_regressors); + fprintf(' contrasts : %d\n', obj.num_contrasts); + + if obj.is_fitted + fprintf(' fitted : YES (betas, t%s populated; dfe = %s)\n', ... + local_tf(~isempty(obj.contrast_estimates), ', contrasts'), num2str(obj.dfe)); + else + fprintf(' fitted : no (run fit(obj, fmri_data_obj))\n'); + end + + if ~isempty(obj.vif) + fprintf(' max VIF : %.2f\n', max(obj.vif)); + end + + if ~isempty(obj.warnings) + fprintf(' warnings : %d (see obj.warnings)\n', numel(obj.warnings)); + end + + fprintf(' %s\n', repmat('-', 1, 60)); + fprintf(' methods(glm_map) for a list of operations.\n\n'); + + end % disp + + end % methods + +end % classdef + + +% ===================================================================== +% Local helper functions +% ===================================================================== +function out = local_collect_U_field(design, fieldname) +% Read-through accessor: collect a field (ons/dur/name) from the nested +% design.Sess(s).U(i) structure into a flat cell array. Returns {} if the +% wrapped design has no session/onset structure yet. + +out = {}; +if isempty(design) || ~isprop(design, 'Sess') || isempty(design.Sess) + return +end + +for s = 1:numel(design.Sess) + if ~isfield(design.Sess(s), 'U') || isempty(design.Sess(s).U) + continue + end + for u = 1:numel(design.Sess(s).U) + if isfield(design.Sess(s).U(u), fieldname) + out{end + 1} = design.Sess(s).U(u).(fieldname); %#ok + end + end +end + +end % local_collect_U_field + + +function s = local_tf(tf, str) +% Return str if tf is true, else ''. +if tf, s = str; else, s = ''; end +end % local_tf diff --git a/CanlabCore/@glm_map/import_SPM.m b/CanlabCore/@glm_map/import_SPM.m new file mode 100644 index 00000000..fe609cd9 --- /dev/null +++ b/CanlabCore/@glm_map/import_SPM.m @@ -0,0 +1,54 @@ +function obj = import_SPM(obj, SPM, varargin) +% Import a first-level design (and optional betas) from an SPM model. +% +% Populates a glm_map object from an SPM structure or an SPM.mat path, +% mapping SPM's design fields into the wrapped fmri_glm_design_matrix object. +% Compatible with SPM12 and SPM25 first-level models. +% +% :Usage: +% :: +% +% obj = import_SPM(glm_map, SPM) % SPM struct already loaded +% obj = import_SPM(glm_map, '/path/SPM.mat') +% +% :Inputs: +% +% **obj:** +% A glm_map object (typically empty: glm_map). +% +% **SPM:** +% Either an SPM structure (as loaded from SPM.mat) or a char/string +% path to an SPM.mat file. +% +% :Optional Inputs: +% +% **'load_betas':** +% Also load the estimated beta_*.nii images into obj.betas. +% +% :Outputs: +% +% **obj:** +% glm_map with .design populated from SPM, level set to 1, and +% is_timeseries set true. +% +% :See also: +% - fmri_glm_design_matrix, build_design +% +% .. +% Programmers' notes: +% SCAFFOLD - not yet implemented. Planned field mapping (SPM -> object): +% SPM.Sess -> obj.design.Sess (onsets .ons, durations .dur, +% names .name, parametric mods .P, covariates .C) +% SPM.xBF -> obj.design.xBF (basis set, T, T0, UNITS, Volterra) +% SPM.xX -> obj.design.xX (design matrix .X, partitions, names) +% SPM.xY.RT -> obj.design.TR +% SPM.nscan -> obj.design.nscan +% obj.level = 1; obj.is_timeseries = true; +% SPM25 deltas vs SPM12 will be reconciled against a real SPM.mat at +% implementation time. +% .. + +error('glm_map:NotImplemented', ... + 'import_SPM() is scaffolded but not yet implemented. Planned for Phase 3.'); + +end % import_SPM diff --git a/CanlabCore/@glm_map/montage.m b/CanlabCore/@glm_map/montage.m new file mode 100644 index 00000000..0b984c46 --- /dev/null +++ b/CanlabCore/@glm_map/montage.m @@ -0,0 +1,57 @@ +function o2 = montage(obj, varargin) +% Display a montage of a chosen result map from a glm_map object. +% +% Selects one of the fitted statistic_image maps (betas, t, contrast +% estimates, or contrast t) and delegates to statistic_image/image_vector +% montage for display. +% +% :Usage: +% :: +% +% montage(obj) % default: thresholded t map +% montage(obj, which_map, ...) % which_map selects the map +% o2 = montage(obj, ...) % return the fmridisplay handle +% +% :Inputs: +% +% **obj:** +% A fitted glm_map object. +% +% :Optional Inputs: +% +% **which_map:** +% One of 'betas' | 't' | 'contrast' | 'contrast_t' (default 't'). +% Any remaining arguments are passed through to statistic_image.montage. +% +% :Outputs: +% +% **o2:** +% The fmridisplay object created by the underlying montage call. +% +% :Examples: +% :: +% +% montage(g, 't'); +% montage(g, 'contrast', 'trans_white'); +% +% :See also: +% - statistic_image, image_vector.montage, table, plot_design +% +% .. +% Programmers' notes: +% 2026 - Initial implementation (thin delegation). +% .. + +[map, which_map, varargin] = select_map(obj, varargin{:}); + +if isempty(map) + error('glm_map:NoMap', 'Requested map ''%s'' is empty. Fit the model first.', which_map); +end + +if nargout > 0 + o2 = montage(map, varargin{:}); +else + montage(map, varargin{:}); +end + +end % montage diff --git a/CanlabCore/@glm_map/plot_design.m b/CanlabCore/@glm_map/plot_design.m new file mode 100644 index 00000000..93986773 --- /dev/null +++ b/CanlabCore/@glm_map/plot_design.m @@ -0,0 +1,74 @@ +function plot_design(obj, varargin) +% Plot the design matrix and per-regressor VIFs for a glm_map object. +% +% Shows the design matrix X as an image (columns = regressors) and, when +% available, a bar plot of the variance inflation factor for each regressor. +% +% :Usage: +% :: +% +% plot_design(obj) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design available (obj.X non-empty). +% +% :Outputs: +% +% A figure with the design matrix and (if computed) VIF bars. +% +% :Examples: +% :: +% +% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2); +% g = diagnostics(g, 'noverbose'); +% plot_design(g); +% +% :See also: +% - diagnostics, fmri_glm_design_matrix.plot +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. +% .. + +X = obj.X; +if isempty(X) + error('glm_map:NoDesign', 'No design matrix available (obj.X is empty).'); +end + +rn = obj.regressor_names; +has_vif = ~isempty(obj.vif); + +create_figure('glm_map design', 1, 1 + has_vif); + +% --- Design matrix image --- +subplot(1, 1 + has_vif, 1); +imagesc(X); +colormap(gray); +colorbar; +xlabel('Regressor'); +ylabel('Image / observation'); +title('Design matrix X'); +if ~isempty(rn) && numel(rn) == size(X, 2) + set(gca, 'XTick', 1:size(X, 2), 'XTickLabel', rn, 'XTickLabelRotation', 45); +end + +% --- VIF bars --- +if has_vif + subplot(1, 2, 2); + bar(obj.vif); + ylabel('Variance inflation factor'); + xlabel('Regressor'); + title('VIF per regressor'); + if ~isempty(rn) && numel(rn) == numel(obj.vif) + set(gca, 'XTick', 1:numel(obj.vif), 'XTickLabel', rn, 'XTickLabelRotation', 45); + end + hold on; + yl = get(gca, 'YLim'); + plot(get(gca, 'XLim'), [4 4], 'r--'); % conventional VIF = 4 reference + set(gca, 'YLim', [0 max(yl(2), 4.5)]); +end + +end % plot_design diff --git a/CanlabCore/@glm_map/private/select_map.m b/CanlabCore/@glm_map/private/select_map.m new file mode 100644 index 00000000..c0bb7a50 --- /dev/null +++ b/CanlabCore/@glm_map/private/select_map.m @@ -0,0 +1,35 @@ +function [map, which_map, varargin] = select_map(obj, varargin) +% Private helper for glm_map: select a result map by keyword. +% +% Pulls an optional leading map-name argument from varargin and returns the +% corresponding statistic_image (or fmri_data) result map stored in the +% glm_map object, the resolved name, and the remaining varargin. +% +% which_map (first arg if char and a recognized name) is one of: +% 'betas' | 't' | 'contrast' | 'contrast_t' +% Default is 't'. + +which_map = 't'; + +valid = {'betas', 't', 'contrast', 'contrast_t'}; + +if ~isempty(varargin) && (ischar(varargin{1}) || isstring(varargin{1})) ... + && any(strcmpi(char(varargin{1}), valid)) + which_map = lower(char(varargin{1})); + varargin(1) = []; +end + +switch which_map + case 'betas' + map = obj.betas; + case 't' + map = obj.t; + case 'contrast' + map = obj.contrast_estimates; + case 'contrast_t' + map = obj.contrast_t; + otherwise + error('glm_map:BadMap', 'Unknown map ''%s''. Use betas | t | contrast | contrast_t.', which_map); +end + +end % select_map diff --git a/CanlabCore/@glm_map/summary.m b/CanlabCore/@glm_map/summary.m new file mode 100644 index 00000000..7c56200c --- /dev/null +++ b/CanlabCore/@glm_map/summary.m @@ -0,0 +1,43 @@ +function summary(obj) +% Print a one-screen summary of a glm_map object. +% +% Shows the design summary (via disp) plus any accumulated warnings and the +% provenance history. +% +% :Usage: +% :: +% +% summary(obj) +% +% :Inputs: +% +% **obj:** +% A glm_map object. +% +% :See also: +% - glm_map, diagnostics +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. +% .. + +disp(obj); + +if ~isempty(obj.warnings) + fprintf(' Warnings:\n'); + for i = 1:numel(obj.warnings) + fprintf(' - %s\n', obj.warnings{i}); + end + fprintf('\n'); +end + +if ~isempty(obj.history) + fprintf(' History:\n'); + for i = 1:numel(obj.history) + fprintf(' %d. %s\n', i, obj.history{i}); + end + fprintf('\n'); +end + +end % summary diff --git a/CanlabCore/@glm_map/table.m b/CanlabCore/@glm_map/table.m new file mode 100644 index 00000000..37451b1e --- /dev/null +++ b/CanlabCore/@glm_map/table.m @@ -0,0 +1,50 @@ +function varargout = table(obj, varargin) +% Print a table of significant regions for a chosen glm_map result map. +% +% Selects one of the fitted statistic_image maps and delegates to +% statistic_image.table (which forms regions and labels them against an atlas). +% +% :Usage: +% :: +% +% table(obj) % default: thresholded t map +% table(obj, which_map, ...) % which_map selects the map +% [results_table, ...] = table(obj, ...) +% +% :Inputs: +% +% **obj:** +% A fitted glm_map object. +% +% :Optional Inputs: +% +% **which_map:** +% One of 'betas' | 't' | 'contrast' | 'contrast_t' (default 't'). +% Remaining arguments pass through to statistic_image.table. +% +% :Outputs: +% +% Whatever statistic_image.table returns (e.g. a results table object). +% +% :Examples: +% :: +% +% table(g, 'contrast'); +% +% :See also: +% - statistic_image.table, region, montage +% +% .. +% Programmers' notes: +% 2026 - Initial implementation (thin delegation). +% .. + +[map, which_map, varargin] = select_map(obj, varargin{:}); + +if isempty(map) + error('glm_map:NoMap', 'Requested map ''%s'' is empty. Fit the model first.', which_map); +end + +[varargout{1:nargout}] = table(map, varargin{:}); + +end % table diff --git a/CanlabCore/@glm_map/threshold.m b/CanlabCore/@glm_map/threshold.m new file mode 100644 index 00000000..dc5c5ddc --- /dev/null +++ b/CanlabCore/@glm_map/threshold.m @@ -0,0 +1,71 @@ +function obj = threshold(obj, varargin) +% Re-threshold the statistic maps of a glm_map object (no refitting). +% +% Applies a statistical threshold to the fitted t and/or contrast_t maps by +% delegating to statistic_image.threshold. The underlying statistic values +% are preserved; only the significance mask changes. +% +% :Usage: +% :: +% +% obj = threshold(obj, pval, thresh_type, varargin) +% +% :Inputs: +% +% **obj:** +% A fitted glm_map object (obj.is_fitted == true). +% +% Remaining inputs are passed through unchanged to +% statistic_image.threshold, e.g. (.001, 'unc', 'k', 10) or (.05, 'fdr'). +% +% :Optional Inputs: +% +% **'which_map', 'contrast' | 't' | 'both':** +% Which map(s) to threshold. Default 'both' (t and contrast_t when +% present). +% +% :Outputs: +% +% **obj:** +% glm_map with the selected statistic_image maps re-thresholded. +% +% :Examples: +% :: +% +% g = threshold(g, .001, 'unc', 'k', 10); +% g = threshold(g, .05, 'fdr', 'which_map', 'contrast'); +% +% :See also: +% - statistic_image.threshold +% +% .. +% Programmers' notes: +% 2026 - Initial implementation (thin delegation). +% .. + +if ~obj.is_fitted + error('glm_map:NotFitted', 'Object is not fitted; run fit(obj, data) before thresholding.'); +end + +% Pull out the which_map keyword if present +which_map = 'both'; +wh = find(strcmpi(varargin, 'which_map')); +if ~isempty(wh) + which_map = varargin{wh(1) + 1}; + varargin(wh(1):wh(1) + 1) = []; +end + +do_t = ismember(lower(which_map), {'t', 'both'}); +do_con = ismember(lower(which_map), {'contrast', 'both'}); + +if do_t && ~isempty(obj.t) + obj.t = threshold(obj.t, varargin{:}); +end + +if do_con && ~isempty(obj.contrast_t) + obj.contrast_t = threshold(obj.contrast_t, varargin{:}); +end + +obj.history{end + 1} = 'threshold: re-thresholded statistic maps'; + +end % threshold From 77acea3da5bb447d27838ae84151746db8e063c8 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Wed, 17 Jun 2026 19:30:06 -0400 Subject: [PATCH 2/8] glm_map: implement fit() over fmri_data.regress fit() is now a working orchestration layer: validates the design, attaches obj.X to the fmri_data, assembles regress() arguments (robust/AR/threshold/contrasts/names), runs regress, and unpacks the outputs into the object's statistic_image maps (betas, t, contrast_estimates, contrast_t) plus sigma, dfe, residuals. Records fit_parameters and runs diagnostics() for the full VIF/cVIF/leverage/ collinearity set. AR models are gated on is_timeseries. Also resolve multi-image result maps for display: select_map now parses an optional image index (bare scalar or 'wh_image'), montage selects it when given, and table auto-selects a single image (required because statistic_image.table handles one image at a time). Verified end-to-end on load_image_set('emotionreg'): fit -> dfe=28, betas [35676x2], diagnostics populated, threshold and atlas-labeled table() both working. Co-Authored-By: Claude Opus 4.8 (1M context) --- CanlabCore/@glm_map/fit.m | 208 ++++++++++++++++++++--- CanlabCore/@glm_map/montage.m | 7 +- CanlabCore/@glm_map/private/select_map.m | 33 +++- CanlabCore/@glm_map/table.m | 19 ++- 4 files changed, 236 insertions(+), 31 deletions(-) diff --git a/CanlabCore/@glm_map/fit.m b/CanlabCore/@glm_map/fit.m index 79f20af8..15f89aeb 100644 --- a/CanlabCore/@glm_map/fit.m +++ b/CanlabCore/@glm_map/fit.m @@ -1,10 +1,10 @@ function obj = fit(obj, data, varargin) % Fit a mass-univariate GLM and populate result maps in a glm_map object. % -% Scikit-learn-style fit: builds the design (if event/1st-level), runs -% fmri_data.regress on the data, and unpacks the outputs into the object's -% statistic_image result maps (betas, t, contrast_estimates, contrast_t), -% along with sigma, dfe, and the design diagnostics. +% Scikit-learn-style fit: builds the design (if event/1st-level and not yet +% built), runs fmri_data.regress on the data, and unpacks the outputs into +% the object's statistic_image result maps (betas, t, contrast_estimates, +% contrast_t), plus sigma, dfe, and design diagnostics. % % :Usage: % :: @@ -15,11 +15,12 @@ % % **obj:** % A glm_map object with a design specified (either a .design event -% model or a direct .X matrix). +% model or a direct .X matrix), and optionally contrasts. % % **data:** -% An fmri_data object whose images correspond to the rows of X -% (timeseries for level 1, contrast/subject images for level 2). +% An fmri_data object whose images (columns of .dat) correspond to the +% rows of X (a within-run timeseries for level 1, contrast/subject +% images for level 2). % % :Optional Inputs: % @@ -30,9 +31,19 @@ % Use an autoregressive error model of the given order. Only valid % when obj.is_timeseries is true. % +% **'pthresh', [p-value]:** +% Initial threshold p-value applied to the t/contrast_t maps +% (default 0.001). +% +% **'thresh_type', 'unc' | 'fdr':** +% Initial threshold type (default 'unc'). +% % **'residuals':** % Also store residuals in obj.residuals. % +% **'display':** +% Show fmri_data.regress orthviews (default off; glm_map suppresses it). +% % **'doverbose' / 'noverbose':** % Toggle verbose output (default true). % @@ -46,28 +57,181 @@ % :: % % dat = load_image_set('emotionreg'); -% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2); +% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2, ... +% 'regressor_names', {'intercept' 'cov'}); +% g = add_contrasts(g, [0 1], {'cov_effect'}); % g = fit(g, dat); +% table(g, 'contrast'); montage(g, 'contrast_t'); % % :See also: -% - fmri_data.regress, build_design, diagnostics +% - fmri_data.regress, build_design, diagnostics, threshold % % .. % Programmers' notes: -% SCAFFOLD - not yet implemented. Planned behavior: -% 1. If event/1st-level and X not built, call build_design(obj). -% 2. Attach obj.X to data.X (and obj.contrasts/contrast_names). -% 3. If obj.is_timeseries and 'AR' requested, pass 'AR' to regress; -% otherwise error if 'AR' requested on non-timeseries data. -% 4. out = regress(data, ...): unpack out.b -> obj.betas, out.t -> obj.t, -% out.contrast_images -> obj.contrast_estimates, out.con_t -> -% obj.contrast_t, out.sigma -> obj.sigma, out.df -> obj.dfe. -% 5. Populate obj.vif/contrast_vif/leverages/warnings from out.diagnostics -% (or call diagnostics(obj)). -% 6. Record options in obj.fit_parameters and append to obj.history. +% 2026 - Initial implementation. Thin orchestration layer over +% fmri_data.regress: glm_map owns the design/diagnostics/result container, +% regress remains the compute engine. % .. -error('glm_map:NotImplemented', ... - 'fit() is scaffolded but not yet implemented. Planned for Phase 2 (wraps fmri_data.regress).'); +% ------------------------------------------------------------------------- +% Parse options +% ------------------------------------------------------------------------- +doverbose = ~any(strcmpi(varargin, 'noverbose')); +do_robust = any(strcmpi(varargin, 'robust')); +do_resid = any(strcmpi(varargin, 'residuals')) || any(strcmpi(varargin, 'residual')); +do_display = any(strcmpi(varargin, 'display')); + +pthresh = 0.001; +wh = find(strcmpi(varargin, 'pthresh')); +if ~isempty(wh), pthresh = varargin{wh(1) + 1}; end + +thresh_type = 'unc'; +wh = find(strcmpi(varargin, 'thresh_type')); +if ~isempty(wh), thresh_type = varargin{wh(1) + 1}; end +if ~ismember(lower(thresh_type), {'unc', 'fdr'}) + error('glm_map:BadThreshType', 'thresh_type must be ''unc'' or ''fdr''.'); +end + +ar_order = 0; +wh = find(strcmpi(varargin, 'AR')); +if ~isempty(wh), ar_order = varargin{wh(1) + 1}; end + +% ------------------------------------------------------------------------- +% Validate inputs +% ------------------------------------------------------------------------- +if ~isa(data, 'fmri_data') + error('glm_map:BadData', 'data must be an fmri_data object.'); +end + +if ar_order > 0 && ~obj.is_timeseries + error('glm_map:ARnotTimeseries', ... + 'AR error models require obj.is_timeseries == true (within-run timeseries data).'); +end + +% Ensure a design matrix is available; build it for event/1st-level models +if isempty(obj.X) + if obj.level == 1 && ~isempty(obj.design) + obj = build_design(obj); + else + error('glm_map:NoDesign', ... + 'No design available. Supply obj.X (direct mode) or an obj.design to build (event mode).'); + end +end + +X = obj.X; + +if size(data.dat, 2) ~= size(X, 1) + error('glm_map:SizeMismatch', ... + 'data has %d images but the design has %d rows. They must match.', size(data.dat, 2), size(X, 1)); +end + +% ------------------------------------------------------------------------- +% Assemble regress() argument list +% ------------------------------------------------------------------------- +data.X = X; + +% Threshold p-value must precede the 'unc'/'fdr' keyword (regress convention) +regress_args = {pthresh, thresh_type}; + +if do_robust, regress_args{end + 1} = 'robust'; end +if ar_order > 0, regress_args = [regress_args, {'AR', ar_order}]; end +if do_resid, regress_args{end + 1} = 'residual'; end + +% Suppress regress's own orthviews unless explicitly requested +if ~do_display, regress_args{end + 1} = 'nodisplay'; end +if ~doverbose, regress_args{end + 1} = 'noverbose'; end + +% Regressor names (length must match number of design columns) +rn = obj.regressor_names; +if ~isempty(rn) && numel(rn) == size(X, 2) + regress_args = [regress_args, {'names', rn}]; +end + +% Contrasts: regress expects C as [regressors x contrasts] (our storage) +if ~isempty(obj.contrasts) + regress_args = [regress_args, {'C', obj.contrasts}]; + if ~isempty(obj.contrast_names) + regress_args = [regress_args, {'contrast_names', obj.contrast_names}]; + end +end + +if ~isempty(obj.analysis_name) + regress_args = [regress_args, {'analysis_name', obj.analysis_name}]; +end + +% ------------------------------------------------------------------------- +% Run the regression +% ------------------------------------------------------------------------- +out = regress(data, regress_args{:}); + +% ------------------------------------------------------------------------- +% Unpack result maps +% ------------------------------------------------------------------------- +obj.betas = out.b; +obj.t = out.t; +obj.sigma = out.sigma; + +% Error degrees of freedom (scalar summary; per-voxel df lives in out.df) +if isfield(out, 'df') && ~isempty(out.df) && ~isempty(out.df.dat) + obj.dfe = double(median(out.df.dat(:))); +else + obj.dfe = size(X, 1) - size(X, 2); +end +obj.t.dfe = obj.dfe; + +if isfield(out, 'contrast_images') && ~isempty(out.contrast_images) + obj.contrast_estimates = out.contrast_images; + obj.contrast_t = out.con_t; + obj.contrast_t.dfe = obj.dfe; +end + +if do_resid && isfield(out, 'resid') + obj.residuals = out.resid; +end + +% Sync names back from regress (it may have generated/added them) +if isempty(obj.design) + obj.regressor_names_direct = out.variable_names; +end +if ~isempty(out.contrast_names) + obj.contrast_names = out.contrast_names; +end + +% ------------------------------------------------------------------------- +% Diagnostics and provenance +% ------------------------------------------------------------------------- +obj.warnings = [obj.warnings(:); out.warnings(:)]'; + +obj.fit_parameters = struct( ... + 'robust', do_robust, ... + 'ar_order', ar_order, ... + 'is_timeseries', obj.is_timeseries, ... + 'pthresh', pthresh, ... + 'thresh_type', thresh_type, ... + 'do_resid', do_resid); + +% Compute the full diagnostic set (adds cVIF, condition number, collinearity +% report; uses canonical VIF/cVIF rather than regress's getvif) +obj = diagnostics(obj, 'noverbose'); + +obj.history{end + 1} = sprintf('fit: regress (%s%s), p<%g %s, dfe=%g', ... + local_iif(do_robust, 'robust', 'OLS'), ... + local_iif(ar_order > 0, sprintf(', AR(%d)', ar_order), ''), ... + pthresh, thresh_type, obj.dfe); + +if doverbose + fprintf('\n glm_map fit complete: %d regressors, %d contrasts, dfe = %g.\n', ... + obj.num_regressors, obj.num_contrasts, obj.dfe); + if ~isempty(obj.warnings) + fprintf(' %d warning(s); see obj.warnings.\n', numel(obj.warnings)); + end +end end % fit + + +% ===================================================================== +function s = local_iif(tf, a, b) +% Inline if: return a if tf else b. +if tf, s = a; else, s = b; end +end % local_iif diff --git a/CanlabCore/@glm_map/montage.m b/CanlabCore/@glm_map/montage.m index 0b984c46..e7b1e593 100644 --- a/CanlabCore/@glm_map/montage.m +++ b/CanlabCore/@glm_map/montage.m @@ -42,12 +42,17 @@ % 2026 - Initial implementation (thin delegation). % .. -[map, which_map, varargin] = select_map(obj, varargin{:}); +[map, which_map, wh_image, varargin] = select_map(obj, varargin{:}); if isempty(map) error('glm_map:NoMap', 'Requested map ''%s'' is empty. Fit the model first.', which_map); end +% Select a single image if requested (montage can also render multi-image maps) +if ~isempty(wh_image) + map = get_wh_image(map, wh_image); +end + if nargout > 0 o2 = montage(map, varargin{:}); else diff --git a/CanlabCore/@glm_map/private/select_map.m b/CanlabCore/@glm_map/private/select_map.m index c0bb7a50..e61f95f8 100644 --- a/CanlabCore/@glm_map/private/select_map.m +++ b/CanlabCore/@glm_map/private/select_map.m @@ -1,24 +1,43 @@ -function [map, which_map, varargin] = select_map(obj, varargin) -% Private helper for glm_map: select a result map by keyword. +function [map, which_map, wh_image, varargin] = select_map(obj, varargin) +% Private helper for glm_map: select a result map (and optional image) by keyword. % -% Pulls an optional leading map-name argument from varargin and returns the -% corresponding statistic_image (or fmri_data) result map stored in the -% glm_map object, the resolved name, and the remaining varargin. +% Pulls an optional leading map-name argument and an optional image index from +% varargin, and returns the corresponding statistic_image result map stored in +% the glm_map object, the resolved map name, the requested image index (or [] +% if none given), and the remaining varargin. % % which_map (first arg if char and a recognized name) is one of: -% 'betas' | 't' | 'contrast' | 'contrast_t' -% Default is 't'. +% 'betas' | 't' | 'contrast' | 'contrast_t' (default 't') +% +% The image index may be given either as a bare numeric scalar following the +% map name, or as a 'wh_image'/'image' keyword-value pair. For betas/t this +% indexes regressors; for contrast/contrast_t it indexes contrasts. which_map = 't'; +wh_image = []; valid = {'betas', 't', 'contrast', 'contrast_t'}; +% Optional leading map name if ~isempty(varargin) && (ischar(varargin{1}) || isstring(varargin{1})) ... && any(strcmpi(char(varargin{1}), valid)) which_map = lower(char(varargin{1})); varargin(1) = []; end +% Optional image index as a bare numeric scalar +if ~isempty(varargin) && isnumeric(varargin{1}) && isscalar(varargin{1}) + wh_image = varargin{1}; + varargin(1) = []; +end + +% Optional image index as a keyword-value pair +wh = find(strcmpi(varargin, 'wh_image') | strcmpi(varargin, 'image')); +if ~isempty(wh) + wh_image = varargin{wh(1) + 1}; + varargin(wh(1):wh(1) + 1) = []; +end + switch which_map case 'betas' map = obj.betas; diff --git a/CanlabCore/@glm_map/table.m b/CanlabCore/@glm_map/table.m index 37451b1e..73ecd549 100644 --- a/CanlabCore/@glm_map/table.m +++ b/CanlabCore/@glm_map/table.m @@ -39,12 +39,29 @@ % 2026 - Initial implementation (thin delegation). % .. -[map, which_map, varargin] = select_map(obj, varargin{:}); +[map, which_map, wh_image, varargin] = select_map(obj, varargin{:}); if isempty(map) error('glm_map:NoMap', 'Requested map ''%s'' is empty. Fit the model first.', which_map); end +% statistic_image.table requires a single image; resolve one. +nimg = size(map.dat, 2); +if nimg > 1 + if isempty(wh_image), wh_image = 1; end + labels = map.image_labels; + if iscell(labels) && numel(labels) >= wh_image && ~isempty(labels{wh_image}) + thislabel = labels{wh_image}; + else + thislabel = sprintf('image %d', wh_image); + end + fprintf(' %s map has %d images; showing %s (%s). Pass an index to choose another, e.g. table(obj, ''%s'', k).\n', ... + which_map, nimg, num2str(wh_image), thislabel, which_map); + map = get_wh_image(map, wh_image); +elseif ~isempty(wh_image) + map = get_wh_image(map, wh_image); +end + [varargout{1:nargout}] = table(map, varargin{:}); end % table From 5c250be90c8174e451f721ecfe3c43b2dc6dcd84 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Wed, 17 Jun 2026 19:55:50 -0400 Subject: [PATCH 3/8] glm_map: implement build_design() for event/1st-level designs build_design() delegates to fmri_glm_design_matrix.build, which convolves onsets/durations with the basis set and assembles design.xX.X. After the call the Dependent obj.X and obj.regressor_names read through to the built design, so event-mode fit() works the same as direct mode. Also harden the Dependent X / regressor_names accessors: the wrapped design seeds xX as a 0x0 struct array (constructor uses 'name',{}), so guard with isscalar() before indexing .X / .name to avoid an "isempty: not enough input arguments" error pre-build. Verified: glm_map(fmri_glm_design_matrix) -> build_design -> X [200x19] from onsets, regressor names read through, event-mode fit() completes. Co-Authored-By: Claude Opus 4.8 (1M context) --- CanlabCore/@glm_map/build_design.m | 63 +++++++++++++++++++++++------- CanlabCore/@glm_map/glm_map.m | 6 +-- 2 files changed, 52 insertions(+), 17 deletions(-) diff --git a/CanlabCore/@glm_map/build_design.m b/CanlabCore/@glm_map/build_design.m index 716086c0..9520e7d0 100644 --- a/CanlabCore/@glm_map/build_design.m +++ b/CanlabCore/@glm_map/build_design.m @@ -1,10 +1,11 @@ function obj = build_design(obj, varargin) % Build the design matrix X for an event/1st-level glm_map from onsets. % -% Delegates to the wrapped fmri_glm_design_matrix object (and ultimately -% onsets2fmridesign) to convolve onsets/durations with the basis set and -% assemble X, regressor names, and any nuisance covariates. After this call, -% the Dependent obj.X and obj.regressor_names read through to the built design. +% Delegates to the wrapped fmri_glm_design_matrix object's build method, +% which convolves onsets/durations with the basis set and assembles the +% design matrix (interest, covariate, and baseline partitions) into +% design.xX.X with column names in design.xX.name. After this call, the +% Dependent obj.X and obj.regressor_names read through to the built design. % % :Usage: % :: @@ -14,36 +15,70 @@ % :Inputs: % % **obj:** -% A glm_map object with a non-empty .design (fmri_glm_design_matrix). +% A glm_map object with a non-empty .design (fmri_glm_design_matrix) +% that has conditions/onsets assigned. % % :Optional Inputs: % -% **'doplot' / 'noplot':** -% Plot the resulting design matrix (default false). +% **'doplot' / 'plot':** +% Plot the resulting design via fmri_glm_design_matrix.plot (default off). % % :Outputs: % % **obj:** % The glm_map object with the wrapped design built (design.xX.X populated). % +% :Examples: +% :: +% +% TR = 2; nscan = 200; +% ons = {[10 40 70 100]', [25 55 85 115]'}; % 2 conditions, 1 session +% d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... +% 'onsets', ons, 'condition_names', {'A' 'B'}); +% g = glm_map(d); +% g = build_design(g); +% size(g.X) % built design matrix +% g.regressor_names +% % :See also: % - fmri_glm_design_matrix, fmri_glm_design_matrix.build, onsets2fmridesign % % .. % Programmers' notes: -% SCAFFOLD - not yet implemented. Planned behavior: -% 1. Require obj.level == 1 and ~isempty(obj.design). -% 2. obj.design = build(obj.design); % existing method -% 3. Sync obj.contrasts/contrast_names defaults if requested. -% 4. Append to obj.history. +% 2026 - Initial implementation (delegates to fmri_glm_design_matrix.build). % .. +% ------------------------------------------------------------------------- +% Parse options +% ------------------------------------------------------------------------- +doplot = any(strcmpi(varargin, 'doplot')) || any(strcmpi(varargin, 'plot')); + +% ------------------------------------------------------------------------- +% Validate +% ------------------------------------------------------------------------- if isempty(obj.design) error('glm_map:NoDesign', ... 'build_design requires an fmri_glm_design_matrix in obj.design (event/1st-level mode).'); end -error('glm_map:NotImplemented', ... - 'build_design() is scaffolded but not yet implemented. Planned for Phase 3 (wraps fmri_glm_design_matrix.build).'); +if isempty(obj.design.Sess) || isempty(obj.design.Sess(1).U) || isempty(obj.design.Sess(1).U(1).name) + error('glm_map:NoOnsets', ... + ['The wrapped design has no conditions/onsets assigned. Add them first, e.g. ' ... + 'fmri_glm_design_matrix(TR, ''nscan'', nscan, ''onsets'', ons, ''condition_names'', names).']); +end + +% ------------------------------------------------------------------------- +% Build +% ------------------------------------------------------------------------- +obj.design = build(obj.design); + +obj.level = 1; % building implies a first-level/event model + +if doplot + plot(obj.design); +end + +obj.history{end + 1} = sprintf('build_design: built design matrix [%d x %d] from onsets', ... + size(obj.X, 1), size(obj.X, 2)); end % build_design diff --git a/CanlabCore/@glm_map/glm_map.m b/CanlabCore/@glm_map/glm_map.m index e63610af..2e1313b5 100644 --- a/CanlabCore/@glm_map/glm_map.m +++ b/CanlabCore/@glm_map/glm_map.m @@ -284,7 +284,7 @@ % Prefer a built design matrix from the wrapped design object; % fall back to the direct-mode matrix. val = []; - if ~isempty(obj.design) && isstruct(obj.design.xX) ... + if ~isempty(obj.design) && isstruct(obj.design.xX) && isscalar(obj.design.xX) ... && isfield(obj.design.xX, 'X') && ~isempty(obj.design.xX.X) val = obj.design.xX.X; elseif ~isempty(obj.Xdirect) @@ -294,7 +294,7 @@ function val = get.regressor_names(obj) val = {}; - if ~isempty(obj.design) && isstruct(obj.design.xX) ... + if ~isempty(obj.design) && isstruct(obj.design.xX) && isscalar(obj.design.xX) ... && isfield(obj.design.xX, 'name') && ~isempty(obj.design.xX.name) val = obj.design.xX.name; elseif ~isempty(obj.regressor_names_direct) @@ -344,7 +344,7 @@ function obj = set.X(obj, val) % Setting X targets the direct/group design backing store. - if ~isempty(obj.design) && isstruct(obj.design.xX) ... + if ~isempty(obj.design) && isstruct(obj.design.xX) && isscalar(obj.design.xX) ... && isfield(obj.design.xX, 'X') && ~isempty(obj.design.xX.X) warning('glm_map:DesignPresent', ... ['This glm_map has a built event design in .design; setting .X stores a direct ' ... From 6cbfc466ffb053caeb0f8a723d44a13406582ba6 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Wed, 17 Jun 2026 19:59:29 -0400 Subject: [PATCH 4/8] glm_map: implement import_SPM() for SPM12/SPM25 first-level models import_SPM() populates a glm_map from an SPM structure, an SPM.mat path, or a directory containing one. Because fmri_glm_design_matrix mirrors SPM's schema, the import is a guarded copy of substructs into the wrapped design: xY.RT -> TR, nscan, xBF, Sess (onsets/durations/names/pmods/ covariates), and xX (matrix + names). Sets level=1 and is_timeseries=true. Optional 'load_betas' reads beta_*.nii/.img from SPM.swd into obj.betas, labeled by xX.name. Whole-substruct copy is resilient to SPM12 vs SPM25 auxiliary-field differences. Verified on a synthetic SPM struct: TR/onsets/condition_names/regressor names read through correctly, and fit() against a matching timeseries yields dfe = nscan - nregressors. Co-Authored-By: Claude Opus 4.8 (1M context) --- CanlabCore/@glm_map/import_SPM.m | 182 +++++++++++++++++++++++++++---- 1 file changed, 162 insertions(+), 20 deletions(-) diff --git a/CanlabCore/@glm_map/import_SPM.m b/CanlabCore/@glm_map/import_SPM.m index fe609cd9..c72f8c74 100644 --- a/CanlabCore/@glm_map/import_SPM.m +++ b/CanlabCore/@glm_map/import_SPM.m @@ -1,15 +1,18 @@ function obj = import_SPM(obj, SPM, varargin) -% Import a first-level design (and optional betas) from an SPM model. +% Import a first-level design (and optionally betas) from an SPM model. % % Populates a glm_map object from an SPM structure or an SPM.mat path, % mapping SPM's design fields into the wrapped fmri_glm_design_matrix object. -% Compatible with SPM12 and SPM25 first-level models. +% fmri_glm_design_matrix deliberately mirrors SPM's field schema +% (xY/nscan/xBF/Sess/xX), so the import is largely a guarded copy of those +% substructures. Compatible with SPM12 and SPM25 first-level models. % % :Usage: % :: % -% obj = import_SPM(glm_map, SPM) % SPM struct already loaded -% obj = import_SPM(glm_map, '/path/SPM.mat') +% obj = import_SPM(glm_map, SPM) % SPM struct already loaded +% obj = import_SPM(glm_map, '/path/SPM.mat') % path to SPM.mat +% obj = import_SPM(glm_map, '/path/to/dir') % directory containing SPM.mat % % :Inputs: % @@ -17,13 +20,17 @@ % A glm_map object (typically empty: glm_map). % % **SPM:** -% Either an SPM structure (as loaded from SPM.mat) or a char/string -% path to an SPM.mat file. +% Either an SPM structure (as loaded from SPM.mat), a path to an +% SPM.mat file, or a directory containing one. % % :Optional Inputs: % % **'load_betas':** -% Also load the estimated beta_*.nii images into obj.betas. +% Also load the estimated beta_*.nii/.img images into obj.betas +% (looked up in SPM.swd), labeled by SPM.xX.name. +% +% **'noverbose':** +% Suppress console output. % % :Outputs: % @@ -31,24 +38,159 @@ % glm_map with .design populated from SPM, level set to 1, and % is_timeseries set true. % +% :Examples: +% :: +% +% g = import_SPM(glm_map, '/data/sub-01/1stlevel/SPM.mat'); +% g.onsets % read-through to the imported design +% g = add_contrasts(g, c, names); +% % (data still needed to fit: g = fit(g, fmri_timeseries_obj)) +% % :See also: -% - fmri_glm_design_matrix, build_design +% - fmri_glm_design_matrix, build_design, fit % % .. % Programmers' notes: -% SCAFFOLD - not yet implemented. Planned field mapping (SPM -> object): -% SPM.Sess -> obj.design.Sess (onsets .ons, durations .dur, -% names .name, parametric mods .P, covariates .C) -% SPM.xBF -> obj.design.xBF (basis set, T, T0, UNITS, Volterra) -% SPM.xX -> obj.design.xX (design matrix .X, partitions, names) -% SPM.xY.RT -> obj.design.TR -% SPM.nscan -> obj.design.nscan -% obj.level = 1; obj.is_timeseries = true; -% SPM25 deltas vs SPM12 will be reconciled against a real SPM.mat at -% implementation time. +% 2026 - Initial implementation. Field mapping (SPM -> wrapped design): +% SPM.xY.RT -> design.TR / design.xY.RT +% SPM.nscan -> design.nscan +% SPM.xBF -> design.xBF (basis set, T, T0, UNITS, Volterra) +% SPM.Sess -> design.Sess (onsets .ons, durations .dur, names .name, +% parametric mods .P, covariates .C) +% SPM.xX -> design.xX (design matrix .X, partition idx, names) +% SPM12 and SPM25 share these first-level fields; whole-substruct copy is +% resilient to extra/renamed auxiliary fields (W, K, xKXs, pKX, ...). % .. -error('glm_map:NotImplemented', ... - 'import_SPM() is scaffolded but not yet implemented. Planned for Phase 3.'); +% ------------------------------------------------------------------------- +% Parse options +% ------------------------------------------------------------------------- +doverbose = ~any(strcmpi(varargin, 'noverbose')); +do_load_betas = any(strcmpi(varargin, 'load_betas')); + +% ------------------------------------------------------------------------- +% Resolve SPM (struct, file path, or directory) +% ------------------------------------------------------------------------- +if ischar(SPM) || isstring(SPM) + spmpath = char(SPM); + if isfolder(spmpath) + spmpath = fullfile(spmpath, 'SPM.mat'); + end + if ~exist(spmpath, 'file') + error('glm_map:NoSPMfile', 'Could not find SPM.mat at: %s', spmpath); + end + loaded = load(spmpath, 'SPM'); + if ~isfield(loaded, 'SPM') + error('glm_map:BadSPMfile', 'File does not contain an SPM variable: %s', spmpath); + end + SPM = loaded.SPM; +end + +if ~isstruct(SPM) || ~isfield(SPM, 'xX') + error('glm_map:BadSPM', 'SPM must be an SPM structure (or path to one) with an xX field.'); +end + +% ------------------------------------------------------------------------- +% Determine TR +% ------------------------------------------------------------------------- +TR = NaN; +if isfield(SPM, 'xY') && isfield(SPM.xY, 'RT') + TR = SPM.xY.RT; +elseif isfield(SPM, 'xBF') && isfield(SPM.xBF, 'dt') && isfield(SPM.xBF, 'T') + TR = SPM.xBF.dt * SPM.xBF.T; +end + +% ------------------------------------------------------------------------- +% Build the wrapped design object and copy SPM substructures +% ------------------------------------------------------------------------- +design = fmri_glm_design_matrix(TR); + +if isfield(SPM, 'xY'), design.xY = SPM.xY; end +if isfield(SPM, 'nscan'), design.nscan = SPM.nscan; end +if isfield(SPM, 'xBF'), design.xBF = SPM.xBF; end +if isfield(SPM, 'Sess'), design.Sess = SPM.Sess; end +design.xX = SPM.xX; % includes .X and .name (read by the Dependent accessors) + +design.build_method = 'Imported from SPM'; +if ~iscell(design.history), design.history = {}; end +design.history{end + 1} = 'Imported from SPM structure'; + +% ------------------------------------------------------------------------- +% Populate the glm_map object +% ------------------------------------------------------------------------- +obj.design = design; +obj.level = 1; +obj.is_timeseries = true; + +obj.history{end + 1} = sprintf('import_SPM: imported %d-session, %d-column design from SPM', ... + local_num_sess(SPM), size(obj.X, 2)); + +% ------------------------------------------------------------------------- +% Optionally load estimated betas +% ------------------------------------------------------------------------- +if do_load_betas + obj = local_load_betas(obj, SPM, doverbose); +end + +% ------------------------------------------------------------------------- +% Report +% ------------------------------------------------------------------------- +if doverbose + fprintf(' import_SPM: TR=%g, %d regressors, %d session(s).%s\n', ... + TR, size(obj.X, 2), local_num_sess(SPM), ... + local_iif(do_load_betas && ~isempty(obj.betas), ' Betas loaded.', '')); +end end % import_SPM + + +% ===================================================================== +% Local helpers +% ===================================================================== +function n = local_num_sess(SPM) +if isfield(SPM, 'Sess') && ~isempty(SPM.Sess) + n = numel(SPM.Sess); +else + n = 0; +end +end % local_num_sess + + +function obj = local_load_betas(obj, SPM, doverbose) +% Load beta_*.nii / beta_*.img images from the SPM working directory. + +swd = ''; +if isfield(SPM, 'swd') && ~isempty(SPM.swd), swd = SPM.swd; end +if isempty(swd) || ~isfolder(swd) + warning('glm_map:NoSWD', 'SPM.swd is not a valid folder; cannot load betas.'); + return +end + +files = dir(fullfile(swd, 'beta_*.nii')); +if isempty(files), files = dir(fullfile(swd, 'beta_*.img')); end +if isempty(files) + warning('glm_map:NoBetas', 'No beta_*.nii/.img images found in %s.', swd); + return +end + +% Sort numerically by the beta index to keep column order aligned with X +names = {files.name}; +fullnames = cellfun(@(f) fullfile(swd, f), names(:), 'UniformOutput', false); +betas = statistic_image(char(fullnames), 'type', 'Beta', 'noverbose'); + +if isfield(SPM, 'xX') && isfield(SPM.xX, 'name') + betas.image_labels = SPM.xX.name; +end + +obj.betas = betas; + +if doverbose + fprintf(' Loaded %d beta image(s) from %s.\n', numel(fullnames), swd); +end + +end % local_load_betas + + +function s = local_iif(tf, a, b) +if tf, s = a; else, s = b; end +end % local_iif From 10257647e9935607bc12a79cc201e5ea18fd9638 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Thu, 18 Jun 2026 17:37:58 -0400 Subject: [PATCH 5/8] Document glm_map in CLAUDE.md and image_vector see-also Add glm_map to the object-architecture class list in CLAUDE.md and to the image_vector :See also: block. Co-Authored-By: Claude Opus 4.8 (1M context) --- CLAUDE.md | 1 + CanlabCore/@image_vector/image_vector.m | 1 + 2 files changed, 2 insertions(+) diff --git a/CLAUDE.md b/CLAUDE.md index 3bfb7dd0..74f5b4a2 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -38,6 +38,7 @@ Other top-level classes (not subclasses of `image_vector`): - **`brainpathway` / `brainpathway_multisubject`** — connectivity / pathway-modeling objects. - **`canlab_dataset`** — generic subject × variable behavioral/clinical data container with its own `glm`, `mediation`, `scatterplot`, etc. - **`fmri_glm_design_matrix`**, **`fmri_timeseries`**, **`predictive_model`** — specialized containers for design matrices, raw timeseries, and ML model artifacts. +- **`glm_map`** — scikit-learn-style estimator for mass-univariate GLM / multiple regression. Bundles the design specification (wrapping an `fmri_glm_design_matrix` in `.design`, or a direct `.X` matrix), the fitted result maps (betas, t, contrast estimates/t as `statistic_image` objects), and design diagnostics (VIF, contrast VIF, leverage, collinearity). Workflow is `g = glm_map(...)`, `g = fit(g, fmri_data_obj)`, then `diagnostics`, `table`, `montage`. Supports 1st-level event designs (`build_design`, `import_SPM`) and 2nd-level/group designs; `fmri_data.regress` is the compute engine. ### MATLAB `@class` directories diff --git a/CanlabCore/@image_vector/image_vector.m b/CanlabCore/@image_vector/image_vector.m index 35c1874f..f519bcaa 100644 --- a/CanlabCore/@image_vector/image_vector.m +++ b/CanlabCore/@image_vector/image_vector.m @@ -227,6 +227,7 @@ % - statistic_image % - atlas % - region +% - glm_map % % .. % Author and copyright information: From 44db014890e837b12413832de950d190913af5e3 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Thu, 18 Jun 2026 23:52:23 -0400 Subject: [PATCH 6/8] Add unit tests for glm_map (canlab_test_glm_map) 12 functiontests covering construction and Dependent property read-through, design diagnostics (incl. rank-deficiency detection), direct-mode (2nd-level) fit over fmri_data.regress with contrasts, threshold delegation, event-mode build_design via the wrapped fmri_glm_design_matrix, import_SPM + event-mode fit on a synthetic SPM struct, and the main input-validation error paths (no design, AR without timeseries, contrast size mismatch). Auto-discovered by canlab_run_all_tests. All 12 pass (~5s). Co-Authored-By: Claude Opus 4.8 (1M context) --- .../Unit_tests/glm_map/canlab_test_glm_map.m | 217 ++++++++++++++++++ 1 file changed, 217 insertions(+) create mode 100644 CanlabCore/Unit_tests/glm_map/canlab_test_glm_map.m diff --git a/CanlabCore/Unit_tests/glm_map/canlab_test_glm_map.m b/CanlabCore/Unit_tests/glm_map/canlab_test_glm_map.m new file mode 100644 index 00000000..2e5d276f --- /dev/null +++ b/CanlabCore/Unit_tests/glm_map/canlab_test_glm_map.m @@ -0,0 +1,217 @@ +function tests = canlab_test_glm_map +%CANLAB_TEST_GLM_MAP Smoke/behavior tests for the glm_map object class. +% +% Covers construction, true Dependent property read-through, design +% diagnostics, direct-mode (2nd-level) fit over fmri_data.regress, event-mode +% build_design via the wrapped fmri_glm_design_matrix, import_SPM, and the +% main input-validation error paths. + +tests = functiontests(localfunctions); +end + + +% ===================================================================== +% Construction and Dependent properties +% ===================================================================== +function test_empty_constructor(tc) +g = glm_map; +tc.verifyClass(g, 'glm_map'); +tc.verifyEqual(g.level, 2); +tc.verifyFalse(g.is_timeseries); +tc.verifyFalse(g.is_fitted); +tc.verifyEqual(g.num_images, 0); +tc.verifyEqual(g.num_regressors, 0); +tc.verifyEqual(g.num_contrasts, 0); +end + + +function test_direct_design_dependent_properties(tc) +X = [ones(20, 1) zscore((1:20)') randn(20, 1)]; +g = glm_map('X', X, 'level', 2, 'regressor_names', {'intercept', 'slope', 'nuis'}); + +tc.verifyEqual(g.num_images, 20); +tc.verifyEqual(g.num_regressors, 3); +tc.verifyEqual(g.X, X); % read-through to Xdirect +tc.verifyEqual(g.regressor_names, {'intercept', 'slope', 'nuis'}); +tc.verifyFalse(g.is_fitted); +end + + +function test_add_contrasts(tc) +X = [ones(20, 1) zscore((1:20)') randn(20, 1)]; +g = glm_map('X', X, 'level', 2); +g = add_contrasts(g, [0 1 0; 0 0 1], {'slope_con', 'nuis_con'}); + +tc.verifyEqual(g.num_contrasts, 2); +tc.verifyEqual(size(g.contrasts), [3 2]); % stored [regressors x contrasts] +tc.verifyEqual(g.contrast_names, {'slope_con', 'nuis_con'}); +end + + +% ===================================================================== +% Diagnostics (design only; no fit required) +% ===================================================================== +function test_diagnostics_on_design(tc) +rng(0); +X = [ones(30, 1) zscore((1:30)') randn(30, 1)]; +g = glm_map('X', X, 'level', 2); +g = add_contrasts(g, [0 1 0], {'slope'}); +g = diagnostics(g, 'noverbose'); + +tc.verifyNumElements(g.vif, 3); +tc.verifyNumElements(g.contrast_vif, 1); +tc.verifyGreaterThan(g.condition_number, 0); +tc.verifyFalse(g.rank_deficient); +tc.verifyTrue(isstruct(g.collinearity_report)); +end + + +function test_diagnostics_flags_rank_deficiency(tc) +X = [ones(20, 1) (1:20)' (1:20)']; % duplicate column +g = glm_map('X', X, 'level', 2); +w = warning('off', 'all'); +c = onCleanup(@() warning(w)); +g = diagnostics(g, 'noverbose'); + +tc.verifyTrue(g.rank_deficient); +tc.verifySize(g.collinearity_report.duplicate_column_pairs, [1 2]); +tc.verifyNotEmpty(g.warnings); +end + + +% ===================================================================== +% Direct-mode (2nd-level / group) fit on real sample data +% ===================================================================== +function test_direct_mode_fit(tc) +dat = canlab_get_sample_fmri_data(); +n = size(dat.dat, 2); +g = glm_map('X', [ones(n, 1) zscore((1:n)')], 'level', 2, ... + 'regressor_names', {'intercept', 'cov'}); +g = add_contrasts(g, [0 1], {'cov_effect'}); +g = fit(g, dat, 'noverbose'); + +tc.verifyTrue(g.is_fitted); +tc.verifyEqual(g.dfe, n - 2); +tc.verifyClass(g.betas, 'statistic_image'); +tc.verifyClass(g.t, 'statistic_image'); +tc.verifyClass(g.contrast_estimates, 'statistic_image'); +tc.verifyClass(g.contrast_t, 'statistic_image'); +tc.verifyEqual(size(g.betas.dat, 2), 2); % one beta image per regressor +tc.verifyEqual(size(g.contrast_estimates.dat, 2), 1); % one image per contrast +tc.verifyNotEmpty(g.vif); +end + + +function test_threshold_returns_glm_map(tc) +dat = canlab_get_sample_fmri_data(); +n = size(dat.dat, 2); +g = glm_map('X', [ones(n, 1) zscore((1:n)')], 'level', 2); +g = fit(g, dat, 'noverbose'); + +w = warning('off', 'all'); +c = onCleanup(@() warning(w)); +g = threshold(g, .01, 'unc', 'which_map', 't'); +tc.verifyClass(g, 'glm_map'); +tc.verifyTrue(g.is_fitted); +end + + +% ===================================================================== +% Event-mode (1st-level) design build via wrapped fmri_glm_design_matrix +% ===================================================================== +function test_event_mode_build_design(tc) +TR = 2; nscan = 120; +ons = {[10 40 70 100]', [25 55 85 115]'}; +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', ons, 'condition_names', {'A', 'B'}); +g = glm_map(d); + +tc.verifyEqual(g.level, 1); % wrapping a design implies 1st-level +tc.verifyEqual(g.TR, TR); % read-through to design.TR +tc.verifyEmpty(g.X); % not built yet + +w = warning('off', 'all'); +c = onCleanup(@() warning(w)); +g = build_design(g); + +tc.verifyEqual(size(g.X, 1), nscan); % built design matrix +tc.verifyGreaterThan(size(g.X, 2), 0); +tc.verifyNotEmpty(g.regressor_names); +tc.verifyNumElements(g.onsets, 2); % onsets read-through to design +end + + +% ===================================================================== +% import_SPM (SPM12/SPM25 first-level) +% ===================================================================== +function test_import_SPM_and_fit(tc) +rng(1); +SPM = local_make_synthetic_spm(); +g = import_SPM(glm_map, SPM, 'noverbose'); + +tc.verifyEqual(g.level, 1); +tc.verifyTrue(g.is_timeseries); +tc.verifyEqual(g.TR, SPM.xY.RT); +tc.verifyEqual(size(g.X), size(SPM.xX.X)); +tc.verifyEqual(g.regressor_names, SPM.xX.name); +tc.verifyNumElements(g.onsets, 2); % onsets read-through +tc.verifyEqual([g.condition_names{:}], {'Cue', 'Pain'}); + +% Fit against a matching synthetic timeseries +nscan = size(SPM.xX.X, 1); +sim = fmri_data; +sim.dat = randn(25, nscan); +g = fit(g, sim, 'noverbose'); + +tc.verifyTrue(g.is_fitted); +tc.verifyEqual(g.dfe, nscan - size(SPM.xX.X, 2)); +end + + +% ===================================================================== +% Input-validation error paths +% ===================================================================== +function test_fit_without_design_errors(tc) +g = glm_map; % no X, no design +dat = fmri_data; dat.dat = randn(10, 5); +tc.verifyError(@() fit(g, dat, 'noverbose'), 'glm_map:NoDesign'); +end + + +function test_AR_requires_timeseries_errors(tc) +X = [ones(12, 1) randn(12, 1)]; +g = glm_map('X', X, 'level', 2); % is_timeseries == false +dat = fmri_data; dat.dat = randn(8, 12); +tc.verifyError(@() fit(g, dat, 'AR', 1, 'noverbose'), 'glm_map:ARnotTimeseries'); +end + + +function test_add_contrasts_size_mismatch_errors(tc) +X = [ones(15, 1) randn(15, 2)]; % 3 regressors +g = glm_map('X', X, 'level', 2); +tc.verifyError(@() add_contrasts(g, [1 0], {'bad'}), 'glm_map:ContrastSize'); +end + + +% ===================================================================== +% Local helpers +% ===================================================================== +function SPM = local_make_synthetic_spm() +% Minimal SPM12/SPM25-shared first-level structure for import tests. +nscan = 60; +TR = 1.5; +SPM = struct(); +SPM.xY.RT = TR; +SPM.nscan = nscan; +SPM.xBF = struct('name', 'hrf', 'length', 32, 'order', 1, 'dt', TR / 16, ... + 'T', 16, 'T0', 1, 'UNITS', 'secs', 'Volterra', 1); + +P = struct('name', 'none', 'P', [], 'h', 0); +U(1) = struct('dt', TR / 16, 'name', {{'Cue'}}, 'ons', [12 30 48]', 'dur', [0 0 0]', 'P', P, 'u', [], 'pst', []); +U(2) = struct('dt', TR / 16, 'name', {{'Pain'}}, 'ons', [18 36 54]', 'dur', [6 6 6]', 'P', P, 'u', [], 'pst', []); +SPM.Sess = struct('U', U, 'C', struct('C', [], 'name', {{}}), 'row', 1:nscan, 'col', 1:3, 'Fc', []); + +SPM.xX.X = [randn(nscan, 2) ones(nscan, 1)]; +SPM.xX.name = {'Sn(1) Cue*bf(1)', 'Sn(1) Pain*bf(1)', 'Sn(1) constant'}; +SPM.xX.iH = [1 2]; SPM.xX.iC = []; SPM.xX.iB = 3; SPM.xX.iG = []; +end From 12cd910763647d68cc31001e4b21fd5a8b4aeb61 Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Fri, 19 Jun 2026 14:30:46 -0400 Subject: [PATCH 7/8] Harden walkthrough CI tests: snapshots + fault-tolerant harness The nightly tests-walkthroughs job had never passed (42/42 red). It ran the CANlab_help_examples tutorials end-to-end via evalc() with no hardening, so the first orthviews/surface/prompt/missing-data error on the headless runner failed the whole test. Decouple the tests from the live tutorials and harden them: - Add verbatim snapshots of the 10 walkthroughs under walkthroughs/private/ (genpath-excluded, so they never shadow the real tutorials on CI). Refresh by overwriting from example_help_files/. - Add helpers/canlab_run_walkthrough_snapshot.m: runs a snapshot %%-cell by cell, headless, each cell in its own try/catch in a shared workspace, so a graphics-only section that fails on a headless runner does not abort the compute sections. - Add helpers/canlab_classify_environment_error.m: buckets caught errors into graphics / input / data / cascade / genuine (centralizes the heuristics previously inlined in canlab_test_help_examples). Environment buckets are skipped (Incomplete); only genuine errors fail, with an informative report naming the section, offending line, and error id. - Rewrite the 10 canlab_test_walkthrough_*.m wrappers to run their snapshot through the harness instead of evalc'ing the external script. - Add walkthroughs/README.md documenting the design, refresh process, and why these stay a separate nightly tier (~5 min local / ~15-20 min CI) rather than folding into the fast per-push suite. Full suite now: 10 passed, 0 failed, 0 incomplete. Two genuine tutorial bugs this surfaced (write-without-overwrite in walkthrough 3; dat_descrip -> metadata_table in 4b) are fixed in CANlab_help_examples and the snapshots refreshed accordingly. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../canlab_classify_environment_error.m | 147 +++++ .../helpers/canlab_run_walkthrough_snapshot.m | 207 +++++++ CanlabCore/Unit_tests/walkthroughs/README.md | 85 +++ ...nlab_test_walkthrough_1_installing_tools.m | 39 +- ...test_walkthrough_2_load_a_sample_dataset.m | 36 +- ...walkthrough_2b_basic_image_visualization.m | 35 +- ...lab_test_walkthrough_2c_loading_datasets.m | 35 +- ...nlab_test_walkthrough_3_voxelwise_t_test.m | 37 +- ...nlab_test_walkthrough_3b_atlases_and_ROI.m | 35 +- ...b_test_walkthrough_4_masking_and_writing.m | 37 +- ...lab_test_walkthrough_4b_3D_visualization.m | 35 +- .../canlab_test_walkthrough_5_regression.m | 35 +- ...st_walkthrough_7_multivariate_prediction.m | 35 +- .../private/canlab_help_1_installing_tools.m | 210 +++++++ .../canlab_help_2_load_a_sample_dataset.m | 199 +++++++ ...canlab_help_2b_basic_image_visualization.m | 149 +++++ .../private/canlab_help_2c_loading_datasets.m | 276 +++++++++ ...nlab_help_3_voxelwise_t_test_walkthrough.m | 224 +++++++ .../canlab_help_3b_atlases_and_ROI_analysis.m | 319 ++++++++++ ...b_help_4_masking_and_writing_nifti_files.m | 163 +++++ .../private/canlab_help_4b_3D_visualization.m | 468 +++++++++++++++ .../canlab_help_5_regression_walkthrough.m | 561 ++++++++++++++++++ ...ab_help_7_multivariate_prediction_basics.m | 391 ++++++++++++ 23 files changed, 3658 insertions(+), 100 deletions(-) create mode 100644 CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m create mode 100644 CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/README.md create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m create mode 100644 CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m diff --git a/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m new file mode 100644 index 00000000..a6e3e587 --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m @@ -0,0 +1,147 @@ +function category = canlab_classify_environment_error(ME, had_env_skip) +%CANLAB_CLASSIFY_ENVIRONMENT_ERROR Bucket a caught error for CI test handling. +% +% category = canlab_classify_environment_error(ME) +% category = canlab_classify_environment_error(ME, had_env_skip) +% +% Classifies an MException caught while running toolbox example/walkthrough +% code into one of a small set of buckets, so that test harnesses can decide +% whether to SKIP (mark Incomplete) or FAIL. The goal is to keep headless CI +% from going red over missing-display / missing-data / interactive-prompt +% conditions that are not real regressions, while still surfacing genuine +% code errors. +% +% This centralizes the heuristics previously hand-coded inside +% canlab_test_help_examples.m (skip_on_environment_error) so the walkthrough +% harness and the per-method help-example tests share one definition. +% +% :Inputs: +% +% **ME:** +% an MException (or any struct with .identifier, .message, and +% optionally .stack fields). +% +% **had_env_skip:** +% logical (default false). When true, "undefined variable/function" +% and "unrecognized field" errors are treated as CASCADE fallout from +% an earlier skipped section rather than genuine failures. Pass the +% running "have we already skipped something upstream" flag here. +% +% :Output: +% +% **category:** +% one of: +% 'graphics' - needs a display / OpenGL / Java / figure window, or +% the error originates inside a graphics-only code path +% (orthviews, surface rendering, etc.). +% 'input' - code tried to prompt for interactive user input, +% unavailable in batch/CI. +% 'data' - an optional data file (signature, atlas, feature set) +% is not on the CI path. +% 'cascade' - an undefined-variable/field error following an earlier +% environment skip (downstream fallout, not a new bug). +% 'genuine' - none of the above; treat as a real failure. +% +% :See also: canlab_run_walkthrough_snapshot, canlab_test_help_examples + +% .. +% Author: CANlab. Part of the CanlabCore Unit_tests headless-CI harness. +% .. + +if nargin < 2 || isempty(had_env_skip) + had_env_skip = false; +end + +msg = lower(ME.message); +id = ''; +if isprop(ME, 'identifier') || isfield(ME, 'identifier') + id = ME.identifier; +end + +% --------------------------------------------------------------------- +% Graphics: explicit graphics error ids, telltale message tokens, or a +% stack frame inside a known graphics-only function. spm_orthviews and the +% surface renderers fail or behave erratically on a headless runner. +% --------------------------------------------------------------------- +gfx_ids = {'MATLAB:graphics:opengl:Unavailable', ... + 'MATLAB:graphics:initialize', ... + 'MATLAB:class:InvalidHandle'}; +is_gfx = any(strcmp(id, gfx_ids)) || ... + strncmp(id, 'MATLAB:hg:', 10) || ... % handle-graphics errors + contains(msg, 'opengl') || contains(msg, 'display') || ... + contains(msg, 'java') || contains(msg, 'jvm') || ... + contains(msg, 'figure window') || contains(msg, 'no display') || ... + contains(msg, 'graphics') || ... + contains(msg, 'invalid or deleted object'); % set/get on a handle + % whose graphics section + % was skipped upstream + +gfx_fns = {'spm_orthviews', 'orthviews', 'spm_check_registration', ... + 'spm_figure', 'surface', 'render_on_surface', 'addbrain', ... + 'cluster_surf', 'isosurface', 'riverplot', 'tor_3d', ... + 'render_blobs', 'cluster_orthviews'}; +if ~is_gfx && isfield_or_prop(ME, 'stack') && ~isempty(ME.stack) + stack_names = {ME.stack.name}; + if any(ismember(stack_names, gfx_fns)) + is_gfx = true; + end +end +if is_gfx + category = 'graphics'; + return +end + +% --------------------------------------------------------------------- +% Interactive input: input() / keyboard prompts are unavailable in -batch. +% MATLAB reports MissingRequiredCapability with a "support for user input" +% message in that case. +% --------------------------------------------------------------------- +is_input = strcmp(id, 'MATLAB:services:MissingRequiredCapability') || ... + contains(msg, 'support for user input') || ... + contains(msg, 'input is not available'); +if is_input + category = 'input'; + return +end + +% --------------------------------------------------------------------- +% Missing optional data files (NPS+ signatures, Neurosynth feature set, +% Bianciardi atlas sources, etc.). load_image_set / annotate_* abort with a +% disp() + bare error('Exiting'), so the message is literally "Exiting". +% --------------------------------------------------------------------- +bare_exiting = isempty(id) && ... + ismember(strtrim(strrep(ME.message, '.', '')), {'Exiting', 'exiting'}); +is_data = contains(msg, 'cannot find images') || ... + contains(msg, 'find and add the file') || ... + contains(msg, 'not found in matlab path') || ... + contains(msg, 'no such file or directory') || ... + contains(msg, 'cannot find') || ... + bare_exiting; +if is_data + category = 'data'; + return +end + +% --------------------------------------------------------------------- +% Cascade: an earlier section was skipped, so a variable/field it would +% have defined is now missing. Treat as fallout, not a new regression. +% --------------------------------------------------------------------- +if had_env_skip && ( ... + strcmp(id, 'MATLAB:UndefinedFunction') || ... + strcmp(id, 'MATLAB:undefinedVarOrClass') || ... + strcmp(id, 'MATLAB:nonExistentField') || ... + contains(msg, 'undefined') || ... + contains(msg, 'unrecognized field')) + category = 'cascade'; + return +end + +category = 'genuine'; + +end + + +function tf = isfield_or_prop(obj, name) +% MException is an object (use isprop); a plain struct uses isfield. +tf = (isobject(obj) && isprop(obj, name)) || (isstruct(obj) && isfield(obj, name)); +end diff --git a/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m b/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m new file mode 100644 index 00000000..e8dcce50 --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m @@ -0,0 +1,207 @@ +function canlab_run_walkthrough_snapshot(tc, wt__snapshot_path) +%CANLAB_RUN_WALKTHROUGH_SNAPSHOT Run a frozen walkthrough snapshot under CI hardening. +% +% canlab_run_walkthrough_snapshot(tc, snapshot_path) +% +% Executes a copied CANlab_help_examples walkthrough script (a "snapshot" +% bundled under Unit_tests/walkthroughs/private/) one %%-cell at a time, in +% a single shared workspace, with headless-graphics and fault-tolerant error +% handling. Reports outcomes to the supplied matlab.unittest TestCase `tc`. +% +% Why snapshots instead of running the real tutorials. The walkthroughs in +% the CANlab_help_examples repo are teaching scripts, full of graphics +% (orthviews, surface), optional data sets, and the occasional interactive +% prompt. They should stay clean tutorials, not be retrofitted with test +% scaffolding, and CanlabCore CI should not depend on an external repo's +% un-hardened scripts. So we keep verbatim copies here and apply the +% hardening at run time, in one reusable place (this function) rather than +% editing each copy. Refresh a snapshot by overwriting the file under +% private/ from example_help_files/ — no re-hardening needed. +% +% Hardening applied: +% * Headless graphics: DefaultFigureVisible is forced 'off' for the run +% (callers also set this; doing it here makes the harness safe to call +% directly). Figures are closed between sections. +% * Section isolation: the script is split at its %% cell boundaries and +% each cell is run in its own try/catch, so a graphics-only section that +% fails on a headless runner does not abort the compute sections. +% * Error classification (see canlab_classify_environment_error): +% - graphics / input / data / cascade errors -> section SKIPPED, +% recorded, execution continues. +% - genuine errors -> the harness stops and the test FAILS with an +% informative message (which section, error id, message, and the +% offending top stack frame), because downstream cells usually +% cascade once a real error occurs. +% * Outcome mapping to matlab.unittest: +% - any genuine failure -> tc.verifyFail (test FAILS) +% - every section skipped (env) -> tc.assumeFail (test INCOMPLETE) +% - at least one section ran ok -> PASS, with a logged skip summary. +% +% Variable hygiene: every internal variable in this function is prefixed +% `wt__` because snapshot code is eval'd in this workspace and would +% otherwise clobber ordinary loop variables (i, n, t, r, ...). +% +% :Inputs: +% **tc:** a matlab.unittest.TestCase (from functiontests). +% **snapshot_path:** absolute path to the snapshot .m file to run. +% +% :See also: canlab_classify_environment_error, canlab_run_all_tests + +% .. +% Author: CANlab. Part of the CanlabCore Unit_tests headless-CI harness. +% .. + +tc.assertTrue(exist(wt__snapshot_path, 'file') == 2, ... + sprintf('Walkthrough snapshot not found: %s', wt__snapshot_path)); + +[~, wt__name] = fileparts(wt__snapshot_path); +wt__src = fileread(wt__snapshot_path); +wt__sections = local_split_sections(wt__src); +wt__nsec = numel(wt__sections); + +% Headless: callers set this too, but be self-sufficient. +wt__prev_vis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); +wt__cleanup = onCleanup(@() set(0, 'DefaultFigureVisible', wt__prev_vis)); + +wt__skipped = {}; % cellstr of human-readable skip notes +wt__ran_real = 0; % count of sections that executed without error +wt__had_env_skip = false; % drives cascade classification downstream +wt__failure = []; % first genuine failure, if any (struct) + +for wt__s = 1:wt__nsec + + wt__code = wt__sections{wt__s}; + + % Skip cells that are pure comment/whitespace, or that define a + % function (snapshots are scripts; a trailing function would not eval). + if local_is_noop(wt__code) + continue + end + + try + evalc(wt__code); % runs in THIS workspace; vars persist + wt__ran_real = wt__ran_real + 1; + catch wt__ME + wt__cat = canlab_classify_environment_error(wt__ME, wt__had_env_skip); + if strcmp(wt__cat, 'genuine') + wt__failure = struct( ... + 'section', wt__s, ... + 'id', wt__ME.identifier, ... + 'message', wt__ME.message, ... + 'where', local_top_frame(wt__ME), ... + 'snippet', local_first_stmt(wt__code)); + break % stop; later cells will cascade + else + wt__had_env_skip = true; + wt__skipped{end+1, 1} = sprintf(' section %d/%d [%s]: %s', ... + wt__s, wt__nsec, wt__cat, local_oneline(wt__ME.message)); %#ok + end + end + + close all force +end + +% --------------------------------------------------------------------- +% Map collected outcomes onto the TestCase. +% --------------------------------------------------------------------- +if ~isempty(wt__failure) + wt__report = sprintf([ ... + '%s: genuine error at section %d of %d.\n' ... + ' identifier: %s\n' ... + ' message: %s\n' ... + ' section starts: %s\n' ... + ' at: %s\n' ... + ' (%d earlier section(s) skipped for environment reasons)'], ... + wt__name, wt__failure.section, wt__nsec, ... + wt__failure.id, local_oneline(wt__failure.message), ... + wt__failure.snippet, wt__failure.where, numel(wt__skipped)); + tc.verifyFail(wt__report); + +elseif wt__ran_real == 0 + tc.assumeFail(sprintf( ... + '%s: all %d section(s) skipped for environment reasons:\n%s', ... + wt__name, wt__nsec, strjoin(wt__skipped, newline))); + +else + if ~isempty(wt__skipped) + tc.log(1, sprintf('%s: %d of %d section(s) skipped (environment):\n%s', ... + wt__name, numel(wt__skipped), wt__nsec, strjoin(wt__skipped, newline))); + end + tc.verifyTrue(true, sprintf('%s ran %d section(s) without genuine error', ... + wt__name, wt__ran_real)); +end + +end + + +% ===================================================================== +% Local helpers +% ===================================================================== + +function sections = local_split_sections(src) +% Split source at %% cell markers. Each section keeps its leading marker +% line (a comment, harmless to eval). Everything before the first marker is +% section 1. +lines = regexp(src, '\r\n|\r|\n', 'split'); +is_head = ~cellfun('isempty', regexp(lines, '^\s*%%', 'once')); +starts = find(is_head); +if isempty(starts) || starts(1) ~= 1 + starts = [1, starts]; +end +starts = unique(starts); +ends = [starts(2:end) - 1, numel(lines)]; +sections = cell(numel(starts), 1); +for k = 1:numel(starts) + sections{k} = strjoin(lines(starts(k):ends(k)), newline); +end +end + + +function tf = local_is_noop(code) +% True if the cell has no executable content: blank, comment-only, or a +% function definition (which cannot be eval'd as a statement). +stripped = regexprep(code, '%[^\n]*', ''); % drop line comments +stripped = strtrim(stripped); +tf = isempty(stripped) || ~isempty(regexp(stripped, '^function\b', 'once')); +end + + +function s = local_oneline(msg) +% Collapse a multi-line error message to a single trimmed line. +s = strtrim(regexprep(msg, '\s+', ' ')); +end + + +function s = local_top_frame(ME) +% "funcname (line N)" for the deepest stack frame that is not this harness +% (errors raised at a section's top level otherwise just point back at the +% evalc call site here, which is unhelpful). Falls back to a marker. +if isempty(ME.stack) + s = ''; + return +end +this_fn = mfilename; +for k = 1:numel(ME.stack) + if ~strcmp(ME.stack(k).name, this_fn) + s = sprintf('%s (line %d)', ME.stack(k).name, ME.stack(k).line); + return + end +end +s = ''; +end + + +function s = local_first_stmt(code) +% First non-comment, non-blank line of a section, trimmed, for the report. +lines = regexp(code, '\r\n|\r|\n', 'split'); +s = ''; +for k = 1:numel(lines) + ln = strtrim(lines{k}); + if ~isempty(ln) && ~strncmp(ln, '%', 1) + s = ln; + if numel(s) > 100, s = [s(1:97) '...']; end + return + end +end +end diff --git a/CanlabCore/Unit_tests/walkthroughs/README.md b/CanlabCore/Unit_tests/walkthroughs/README.md new file mode 100644 index 00000000..4bfbbd52 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/README.md @@ -0,0 +1,85 @@ +# Walkthrough integration tests (Tier B) + +These tests smoke-test the end-to-end CANlab tutorials ("walkthroughs") from +the [CANlab_help_examples](https://github.com/canlab/CANlab_help_examples) +repository, to catch the case where a tutorial silently bit-rots against the +current toolbox. + +They are **slow** and depend on more sibling repos and optional data than the +fast per-push unit suite, so they are a separate tier, run on the nightly +`tests-walkthroughs` GitHub Actions workflow, **not** on every push. They are +skipped by `canlab_run_all_tests` unless you ask for them: + +```matlab +canlab_run_all_tests('Walkthroughs', 'only') % just the walkthroughs +canlab_run_all_tests('Walkthroughs', 'include') % unit suite + walkthroughs +canlab_run_all_tests % default: walkthroughs skipped +``` + +## Design: snapshots + a hardening harness + +The tutorials in CANlab_help_examples are teaching scripts. They should stay +clean tutorials — full of `orthviews`, `surface`, optional signature/atlas +data, and the occasional interactive prompt — and **not** be retrofitted with +test scaffolding. CanlabCore CI should also not depend on an external repo's +un-hardened scripts at run time. So instead of executing the tutorials in +place, this folder works from frozen copies: + +``` +walkthroughs/ + canlab_test_walkthrough_*.m <- thin matlab.unittest wrappers (discovered) + private/ + canlab_help_*.m <- VERBATIM snapshots of the tutorials + README.md +``` + +- `private/` is excluded from `genpath`, so the snapshots never go on the + MATLAB path and never collide with the real tutorials when both repos are + checked out on CI. The harness reads them by absolute path. +- Each `canlab_test_walkthrough_*.m` wrapper just points at its snapshot and + calls the shared harness. The hardening lives in one place + (`Unit_tests/helpers/canlab_run_walkthrough_snapshot.m`), not duplicated into + each copy — so refreshing a snapshot is a clean file overwrite. + +### What the harness does + +`canlab_run_walkthrough_snapshot(tc, snapshot_path)`: + +1. Forces `DefaultFigureVisible 'off'` (headless). +2. Splits the snapshot at its `%%` cell boundaries and runs each cell in its + own `try/catch`, in one shared workspace, so a graphics-only section that + fails on a headless runner does not abort the compute sections. +3. Classifies any caught error via + `Unit_tests/helpers/canlab_classify_environment_error.m`: + - **graphics / input / data / cascade** → section **skipped** (recorded), + execution continues. + - **genuine** → the harness stops and the test **fails** with an informative + report (which section, the offending line, error id and message). +4. Maps outcomes to `matlab.unittest`: + - any genuine failure → **Failed** + - every section skipped for environment reasons → **Incomplete** + - at least one section ran → **Passed** (with a logged skip summary) + +So the nightly goes red only on a *genuine* regression, and tells you exactly +where; missing-display / missing-data / prompt conditions become skips. + +## Refreshing a snapshot + +When a tutorial changes upstream, re-copy it verbatim: + +```bash +cp ../../../../CANlab_help_examples/example_help_files/canlab_help_5_regression_walkthrough.m \ + private/canlab_help_5_regression_walkthrough.m +``` + +No re-hardening needed — the harness handles it. Keep the copies verbatim so +they can be diffed against upstream. + +## Why a separate nightly tier (not folded into the per-push suite) + +Measured wall-time for the full set is ~5 min on a fast workstation and an +estimated ~15–20 min on the 2-core GitHub Actions Linux runner. The per-push +`tests` suite is the fast gate that every PR blocks on; adding 15–20 min plus +graphics/data-dependent flakiness to it would slow iteration and make the +required check unreliable. End-to-end tutorial runs are exactly the kind of +slow, environment-sensitive integration test that belongs in a nightly tier. diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m index 5398e0e4..01db58ce 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m @@ -1,34 +1,45 @@ function tests = canlab_test_walkthrough_1_installing_tools -%CANLAB_TEST_WALKTHROUGH_1_INSTALLING_TOOLS End-to-end smoke test of canlab_help_1_installing_tools. +%CANLAB_TEST_WALKTHROUGH_1_INSTALLING_TOOLS Hardened smoke test: Installing tools walkthrough. % -% Wrapper around CANlab_help_examples/example_help_files/canlab_help_1_installing_tools.m. -% Acts as a Tier B integration test: confirms the walkthrough runs end-to-end -% on the current toolbox without erroring. Does not check pixel content, -% statistics, or output values; reaching the end of the script counts as a pass. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_1_installing_tools" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. % -% Skipped (Incomplete) if CANlab_help_examples is not on the MATLAB path. +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> -% Run each walkthrough in its own temp working directory so output files -% (NIfTI writes, saved figures) don't collide with prior runs. +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_1_installing_tools.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_1_installing_tools'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m index 0ea0f7cd..ba5d6a47 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m @@ -1,29 +1,45 @@ function tests = canlab_test_walkthrough_2_load_a_sample_dataset -%CANLAB_TEST_WALKTHROUGH_2_LOAD_A_SAMPLE_DATASET Smoke test of canlab_help_2_load_a_sample_dataset. +%CANLAB_TEST_WALKTHROUGH_2_LOAD_A_SAMPLE_DATASET Hardened smoke test: Load a sample dataset walkthrough. % -% Tier B integration test: runs the walkthrough end-to-end and counts a -% return without error as a pass. See canlab_test_walkthrough_1_installing_tools -% for shared rationale. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2_load_a_sample_dataset" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2_load_a_sample_dataset.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2_load_a_sample_dataset'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m index a836e12a..09d78caa 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_2b_basic_image_visualization -%CANLAB_TEST_WALKTHROUGH_2B_BASIC_IMAGE_VISUALIZATION Smoke test of canlab_help_2b_basic_image_visualization. +%CANLAB_TEST_WALKTHROUGH_2B_BASIC_IMAGE_VISUALIZATION Hardened smoke test: Basic image visualization walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2b_basic_image_visualization" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2b_basic_image_visualization.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2b_basic_image_visualization'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m index f1d3a208..17258aae 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_2c_loading_datasets -%CANLAB_TEST_WALKTHROUGH_2C_LOADING_DATASETS Smoke test of canlab_help_2c_loading_datasets. +%CANLAB_TEST_WALKTHROUGH_2C_LOADING_DATASETS Hardened smoke test: Loading datasets walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2c_loading_datasets" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2c_loading_datasets.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2c_loading_datasets'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m index dc36fced..3e0f3b72 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m @@ -1,30 +1,45 @@ function tests = canlab_test_walkthrough_3_voxelwise_t_test -%CANLAB_TEST_WALKTHROUGH_3_VOXELWISE_T_TEST Smoke test of canlab_help_3_voxelwise_t_test_walkthrough. +%CANLAB_TEST_WALKTHROUGH_3_VOXELWISE_T_TEST Hardened smoke test: Voxelwise t-test walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. The WorkingFolderFixture in setup -% provides a clean cwd, which matters here because the walkthrough writes -% NIfTI output files. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_3_voxelwise_t_test_walkthrough" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_3_voxelwise_t_test_walkthrough.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_3_voxelwise_t_test_walkthrough'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m index bbf8135c..8aa9a9be 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_3b_atlases_and_ROI -%CANLAB_TEST_WALKTHROUGH_3B_ATLASES_AND_ROI Smoke test of canlab_help_3b_atlases_and_ROI_analysis. +%CANLAB_TEST_WALKTHROUGH_3B_ATLASES_AND_ROI Hardened smoke test: Atlases and ROI analysis walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_3b_atlases_and_ROI_analysis" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_3b_atlases_and_ROI_analysis.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_3b_atlases_and_ROI_analysis'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m index 7958e2eb..b3c54ebb 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m @@ -1,30 +1,45 @@ function tests = canlab_test_walkthrough_4_masking_and_writing -%CANLAB_TEST_WALKTHROUGH_4_MASKING_AND_WRITING Smoke test of canlab_help_4_masking_and_writing_nifti_files. +%CANLAB_TEST_WALKTHROUGH_4_MASKING_AND_WRITING Hardened smoke test: Masking and writing NIfTI files walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. The WorkingFolderFixture in setup -% provides a clean cwd, which matters here because the walkthrough writes -% NIfTI output files. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_4_masking_and_writing_nifti_files" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_4_masking_and_writing_nifti_files.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_4_masking_and_writing_nifti_files'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m index cc52a045..83f8d922 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_4b_3D_visualization -%CANLAB_TEST_WALKTHROUGH_4B_3D_VISUALIZATION Smoke test of canlab_help_4b_3D_visualization. +%CANLAB_TEST_WALKTHROUGH_4B_3D_VISUALIZATION Hardened smoke test: 3D visualization walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_4b_3D_visualization" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_4b_3D_visualization.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_4b_3D_visualization'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m index b0d2ca7b..a4dcfc25 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_5_regression -%CANLAB_TEST_WALKTHROUGH_5_REGRESSION Smoke test of canlab_help_5_regression_walkthrough. +%CANLAB_TEST_WALKTHROUGH_5_REGRESSION Hardened smoke test: Regression walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_5_regression_walkthrough" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_5_regression_walkthrough.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_5_regression_walkthrough'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m index f8d86804..64e9cdb3 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_7_multivariate_prediction -%CANLAB_TEST_WALKTHROUGH_7_MULTIVARIATE_PREDICTION Smoke test of canlab_help_7_multivariate_prediction_basics. +%CANLAB_TEST_WALKTHROUGH_7_MULTIVARIATE_PREDICTION Hardened smoke test: Multivariate prediction basics walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_7_multivariate_prediction_basics" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_7_multivariate_prediction_basics.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_7_multivariate_prediction_basics'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m new file mode 100644 index 00000000..2a55efe9 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m @@ -0,0 +1,210 @@ +%% Installing CANlab Core Tools + +%% About the CANlab Core Tools repository and CANlab toolboxes +% +% The CANlab Core Tools repository is on https://github.com/canlab/CanlabCore +% It contains core tools for MRI/fMRI/PET analysis from the Cognitive and +% Affective Neuorscience Lab (Tor Wager, PI) and our collaborators. Many of +% these functions are needed to run other toolboxes, e.g., the CAN lab?s +% multilevel mediation and Martin Lindquist?s hemodynamic response estimation +% toolboxes. +% +% The tools include object-oriented tools for doing neuroimaging analysis with +% simple commands and scripts that provide high-level functionality for +% neuroimaging analysis. For example, there is an "fmri_data" object type +% that contains neuroimaging datasets (both PET and fMRI data are ok, despite +% the name). If you have created and object called my_fmri_data_obj, then +% plot(my_fmri_data_obj) will generate a series of plots specific to neuroimaging +% data, including an interactive brain viewer (courtesy of SPM software). +% predict(my_fmri_data_obj) will perform cross-validated multivariate prediction +% of outcomes based on brain data. ica(my_fmri_data_obj) will perform independent +% components analysis on the data, and so forth. +% +% There are a wide range of functions and techniques that we've used in the +% CANlab. Most of these are supported by specific functionality in our +% object-oriented tools. You can download other toolboxes from: +% https://github.com/canlab +% +% And you can find an overview of core object-oriented tool functionality +% and demos/walkthroughs on https://canlab.github.io +% +% <> +% +% + +%% CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> + +%% Dependencies +% +% * Matlab statistics toolbox +% * Matlab signal processing toolbox +% * Statistical Parametric Mapping (SPM) software https://www.fil.ion.ucl.ac.uk/spm/ +% (this is used for image reading/writing and the orthviews function.) +% +% For full functionality, the other toolboxes below are recommended: +% + +%% Installing SPM +% +% *!* Important: You will also need to install spm12 separately, and add it +% to your matlab path. +% The latest version at the time of writing is SPM12, which can be downloaded here: +% +% +% + +%% Quick start instructions +% +% # Download SPM12 from https://www.fil.ion.ucl.ac.uk/spm/ +% # Find the SPM12 folder and drag it into your command window. +% Matlab will go to that directory (execute the cd command). +% # In Matlab, type |addpath pwd|. This adds the main SPM12 folder to your path. +% # Type |addpath canonical|. This adds the |canonical| subfolder to your +% path. +% # Go to https://github.com/canlab/ and click on CanlabCore +% # Sign into Github (top right corner) +% # Click "Clone or Download" +% # Click "Open in Desktop" (if you have Github Desktop) to clone the +% repository +% # Find the folder and drag it into your Matlab command window. +% # In Matlab, type |g = genpath(pwd)|. This lists all the subfolders in a string |g| +% # Type |addpath(g)|. This adds them to your path. +% # Type |savepath|. This saves your path for future use. +% +% Repeat the Clone/genpath/addpath steps for the Neuroimaging_Pattern_Masks +% repository, and save the path again. + + +%% +% Quick start: Batch install script +% +% The main goal of the code below is to run: +% |canlab_toolbox_setup| +% This is a script that attempts to: +% # Check for CANlab repositories on your computer +% # Download any that are missing +% # Add them to your Matlab path with subfolders. +% +% _Note_: Unfortunately as of Jan 2020 this does not work well on +% everyone's computers, due to OS-related variability. The workaround is to +% install them manually as per the instructions above. +% +% |canlab_toolbox_setup| is in the CanlabCore repository. So before you run it +% the CANlab_Core_Tools must be added to your path with +% subfolders. Otherwise, you will get errors. +% +% The script canlab_toolbox_setup can help download and install the +% toolboxes you need. +% +% First, go to a folder where you want to install +% toolboxes/repositories on your hard drive. Mine are in "Github" and I've +% already installed CanlabCore, so let's find it and go there: + +% Locate the CanlabCore Github files on your local computer: +% ------------------------------------------------------------------------ + +mypath = what(fullfile('CanlabCore', 'CanlabCore')); + +% Check if we've got it +if isempty(mypath) + disp('Download CanlabCore from Github, and go to that folder in Matlab') + disp('by dragging and dropping it from Finder or Explorer into the Matlab Command Window') + return +else + % Add CanlabCore to Matlab path with subfolders + g = genpath(mypath(1).path); + addpath(g); + +end + +mypath = mypath(1).path; + +% Set the folder in which the repositories will be installed +% ------------------------------------------------------------------------ + +% This is the location (folder) in which the repositories will be installed: +base_dir_for_repositories = fileparts(fileparts(mypath)); + +cd(base_dir_for_repositories) +fprintf('\nInstalling repositories in %s\n', base_dir_for_repositories); + +% Find the setup file, canlab_toolbox_setup.m +% ------------------------------------------------------------------------ + +% This is the file +setupfile = fullfile(mypath, 'canlab_toolbox_setup.m'); + +% Make sure we've got it +if ~exist(setupfile, 'file') + disp('Something went wrong. I can''t find canlab_toolbox_setup.m'); + disp('This file should be included in the CanlabCore Github repository.') +end + +% Run the setup script: +% ------------------------------------------------------------------------ +% +% Note: This can be system-dependent and may not work on all computers. +% If it doesn't, you can default to cloning the CANLab repositories and +% adding them with subfolders to your Matlab path. The most important ones +% starting out are CanlabCore and Neuroimaging_Pattern_Masks +% e.g., +% |git clone https://github.com/canlab/CanlabCore.git| +% |git clone https://github.com/canlab/Neuroimaging_Pattern_Masks.git| + +canlab_toolbox_setup + +% This will attempt to locate toolboxes, add them to your path, and give +% you the option to download them from Github if it can't find them. + +% It looks for and installs toolboxes in the current directory, +% which (thanks to the code above) is the path name in base_dir_for_repositories + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m new file mode 100644 index 00000000..8c8b169f --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m @@ -0,0 +1,199 @@ +%% Load and explore a sample dataset +% This example shows how to load a pre-set example fMRI dataset into +% an fmri_data object. + +%% General instructions +% +% Before you start, the CANlab_Core_Tools must be added to your path with +% subfolders. Otherwise, you will get errors. +% +% The script canlab_toolbox_setup can help download and install the +% toolboxes you need. +% +% Sample datasets are in the "Sample_datasets" folder in CANlab_Core_Tools. +% Many tutorials apply pre-trained patterns and masks. +% These are stored in this Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% This example will use emotion regulation data in the folder: +% "Wager_et_al_2008_Neuron_EmotionReg" +% The dataset is a series of contrast images from N = 30 participants. +% Each image is a contrast image for [reappraise neg vs. look neg] +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% Here are a couple of helpful functions we will use for display: +% (you can ignore these.) +dashes = '----------------------------------------------'; +printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes); + +%% The fmri_data object +% +% The fmri_data object class is one of the most important, basic types of +% objects in the CANlab object-oriented toolbox. It stores image data in a +% matrix form, which is more space efficient and friendly for analysis with +% various software packages/algorithms. +% +% For philosophy, see: +% +% +% For more info on the object-oriented approach, see: +% +% +% When you call the _class constructor_ fmri_data(), this is what it does: +% +% <> + +%% Section 1: The quick and easy way to load a pre-specified dataset +% +% The function load_image_set has a number of pre-defined image sets that +% you can load with one simple command. This is the easiest way to load +% sample data. The images must be on your Matlab path for this to work. + +[data_obj, subject_names, image_names] = load_image_set('emotionreg'); + +% data_obj is an fmri_data object containing all 30 images. +% subject_names is a list of short names for each image. +% image_names is a list of the full image names with their path. + +%% Section 2: Manual load. Use filenames to find file names +% +% First, we need to list the file names in a string matrix. +% Then, we can load them into an fmri_data object +% We will use the filenames function to get the names, and the fmri_data +% object constructor function to load the data. + +% First, check whether images are on your path: +% We will search for one image and save the path name. + +printhdr('Check that we can find data images:') +myfile = which('con_00810001.img'); +mydir = fileparts(myfile); +if isempty(mydir), disp('Uh-oh! I can''t find the data.'), else disp('Data found.'), end + +% Now we can list all the file names. + +printhdr('Find files and get their names:') +image_names = filenames(fullfile(mydir, 'con_008100*img'), 'absolute'); +disp('Done.'); + +% Now load them into an fmri_data object. + +printhdr('Loading the image data into the object:') +data_obj = fmri_data(image_names); + +% This is the gateway to doing many other things, which are explained in +% other help files. But just to get us started, let's run through a few +% basic things we can do. We'll mainly just look at some standard plots of +% the data. + +%% Section 3: Plot the data we just loaded +% +% Operations that we can perform on fmri_data objects are called methods. +% You can see a list of methods by typing methods(data_obj). +% Here, we'll call the plot method to visualize the data. + +plot(data_obj) +snapnow + +%% Section 4: Get help on an object +% +% Objects have associated help files. +% Type: + +help fmri_data + +% ...to get help for the fmri_data() class constructor and object. +% +% Some objects also have doc files that contain additional information. +% Try: + +doc fmri_data + +%% Section 5: Explore methods +% +% Objects have methods associated with them, or things you can do with +% them. To see a list of methods for the fmri_data object, type: + +methods(fmri_data) + +%% +% You can also use the name of an _instance_ of an object, i.e., a variable +% in the workspace. |methods(data_obj)| produces the same list. +% In addition, you can get help on any of the object's methods by typing: +% |object class name_dot_method name|. e.g., + +help fmri_data.mean + +% Note: There may be trouble accessing these help files in Matlab R2019b. +% Works in 2018 and earlier. But maybe this is not a general problem...give +% it a try. + +%% Section 6: Explore on your own +% +% 1. Try to run a couple of other methods on your fmri_data object. Some +% require additional inputs, but many do not. Look at the help for more +% info. What do you see? Can you return any meaningful output, and if so, +% what did the method do? +% + +% That's it for this section!! + +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m new file mode 100644 index 00000000..b257361c --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m @@ -0,0 +1,149 @@ +%% Basic Image Visualization + +%% +% Using Basic Plot Methods (plot, orthviews, and montage) to Examine a Dataset +% +% Objects can have methods with intuitive names, some of which overlap +% with names of functions in Matlab or other toolboxes. fmri_data.plot() +% is one of these. When you call plot() and pass in an fmri_data object, +% you invoke the fmri_data object method, and a special plot for fmri_data +% objects is produced. +% +% You can list methods for an object class (e.g., fmri_data) by typing: +% +% |methods(fmri_data)| +% +% You can get help for a method by typing +% |help .| +% e.g., |help fmri_data.plot| +% +% The plot() method takes an fMRI data object as a parameter and produces +% an SPM Orthview presentation along with 6 plots of the data. +% +% The plot method expects an fMRI data object to be passed in. We can +% create an fMRI data object using the emotion regulation dataset +% via the following code: + +[image_obj, networknames, imagenames] = load_image_set('emotionreg', 'noverbose'); + +%% +% Once created, we can pass this data object to the plot function to get +% the entire set of outputs, including Matlab console output regarding +% outliers and corresponding data visualizations, using the simple command: + +plot(image_obj); + +%% Explaining the output +% +% The help file for fmri_data.plot did object has more information about the plots: +% +% e.g., |help fmri_data.plot| +% +% e.g., |help image_obj.plot| + +help fmri_data.plot + +%% Use Orthviews to visualize the mean image +% +% Some other commonly used methods to display images are the orthviews() +% and montage() methods. fmri_data objects have these methods, and other +% object classes do to, like the statistic_image, atlas, and region +% classes. +% +% First, let's create a mean image for the dataset: + +m = mean(image_obj); +orthviews(m); +drawnow, snapnow + +%% Threshold and display +% let's threshold it using the threshold() method. +% Here we'll exclude values between -1 and 1 and view the extreme values: + +m = threshold(m, [-1 1], 'raw-outside'); +orthviews(m); +drawnow, snapnow + +%% Create and display a region object +% We can create a region class object, another type of object, from the +% thresholded image. This has additional info and options about each +% contiguous 'blob' in the suprathreshold map: + +r = region(m); +orthviews(r); +drawnow, snapnow + +%% Orthviews options +% Orthviews methods have a range of options. +% They use the cluster_orthviews function, which uses spm_orthviews +% See |help cluster_orthviews| for options. +% +% Let's try one that visualizes each contiguous blob in a different, solid +% color: + +orthviews(r, 'unique', 'solid'); +drawnow, snapnow + +%% Use montage to visualize the thresholded mean image +% +% Sometimes, we want to view map that shows a canonical range of slices. +% This is really useful for producing standard output for papers +% Arguably, one should *always* view and publish montage maps showing all +% slices, so as to show the "whole picture" and not omit any results. +% +% You can customize this a lot, as it uses the fmridisplay() object class, +% which allows you to add custom montages (in axial, saggital, and coronal +% orientations, add blobs of various types, and remove them and re-plot. +% See |help fmridisplay| and |help fmridisplay.addblobs| for more details. +% +% For now, we'll just stick to a basic plot. We'll first create an empty +% figure,then plot the montage on it. + +create_figure('montage'); axis off; +montage(m); +drawnow, snapnow + +% We've already thresholded it, so it'll use the previous threshold. +% however, we can re-threshold the image and redisplay it as well. + +%% Use montage to visualize each blob in a thresholded map +% +% A really useful thing to do is to take a region object, often from a +% thresholded map, and visualize each region. +% the montage() methods also have a number of options. +% Let's try one, |'regioncenters'|, that plots each blob (region) on a +% separate slice. +% +% Furthermore, we can use the |'colormap'| option to view the regions with +% colors mapped to their associated values, e.g., hot colors for positive +% values and cool colors for negative values. +% +% Finally, we might want to assign names to each region based on an atlas. +% We'll do that before plotting, so that the names appear on the plots. +% These names are saved in the r(i).shorttitle field, for each region i. +% The region.table() method automatically labels them as well. +% You can customize the atlas used; the default is the 'canlab2018_2mm' +% atlas (see |help load_atlas| for more info.) + +r = autolabel_regions_using_atlas(r); + +montage(r, 'regioncenters', 'colormap'); +drawnow, snapnow + +% Note that some large regions may span multiple areas. +% This can happen if the various regions are connected by contiguous +% suprathreshold voxels. + +%% Explore on your own +% +% 1. Try to re-threshold the image using some values you choose and re-plot. +% Look at the help for more options on thresholding. and pick one. What do you see? +% Try to plot only voxels with positive values in the mean image. +% +% 2. Try to bring up only *one* region from the region object (r) in orthviews. +% Can you visualize it in all three views? +% +% 3. Try using a couple of other display options in the montage() and orthviews() +% methods. What do you see? + +% That's it for this section!! \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m new file mode 100644 index 00000000..11aed9bd --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m @@ -0,0 +1,276 @@ +%% Datasets used in CANlab tutorials +% +% Note: this report was generated from |canlab_help_2c_loading_datasets.m| +% in the repository +% + +%% The Neuroimaging_pattern_masks Github repository and website +% +% For these walkthroughs, we'll load several image datasets. +% Some are stored in Github, if the files are small enough. +% The Wager 2008 emotion regulation dataset, for example, is +% in the CANlab Core repository. Others are on figshare or Neurovault. +% +% Many tutorials apply pre-trained patterns and masks. +% These are stored in this Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% To run this tutorial and many others in this series, you'll need to +% download the Neuroimaging_Pattern_Masks Github repository and add it to +% your Matlab path with subfolders. +% The Github site has three main types of datasets, shown here: +% +% <> +% +%% Registries for easy loading +% There are several functions that contain sets of images that you can load +% by name. For example: +% +% |load_image_set()| includes a registry of datasets that you can access by +% name. Try |help load_image|set| for a list of images you can load by +% keyword. +% +% |load_atlas()| also has a named set of atlases/brain parcellations you can load. +% +% |canlab_load_ROI| has a registry of many named regions derived from previous +% studies. This is particularly useful for loading a subcortical ROI and +% visualizing it or applying it to data. + +%% Emotion regulation dataset +% "Wager_et_al_2008_Neuron_EmotionReg" +% The dataset is a series of contrast images from N = 30 participants. +% Each image is a contrast image for [reappraise neg vs. look neg] +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% To load, try the following: + +[data_obj, subject_names, image_names] = load_image_set('emotionreg', 'noverbose'); + +montage(mean(data_obj), 'trans', 'compact2'); +drawnow, snapnow + +% data_obj is an fmri_data object containing all 30 images. +% subject_names is a list of short names for each image. +% image_names is a list of the full image names with their path. + +%% Buckner lab resting-state connectivity maps +% +% A popular series of masks is a set of "networks" developed from resting-state fMRI +% connectivity in 1,000 people. For our purposes, the "network" maps +% consist of a set of voxels that load most highly on each of 7 ICA +% components from the 1,000 Functional Connectomes project. +% These were published by Randy Buckner's lab in 2011 in three papers: +% Choi et al., Buckner et al., and Yeo et al. We'll use the cortical +% networks from Yeo et al. here. +% + +[bucknermaps, mapnames] = load_image_set('bucknerlab', 'noverbose'); % loads 7 masks from Yeo et al. +disp(char(mapnames)) + +% Create a montage plot +o2 = montage(get_wh_image(bucknermaps, 1), 'trans', 'compact2', 'color', rand(1, 3)); +for i = 2:7, o2 = addblobs(o2, region(get_wh_image(bucknermaps, i)), 'trans', 'color', rand(1, 3)); end +drawnow, snapnow + +%% Pain, Cognition, Emotion balanced N = 270 dataset +% +% This dataset is an fmri_data object class object created using CANlab tools (canlab.github.io). +% It contains 270 single-participant fMRI contrast maps across 18 studies (with 15 participants per study). +% Studies are grouped into three domains: Pain, Cognitive Control, and Negative Emotion, with 9 studies each. +% Each domain includes 3 subdomains, with 3 studies each. +% +% The dataset is from Kragel et al. 2018, Nature Neuroscience. +% It is too large for Github, and it's stored on Neurovault.org +% as collection #504. You could get it using CANlab tools like this: +% +% |[files_on_disk, url_on_neurovault, mycollection, myimages] = retrieve_neurovault_collection(504);| +% +% It is also available on Figshare, with DOI 10.6084/m9.figshare.24033402 +% https://figshare.com/ndownloader/files/42143352 +% +% If you download from Neurovault, you'd +% have to add metadata for the study category labels to use it. +% Therefore, we suggest you use the CANLab load_image_set() function, as below. +% It saves a file on your hard drive the first time you run it: +% +% |kragel_2018_nat_neurosci_270_subjects_test_images.mat| +% +% This file also has metadata that is not necessarily on Neurovault. +% If you save this file somewhere on your Matlab path, you'll be able to +% reload and reuse the dataset easily. + +[test_images, names] = load_image_set('kragel18_alldata', 'noverbose'); + +% This field contains a table object with metadata for each image: +metadata = test_images.metadata_table; + +disp('Metadata available in test_images.metadata_table:') +metadata(1:5, :) + +% Make a plot of the spatial correlation of the average image for each study + +imgs = cellstr(test_images.image_names); +m = mean(test_images); + +for i = 1:length(names) + + % Create a mean image for each study and store in "m" object. + + wh = metadata.Studynumber == i; + studymean = mean(get_wh_image(test_images, find(wh))); + m.dat(:, i) = studymean.dat; + +end + +disp('Map of spatial correlations across the mean images for each study'); + +plot_correlation_matrix(m.dat, 'names', names, 'partitions', [ones(1, 6) 2*ones(1, 6) 3*ones(1, 6)], 'partitionlabels', {'Pain', 'Cognition', 'Emotion'}); +drawnow, snapnow + +%% Pain across six intensity levels per person (BMRK3) +% +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). Each image is the average over several +% (4-8) trials of heat delivered at a single stimulus intensity, ranging +% from 44.3 - 49.3 degrees C in one-degree increments. Each image is also +% paired with an average reported pain value for that set of trials, rated +% immmediately after heat experience. +% +% This dataset is interesting for mixed-effects and predictive analyses, as +% it has both within-person and between-person sources of variance. +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj, an fmri_data class object. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + +subject_id = image_obj.additional_info.subject_id; +ratings = reshape(image_obj.Y, 6, 33)'; +temperatures = image_obj.additional_info.temperatures; + +% Plot the ratings + +create_figure('ratings'); +title('Pain ratings by stimulus intensity in BMRK3 dataset') +hold on; plot(ratings', '-', 'Color', [.7 .7 .7], 'LineWidth', .5); +lineplot_columns(ratings, 'color', [.7 .3 .3], 'markerfacecolor', [1 .5 0]); +xlabel('Temperature'); +ylabel('Rating'); +set(gca, 'XTickLabel', [44.3:49.3], 'FontSize', 18); +drawnow, snapnow + +%% CANlab 2018 Combined Atlas object +% An atlas-class object is a specialized subclass of fmri_data that stores +% information about a series of parcels, or pre-defined regions, and in +% some cases the probalistic maps that underlie the final parcellation. +% +% A list of pre-defined atlases is contained in the function |load_atlas|. +% The default atlas for some CANlab functions is the "CANlab combined 2018" +% This was created by Tor Wager from other published atlases. It includes: +% +% * Glasser 2016 Nature 180-region cortical parcellation (in MNI space, not indivdidualized) +% * Pauli 2016 PNAS basal ganglia subregions +% * Amygdala/hippocampal and basal forebrain regions from SPM Anatomy Toolbox +% * Thalamus regions from the Morel thalamus atlas +% * Subthalamus/Basal forebrain regions from Pauli "reinforcement learning" HCP atlas +% * Diedrichsen cerebellar atlas regions +% * Multiple named brainstem nuclei localized based on individual papers +% * Additional Shen atlas parcels to cover areas (esp. brainstem) not otherwise named +% +% References for the corresponding papers are stored in the atlas object, and +% printed in tables generated with the region.table() method. + +atlas_obj = load_atlas('canlab2018'); + +% Create a figure showing cortical parcels +create_figure('iso'); isosurface(atlas_obj, 'alpha', .5); +drawnow, snapnow + +% Create a figure of the thalamus regions only +thal = select_atlas_subset(atlas_obj, {'Thal'}); +create_figure('iso'); isosurface(thal); +hh = addbrain('hires'); +set(hh, 'FaceAlpha', 0.05); +set(gca, 'XLim', [-40 40], 'YLim', [-40 20]); +drawnow, snapnow + + +%% Kragel 2015 emotion category-predictive patterns +% +% This loads a series of 7 multivariate patterns from Kragel et al. 2015. +% These patterns were developed to predict emotion categories. +% For our purposes here, we'll just treat them as data images. + +[test_dat, names] = load_image_set('kragelemotion', 'noverbose'); % loads 7 masks from Kragel et al. +disp(char(names)) + +% Plot spatial correlations +plot_correlation_matrix(test_dat.dat, 'names', names); +drawnow, snapnow + +%% Subcortical and brainstem regions +% |canlab_load_ROI| has a registry of many named regions derived from previous +% studies. This is particularly useful for loading a subcortical ROI and +% visualizing it or applying it to data. + +help canlab_load_ROI + +vmpfc = canlab_load_ROI('vmpfc'); +nac = canlab_load_ROI('nac'); +pag = canlab_load_ROI('pag'); + +create_figure('brain'); +isosurface(vmpfc, 'color', {[1 0 1]}); +isosurface(nac, 'color', {[1 .5 0]}); +isosurface(pag, 'color', {[1 0 0]}); +bhan = addbrain('hires right'); +view(-119, 5); +drawnow, snapnow + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m new file mode 100644 index 00000000..c3658b4a --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m @@ -0,0 +1,224 @@ +%% Perform a voxel-wise t-test +% +% In this example we perform a second level analysis on first level +% statistical parametric maps. Specifically we use a t-test to obtain an +% ordinary least squares estimate for the group level parameter. +% +% The example uses the emotion regulation data provided with +% CANlab_Core_Tools. This dataset consists of a set of contrast images for +% 30 subjects from a first level analysis. The contrast is [reappraise neg vs. look neg], +% reappraisal of negative images vs. looking at matched negative images. +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% By the end of this example we will have regenerated the results of figure +% 2A of this paper. + +%% Executive summary of the whole analysis +% +% As a summary, here is a complete set of commands to load the data and run +% the entire analysis: +% +% img_obj = load_image_set('emotionreg'); % Load a dataset +% t = ttest(img_obj); % Do a group t-test +% t = threshold(t, .05, 'fdr', 'k', 10); % Threshold with FDR q < .05 and extent threshold of 10 contiguous voxels +% +% % Show regions and print a table with labeled regions: +% montage(t); drawnow, snapnow; % Show results on a slice display +% r = region(t); % Turn t-map into a region object with one element per contig region +% table(r); % Print a table of results using new region names +% +% Here a graphical description of what it does: +% +% <> +% +% Now, let's walk through it step by step, with a few minor additions. + +%% Load sample data + +% We can load this data using load_image_set(), which produces an fmri_data +% object. Data loading exceeds the scope of this tutorial, but a more +% indepth demonstration may be provided by load_a_sample_dataset.m + +[image_obj, networknames, imagenames] = load_image_set('emotionreg'); + +%% Do a t-test + +% Voxel-wise tests estimate parameters for each voxel independently. The +% fmri_data/ttest function performs this just like the native matlab +% ttest() function, and similarly returns a p-value and t-stat. Unlike the +% matlab ttest() function, fmri_data/ttest performs this for every voxel in +% the fmri_data object. All voxels are evaluated independently, but all are +% evaluated nonetheless. A statistic_image is returned. + +t = ttest(image_obj); + +%% Visualize the results +% There are many options. See methods(statistic_image) and methods(region) + +orthviews(t) +drawnow, snapnow; + +% As can be seen, the ttest() result is unthresholded. The threshold +% function can be used to apply a desired alpha level using any of a number +% of methods. Here FDR is used to control for alpha=0.05. Note that no +% information is erased when performing thresholding on a statistic_image. + +t = threshold(t, .05, 'fdr'); +orthviews(t) +drawnow, snapnow; + +% Many neuroimaging packages (e.g., SPM and FSL) do one-tailed tests +% (with one-tailed p-values) and only show you positive effects +% (i.e., relative activations, not relative deactivations). +% All the CANlab tools do two-sided tests, report two-tailed p-values. +% By default, hot colors (orange/yellow) will show activations, and cool +% colors (blues) show deactivations. + +% We can also apply an "extent threshold" of > 10 contiguous voxels: + +t = threshold(t, .05, 'fdr', 'k', 10); +orthviews(t) +drawnow, snapnow; + +% montage is another visualization method. This function may require a +% relatively large amount of memory, depending on the resolution of the +% anatomical underlay image you use. We recommend around 8GB of free memory. +% +create_figure('montage'); axis off +montage(t) +drawnow, snapnow; + +%% Print a table of results + +% First, we'll have to convert to another object type, a "region object". +% This object groups voxels together into "blobs" (often of contiguous +% voxels). It does many things that other object types do, and inter-operates +% with them closely. See methods(region) for more info. + +% Create a region object called "r", with contiguous blobs from the +% thresholded t-map: + +r = region(t); + +% Print the table: + +table(r); + +%% More on getting help + +% Help info is available on https://canlab.github.io/ +% and in the CANlab help examples repository: https://github.com/canlab/CANlab_help_examples +% This contains a series of walkthroughs, including this one. + +% You can get help for the main object types by typing this in Matlab: +% doc fmri_data (or doc any other object types) + +% Also, more detailed help is available for each function using the "help" command. +% You can get help and options for any object method, like "table". But +% because the method names are simple and often overlap with other Matlab functions +% and toolboxes (this is OK for objects!), you will often want to specify +% the object type as well, as follows: + +% help region.table + +%% Write the t-map to disk + +% Now we need to save our results. You can save the objects in your +% workspace or you can write your resulting thresholded map to an analyze +% file. The latter may be useful for generating surface projections using +% Caret or FreeSurfer for instance. +% +% Thresholding did not actually eliminate nonsignificant voxels from our +% statistic_image object (t). If we simply write out that object, we will +% get t-statistics for all voxels. + +t.fullpath = fullfile(pwd, 'example_t_image.nii'); +write(t, 'overwrite') + +% If we use the 'thresh' option, we'll write thresholded values: +% (write() does not overwrite an existing file unless asked, and we just +% wrote this path above, so pass 'overwrite' here too.) +write(t, 'thresh', 'overwrite') + +t_reloaded = statistic_image(t.fullpath, 'type', 'generic'); +orthviews(t_reloaded) + + +%% Some useful optional extensions + +figure; surface(t); +figure; surface(t, 'foursurfaces'); +figure; surface(t, 'coronal_slabs_5'); + +table_of_atlas_regions_covered(t); + + +%% Explore on your own +% +% 1. Click around the thresholded statistic image. What lobes of the brain +% are most of the results in? What can we say about areas that are not +% active, if anything? +% +% 2. Why might some of the results appear to be outside the brain? What +% does this mean for the validity of the analysis? Should we consider these +% areas in our interpretation, or should they be "masked out" (excluded)? +% +% That's it for this section!! + +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> +%% +% % Bogdan Petre on 9/14/2016, updated by Tor Wager July 2018, Jan 2020 +% \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m new file mode 100644 index 00000000..4574e9c7 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m @@ -0,0 +1,319 @@ +%% Using the atlas object for Region of Interest (ROI) analyses + +%% About the atlas object +% An atlas-class object is a specialized subclass of fmri_data that stores +% information about a series of parcels, or pre-defined regions, and in +% some cases the probalistic maps that underlie the final parcellation. +% +% Some common uses of atlas objects include: +% * Labeling regions by best-matching atlas regions, as in region.table() +% * Analysis within specific ROIs, or on ROI averages +% * Defining regions for connectivity and graph theoretic analyses +% +% A list of pre-defined atlases is contained in the function |load_atlas|. +% The default atlas for some CANlab functions is the "CANlab combined 2018" +% This was created by Tor Wager from other published atlases. It includes: +% +% * Glasser 2016 Nature 180-region cortical parcellation (in MNI space, not indivdidualized) +% * Pauli 2016 PNAS basal ganglia subregions +% * Amygdala/hippocampal and basal forebrain regions from SPM Anatomy Toolbox +% * Thalamus regions from the Morel thalamus atlas +% * Subthalamus/Basal forebrain regions from Pauli "reinforcement learning" HCP atlas +% * Diedrichsen cerebellar atlas regions +% * Multiple named brainstem nuclei localized based on individual papers +% * Additional Shen atlas parcels to cover areas (esp. brainstem) not otherwise named +% +% References for the corresponding papers are stored in the atlas object, and +% printed in tables generated with the region.table() method. +% +% There are many methods for atlas, including all the fmri_data methods. +% You can see those with |methods(atlas)|. For example: +% +% Using the |select_atlas_subset()| method: +% You can select any subset of atlas regions by name or number, and return +% them in a new atlas object. You can also 'flatten' regions, combining +% them into a single new region. +% +% % Using the |extract_data()| method: +% You can extract average data from every image in an atlas, for a series of target images. +% +% We will explore these here. + +%% load an atlas +% 'atlas' objects are a class of objects specially designed for brain +% atlases. Here is more information on this class (also try |doc atlas|) + +help atlas + +% The function load_atlas in the CANlab toolbox loads a number of named +% atlases included with the toolbox. Here is a list of named atlases: + +help load_atlas + +% Now load the "CANlab combined 2018" atlas: +atlas_obj = load_atlas('canlab2018_2mm'); + + +%% visualize the atlas regions + +orthviews(atlas_obj); + +o2 = montage(atlas_obj); + +drawnow, snapnow +%% Selecting regions of interest +% There are several ways of selecting regions of interest. You can do it: +% * By name or part of name (in the .labels, .labels_2, or other fields) +% * By number +% * Graphically, or by entering a coordinate and a search radius + +%% +% *Select regions graphically or by coordinate* + +% First bring up the orthviews display: +orthviews(atlas_obj); + +% Now click on a coordinate in the center of vmPFC +% We'll use the SPM code below to do this here: + +spm_orthviews('Reposition', [0 38 -11]) + +% Running this will create a subset atlas containing only regions with centers +% within 20 mm of that coordinate. + +[obj_within_20mm,index_vals] = select_regions_near_crosshairs(atlas_obj, 'thresh', 20); + +orthviews(obj_within_20mm); +drawnow, snapnow + +% To reproduce this later in a script, you can enter the coordinates as +% input instead: + +[obj_within_20mm,index_vals] = select_regions_near_crosshairs(atlas_obj, 'coords', [0 38 -11], 'thresh', 20); + +% This function returns index values, so you can use these directly to +% specify the regions in select_atlas_subset as well. + +find(index_vals) + +test_obj = select_atlas_subset(atlas_obj, find(index_vals)); + +% % This should yield the same map of vmPFC regions: orthviews(test_obj) + +%% +% *Select regions by name* +% Let's select all the regions in the thalamus. All regions are labeled in +% the atlas object, so we can select them by name. + +% Select all regions with "Thal" in the label: +thal = select_atlas_subset(atlas_obj, {'Thal'}) + +% Print the labels: +thal.labels + +% Select a few thalamus/epithalamus regions of interest: +thal = select_atlas_subset(atlas_obj, {'Thal_Intra', 'Thal_VL', 'Thal_MD', 'Thal_LGN', 'Thal_MGN', 'Thal_Hb'}); +thal.labels + +% Select all the regions with "Thal" in the label, and collapse them into a single region: +whole_thal = select_atlas_subset(atlas_obj, {'Thal'}, 'flatten'); + +%% Extract and analyze data from atlas regions +% + +%% +% First, we'll load a dataset to extract data from for an ROI analysis +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj +% This is an fmri_data object from the CANlab Core Tools repository for neuroimaging data analysis. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed +% +% Alternative sample datasets: +% -------------------------------------------- +% This dataset will take time to download from figshare if you haven't yet +% downloaded it. An alternative is to use the dataset included with the +% CANlab Core toolbox (though the results are less interesting for this +% example). +% image_obj = load_image_set('emotionreg'); + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + + +%% +% *Now, we'll extract data from each atlas region* +% - "r" is a region-class object. see >> help region +% - r(i).dat contains averages over voxels within region i. for n images, it contains an n x 1 vector with average data. +% - r(i).data contains an images x voxels matrix of all data within region i. + +r = extract_roi_averages(image_obj, thal); + +%% +% You could also extract data from the whole thalamus into a single +% variable using: +% |r = extract_roi_averages(data_obj, whole_thal);| + +%% +% An alternative is the extract_data method: +% |region_avgs = extract_data(thal, image_obj);| +% The values may differ slightly based on which image is being resampled to +% which. Resampling matches the voxels in two sets of images that may have different +% voxel sizes and voxel boundaries by interpolating values from one to +% match the space of the other. + +%% Do a group ROI analysis and make a "violin plot" (or barplot) for each region +% +% We'll use the |barplot_columns| function, which is a versatile function +% for plotting columns of matrices. It has many options for displaying or +% hiding individual data points, "violins" (distribution estimates), +% standard errors, and more. +% +% The |barplot_columns| function also does a basic one-sample t-test +% analysis on each column, which corresponds to the average of each region +% here. This ROI analysis can be used for inference about a region. +% +% With the pain dataset, for a real analysis we might worry that multiple +% images come from the same subject. We can't do stats across all the +% iamges without accounting for the correlations among images belonging to +% the same subject. But first, let's plot the averages and see waht we get: + +% Concatenate the region averages into a matrix: +roi_averages = cat(2, r.dat); + +% Get the labels for each region from the atlas object, \ +% and replace some characters to make them look a bit nicer: +thal_labels = format_strings_for_legend(thal.labels); + +% Now make the figure: +create_figure('Thalamus regions', 1, 2); + +barplot_columns(roi_averages, 'nofig', 'colors', scn_standard_colors(length(r)), 'names', thal_labels); +xlabel('Region') +ylabel('ROI activity') + +subplot(1, 2, 2) + +barplot_columns(roi_averages, 'nofig', 'colors', scn_standard_colors(length(r)), 'names', thal_labels, 'noind', 'noviolin'); +xlabel('Region') +ylabel('ROI activity') + +%% +% Look at the results. Do they make sense for a pain dataset? +% The sensory/nociceptive regions of the thalamus (VL and intralaminar) respond to painful +% stimuli, but not the regions in visual and auditory pathways (LGN, MGN). The +% habenula also responds. + +%% +% To do a proper ROI analysis, let's just select a subset of the images +% responding to one stimulus intensity, with one image per person. +% the stimulus temperatures are stored in +% |image_obj.additional_info.temperatures| +% So we can use basic Matlab code to select one temperature. + +wh = image_obj.additional_info.temperatures == 48; % A logical vector for 48 degrees +indx_list = find(wh); % a list of which images + +image_obj_48 = get_wh_image(image_obj, indx_list); % Get the images we want + +% You can try to extract data from the thalamus and repeat the t-test. + +%% Selecting multiple regions using custom label fields +% In addition to the .labels field, atlas objects have fields for additional +% labels (.labels_2, etc.). These can be used to group, and select, regions +% at multiple levels of a spatial hierarchy. For example, ?canlab2018_2mm? +% atlas has macro-network and anatomical group labels in the ?labels_2? field. +% You can use these to select all the regions associated with a particular +% resting-state network or group like the brainstem or basal ganglia. +% Here?s a code snippet to locate/pull out all the default mode network +% regions using the select_atlas_subset() method. This method can also return +% a vector that you can use to select elements of associated connectivity matrices. +% +% You can create and/or annotate your own atlas labels and add them to the +% .labels, labels_2, or other fields, and use select_atlas_subset to +% extract and manipulate them. + +atlas_obj = load_atlas('canlab2018_2mm'); + +% List all the macro-level labels in .labels_2 field +disp(unique(atlas_obj.labels_2)') + +% Pull out all regions with ?Def? (Default) as part of the label in labels_2 +atlas_obj = select_atlas_subset(atlas_obj, {'Def'}, 'labels_2'); +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); +% let's add a cortical surface for context +% We'll make the back surface (right) opaque, and the front (left) transparent +han = addbrain('hires right'); +set(han, 'FaceAlpha', 1); +han2 = addbrain('hires left'); +set(han2, 'FaceAlpha', 0.1); +axis image vis3d off +material dull +view(-88, 31); + +lightRestoreSingle +drawnow, snapnow + +%% Explore on your own +% +% 1. One of the ideas about emotion regulation is that positive appraisal involves +% activity increases in the nucleus accumbens (NaC) in the basal ganglia. +% Load the emotion regulation dataset from previous walkthroughs. +% Try to locate and extract data from this region, and do an ROI analysis. +% What do you see? +% Hint: The NAc is part of the ventral striatum. +% To extract the ventral striatum from the canonical CANlab atlas, try: +% NAc = select_atlas_subset(atlas_obj, {'V_Str'}); +% +% 2. try the t-test on 48 degree only images above. What do you see? Which +% regions show sigificant activation? +% +% 3. Try to identify which prefrontal regions show the greatest emotion +% regulation effect, and which show the greatest response during pain. +% Are they the same? + +% That's it for this section!! + + + + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m new file mode 100644 index 00000000..d05620fd --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m @@ -0,0 +1,163 @@ +%% Applying masks and writing image data to disk + +% One common operation in neuroimaging is masking. Masking is defining a set +% of voxels to include or exclude from analysis or presentation. Some uses +% include: +% +% * Apply a generous brain mask to store only in-brain voxels, saving space, +% (CANlab data objects do this by default, with a liberal mask) +% * Identify gray matter voxels to calculate global mean signal or other metrics +% * Identify voxels in ventricles or out-of-brain to try to isolate sources +% of noise or artifacts to remove +% * Identify gray matter voxels of interest to limit the space within which +% you correct for multiple comparisons, thus increasing power +% * Limiting presentation of results to an a priori set of ROIs +% * Limiting the scope of training or testing a multivariate pattern or +% predictive model, to help test hypotheses about which systems/regions +% contribute to prediction. +% +%% The Neuroimaging_pattern_masks Github repository and website +% +% For this walkthrough, we'll load some new image datasets to use as masks. +% A popular one is a set of "networks" developed from resting-state fMRI +% connectivity in 1,000 people. For our purposes, the "network" maps +% consist of a set of voxels that load most highly on each of 7 ICA +% components from the 1,000 Functional Connectomes project. +% These were published by Randy Buckner's lab in 2011 in three papers: +% Choi et al., Buckner et al., and Yeo et al. We'll use the cortical +% networks from Yeo et al. here. +% +% This set of images and other broadly useful sets of images are stored in +% the CANlab function |load_image_set| in a registry that you can access by +% name. Try |help load_image|set| for a list of images you can load by +% keyword. These datasets, if they can be stored on Github, are in this +% Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% To run this tutorial and many others in this series, you'll need to +% download the Neuroimaging_Pattern_Masks Github repository and add it to +% your Matlab path with subfolders. +% The Github site has three main types of datasets, shown here: +% +% <> +% + +%% Load and select a sample mask +% You will need the Neuroimaging_Pattern_Masks repository installed to find +% these images. + +[bucknermaps, mapnames] = load_image_set('bucknerlab'); % loads 7 masks from Yeo et al. + +%% +% This loads 7 images. try |orthviews(dat)| to see them all. +% What values do the images have? +% You can also see descriptive statistics for the images like this: + +desc = descriptives(bucknermaps); + +%% +% The variable |desc| is a structure that contains various descriptive +% stats. +% + +%% +% load_image_set will generally store and return names for each image, +% returned as outputs if requested (here in |mapnames|). Let's print +% those names: + +disp(char(mapnames)) + +%% +% Now we can select one image from those 7 to use as a mask. + +mask_image = bucknermaps.get_wh_image(1); + +%% +% Now we'll load a sample dataset to apply the mask to. +% This loads a series of 7 multivariate patterns from Kragel et al. 2016. +% These patterns were developed to predict emotion categories. +% For our purposes here, we'll just treat them as data images. + +test_dat = load_image_set('kragelemotion', 'noverbose'); % loads 7 masks from Kragel et al. + +% Try |orthviews(dat)| or use |montage| to see these patterns. +% This command shows the first 4 (a limit in |montage| to avoid messy +% displays): + +montage(test_dat); + +%% Apply a mask +% +% We can apply the mask using the apply_mask() method for fmri_data +% objects. It has various options, but the simplest one applies a single +% mask to an fmri_data object. As a rule, you should pass in the object you're +% operating on first, followed by additional modifying arguments/objects. +% Here, this means data object followed by mask object. + +test_dat = apply_mask(test_dat, mask_image); + +% Now let's view the new, masked dataset: + +montage(test_dat); + +% Now, the dataset contains only voxels in the "visual network", and +% subsequent analyses or operations would test only that. + +%% A few exercises +% +% 1. Here's a question to answer: What are some potential uses of the masking +% operation we did above? Why would you want to do it? +% +% 2. Try loading the emotion regulation dataset from previous tutorials, +% and do a one-sample t-test group analysis only within the +% "frontoparietal" network. Try thresholding the results with FDR +% correction for both the full dataset and the analysis within the +% restricted mask. How are the results different? Why? + +%% Writing image data in objects to Nifti files on disk +% +% canlabcore fmri data objects can be easily written to standard image file +% formats. The main function for this is write() +% +% This walkthrough shows how to do this for 1) a mask image, and 2) a +% statistic image. + +%% +% Create a Nifti image called dummy.nii in the current folder. Write the +% mask image to this file. See help for write() for more options + +fname = 'dummy.nii'; % outname file name. the extension (.nii, .img) determines the format. + +% |write(mask_image, 'fname', fname);| + +%% +% If you've already done this, |write()| will return an error. This is a +% safeguard to keep users from inadvertently overwriting important files. +% To force overwriting the existing file, try: + +write(mask_image, 'fname', fname, 'overwrite'); + +%% +% Now let's read it back in from disk and view it again: + +orthviews(fmri_data(fname)) + + +%% Exercises +% +% 1. Select a sample image from |test_dat|, threshold it at an arbitrary threshold +% of absolute values > .0005, and write it to disk. Load it again from +% disk and display it. See the thresholding walkthrough for more help if +% needed. +% Hint: help fmri_data.threshold will give you help on how to do this +% e.g., Retain voxels with absolute value > .0001 +% obj = threshold(obj, [-.0001 0001], 'raw-outside') +% +% 2. Do the one-sample t-test on the emotion regulation dataset (see the +% t-test walkthrough if needed), and threshold the results at 0.05 +% FDR-corrected. Write the thresholded results map to disk. Reload it from +% disk and display it using the montage() method. diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m new file mode 100644 index 00000000..ad99f6bf --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m @@ -0,0 +1,468 @@ +%% 3-D volume visualization + +%% About volume visualisation +% +% Creating 3-D renderings of brain data and results can help localize brain +% activity and show patterns in a way that allows meaningful comparisons +% across studies. +% +% The CANlab object-oriented tools have methods built in for rendering +% volume data on brain surfaces. They generally use Matlab's isosurface and +% isocaps functions. +% + +%% Prepare sample statistical results maps +% +% For this walkthrough, we'll use the "Pain, Cognition, Emotion balanced N = 270 dataset" +% from Kragel et al. 2018, Nature Neuroscience. +% More information about this dataset and how to download it is in an +% earlier walkthrough on "loading datasets" +% +% The code below loads it; but you could also use any other brain dataset +% for this purpose (e.g., the 'emotionreg' dataset in load_image_set. +% +% Our goal is to select subject-level images corresponding to a particular task type and +% do a t-test on these, and visualize the results in various 3-D ways. + +[test_images, names] = load_image_set('kragel18_alldata', 'noverbose'); + +% This field contains a table object with metadata for each image: +metadata = test_images.metadata_table; +metadata(1:5, :) % Show the first 5 rows + +% Here are the 3 domains: +disp(unique(metadata.Domain)) + +% Find all the images of "Neg_Emotion" and pull those into a separate +% fmri_data object: + +wh = find(strcmp(metadata.Domain, 'Neg_Emotion')); +neg_emo = get_wh_image(test_images, wh); + +% Do a t-test on these images, and store the results in a statistic_image +% object called "t". Threshold at q < 0.001 FDR. This is an unusual +% threshold, but we have lots of power here, and want to restrict the +% number of significant voxels for display purposes below. + +t = ttest(neg_emo, .001, 'fdr'); + +% Find the images for other conditions and save them for later: + +t_emo = t; + +wh = find(strcmp(metadata.Domain, 'Cog_control')); +cog_control = get_wh_image(test_images, wh); +t_cog = ttest(cog_control, .001, 'fdr'); + +wh = find(strcmp(metadata.Domain, 'Pain')); +pain = get_wh_image(test_images, wh); +t_pain = ttest(pain, .001, 'fdr'); + + +%% The surface() method +% Most data objects, including fmri_data, statistic_image, and region objects, +% have a surface() method. Entering it with no arguments creates a surface +% or series of surfaces (depending on the object) +% +% The surface() method lets you easily render blobs +% on a number of pre-set choices for brain surfaces or 3-D cutaway surfaces. + +% This generates a series of 6 surfaces of different types, including a +% subcortical cutaway +surface(t); +drawnow, snapnow + +% This generates a different plot +create_figure('lateral surfaces'); +surface_handles = surface(t, 'foursurfaces', 'noverbose'); +snapnow + +% You can see more options with: +% |help statistic_image.surface| or |help t.surface| + +% You can also use the |render_on_surface()| to method to change the colormap: + +render_on_surface(t, surface_handles, 'colormap', 'summer'); +snapnow + +%% Using the fmridisplay object to create a registry of montages/surfaces +% The fmridisplay object is an object class that stores information about a +% set of display items you have created. It stores information about a +% series of montages (in the same figure or different figures) and surfaces +% along with their handles. +% +% You can use this to render activations on a whole set of slices/surfaces +% at once. You can create your own custom views, or use pre-packaged sets +% of montages/surfaces created by the function |canlab_results_fmridisplay| +% +% The code snippet below displays our t-map on a relatively complete series +% of slices and surfaces, so you can see the "whole picture". + +my_display_obj = canlab_results_fmridisplay([], 'montagetype', 'full'); +my_display_obj = addblobs(my_display_obj, region(t)); +snapnow + +%% +% You can see some prepackaged sets in |help canlab_results_fmridisplay| +% These include: +% 'full' Axial, coronal, and saggital slices, 4 cortical surfaces +% 'compact' Midline saggital and two rows of axial slices [the default] +% 'compact2' A single row showing midline saggital and axial slices +% 'multirow' A series of 'compact2' displays in one figure for comparing different images/maps side by side +% 'regioncenters' A series of separate axes, each focused on one region +% +% In addition, you can pass in a series of optional keywords that will be passed onto the +% |addblobs| object method that controls rendering of the blobs on slices. +% (They don't all work for surfaces). These control colors, inclusion of an +% *outline* around each activation blob/region, and more. +% See |help fmridisplay.addblobs| for more information. +% Here is an example using a few of the options. + +my_display_obj = canlab_results_fmridisplay(t, 'compact2', 'color', [1 0 0], 'nooutline', 'trans'); + +%% +% After creating the display, other maps can be added. +% The addblobs() method requires a region object, so we use region(t) to convert to a region object. +% This shows us regions activated during negative emotion in red and pain +% in blue. + +my_display_obj = addblobs(my_display_obj, region(t_pain), 'color', [0 0 1], 'nooutline', 'trans'); +snapnow + +%% +% my_display_obj is an fmridisplay-class object. It has its own methods: + +methods(my_display_obj) + +% One particularly useful feature is that you can add or _remove_ activation blobs. +% This allows you to create a canonical set of views and then render one +% activation map on them, remove the blobs, and render other maps on the +% same set of display items. You can add points, spheres (e.g., from +% different studies), and more. +% +% When you operate on an fmridisplay object, pass it back out as an output +% argument so that the information you updated will be available. + +my_display_obj = removeblobs(my_display_obj); +snapnow + +%% Using addbrain to load or build surfaces +% You can also build any surface object you want and render colors onto it. +% This is a very versatile method. +% +% Matlab creates handles to surfaces (and +% other graphics objects). You can use the handles to manipulate the +% objects in a many ways. They have, for example, color, line, face color, +% edge color, and many more properties. If |h| is a figure handle, |get(h)| +% shows a list of its properties. |set(h, 'MyPropertyName', myvalue)| sets +% the property |'MyPropertyName'| to |myvalue|. +% +% Here, we'll use |addbrain| to build up a set of surfaces with handles, +% and then render colored blobs on them using |surface()|. Instead of the +% object method, you can also use its underlying function: |cluster_surf| +% +% Try |help addbrain| for a list of surfaces, and |help cluster_surf| for +% other color/display/scaling options. + +% Let's build a surface by starting with a group of thalamic nuclei, and +% adding the parabrachial complex and the red nucleus. +% Then, we'll render our t-statistics in colors onto those surfaces. + +create_figure('cutaways'); axis off +surface_handles = addbrain('thalamus_group'); + +surface_handles = [surface_handles addbrain('pbn')]; +surface_handles = [surface_handles addbrain('rn')]; +surface_handles = [surface_handles addbrain('pag')]; + +drawnow, snapnow + +%% +% Now render the statistic image stored in _t_ onto those surfaces. +% All non-activated areas turn gray. + +render_on_surface(t, surface_handles, 'colormap', 'summer'); +set(surface_handles, 'FaceAlpha', .85); +set(surface_handles(5:6), 'FaceAlpha', .15); % Make brainstem and thalamus shell more transparent +drawnow, snapnow + +% Emphasize positive values, in a hot colormap +render_on_surface(t, surface_handles, 'clim', [0.01 3]); +drawnow, snapnow + +% Show positive and negative values, in a hot colormap +render_on_surface(t, surface_handles, 'clim', [-3 3]); +drawnow, snapnow + +%% +% Now render the statistic image stored in _t_ onto those surfaces. +% All non-activated areas turn gray. + +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +%% +% Transparent surfaces are a little hard to see. +% We can use Matlab's powerful handle graphics system to inspect and change +% all kinds of properties. See the Matlab documentation for more details. +% Let's just make the surfaces solid: + +set(surface_handles, 'FaceAlpha', 1); +drawnow, snapnow + +%% +% These surfaces can be rotated too (or zoom in/out, pan, set lighting, +% etc.) Let's shift the angle. + +view(222, 15); % rotate +camzoom(0.8); % zoom out + +% A helpful CANlab function to re-set the lighting after rotating is: + +lightRestoreSingle; + +drawnow, snapnow + +%% +% addbrain also has keywords for composites of multiple surfaces + +create_figure('cutaways'); axis off +surface_handles = addbrain('limbic'); +t.surface('surface_handles', surface_handles, 'noverbose'); + +%% Removing blobs and re-adding new colors +% There is a special command to re-set the surface colors to gray. +% Then we can add a t-test for different contrasts without re-drawing +% surfaces. + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Cognitive Control'); +t = ttest(cog_control, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Pain'); +t = ttest(pain, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Negative Emotion'); +t = ttest(neg_emo, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + + +%% Prepackaged cutaway surfaces +% Addbrain has many 3-D surfaces, including cutaways that show various 3-D +% sections. The code below visualizes the various options for cutaways. +% +% You can also pass any of these keywords into surface() to generate the +% surface and render colored blobs. +% Let's render these with and without an unthresholded pain map + +t = ttest(pain, 1, 'unc'); +t_thr = threshold(t, .05, 'fdr'); + +keywords = {'left_cutaway' 'right_cutaway' 'left_insula_slab' 'right_insula_slab' 'accumbens_slab' 'coronal_slabs' 'coronal_slabs_4' 'coronal_slabs_5'}; + +for i = 1:length(keywords) + + create_figure('cutaways'); axis off + + surface_handles = surface(t_thr, keywords{i}); + + % This essentially runs the code below: + + % surface_handles = addbrain(keywords{i}, 'noverbose'); + % render_on_surface(t, surface_handles, 'clim', [-7 7]); + + % Alternative: This command creates the same surfaces: + % surface_handles = canlab_canonical_brain_surface_cutaways(keywords{i}); + % render_on_surface(t, surface_handles, 'clim', [-7 7]); + + drawnow, snapnow + +end + +%% +% Plot thresholded and unthresholded maps side by side + +create_figure('cutaways',1, 2); axis off + +surface_handles = surface(t, 'right_cutaway'); +title('Pain, Unthresholded') + +subplot(1, 2, 2); +surface_handles = surface(t_thr, 'right_cutaway'); +title('Pain, thresholded') + +drawnow, snapnow + + +%% +% *Controlling options in render_on_surface()* +% +% When the t map has positive and negative values, render_on_surface creates a special +% bicolor split colormap that has warm colors for positive vals, and cool colors +% for negative vals. You can set the color limits (here, in t-values) + +create_figure('cutaways'); axis off +surface_handles = surface(t_thr, 'coronal_slabs'); + +render_on_surface(t_thr, surface_handles, 'clim', [-4 4]); +drawnow, snapnow + +% You can set the colormap to any Matlab colormap: +render_on_surface(t_thr, surface_handles, 'colormap', 'summer'); +drawnow, snapnow + +% If your image is positive-valued only (or negative-valued only), a +% bicolor split map will not be created: + +t = threshold(t_thr, [2 Inf], 'raw-between'); +render_on_surface(t, surface_handles, 'colormap', 'winter', 'clim', [2 6]); +drawnow, snapnow + +%% Creating isosurfaces +% The fmri_data.isosurface() method allows you to create and save +% isosurfaces for any image, parcel, or blob. +% +% You can also save the meshgrid output and surface structure that lets you +% render this surface easily in the future. +% +% Here's a simple example visualizing all of the parcels in an atlas as 3-D +% blobs: + +t = ttest(pain, .01, 'unc'); + +atlas_obj = load_atlas('thalamus'); + +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); + +% Now we'll set some lighting and figure properties +axis image vis3d off +material dull +view(210, 20); +lightRestoreSingle + +drawnow, snapnow + +% You can render blobs on these surfaces too. + +render_on_surface(t, surface_handles, 'colormap', 'hot'); +drawnow, snapnow + + +%% +% Let's load another atlas, and pull out all the "default mode" regions +% We'll use isosurface to visualize these + +atlas_obj = load_atlas('canlab2018_2mm'); +atlas_obj = select_atlas_subset(atlas_obj, {'Def'}, 'labels_2'); + +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); + +% let's add a cortical surface for context +% We'll make the back surface (right) opaque, and the front (left) transparent +han = addbrain('hires right'); +set(han, 'FaceAlpha', 1); + +han2 = addbrain('hires left'); +set(han2, 'FaceAlpha', 0.1); + +axis image vis3d off +material dull +view(-88, 31); +lightRestoreSingle + +drawnow, snapnow + +%% Using isosurface to create a custom surface +% You can build a cutaway surface or a series of them with the isosurface method, from +% any suitable image. The image used here is a mean anatomical image that +% renders nicely. +% As always, you can use Matlab's rendering to change the look of the surfaces + +anat = fmri_data(which('keuken_2014_enhanced_for_underlay.img'), 'noverbose'); + +create_figure('cutaways'); axis off + +surface_handles = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [-Inf 0], 'YLim', [-30 Inf], 'ZLim', [-Inf 40]); + +render_on_surface(t, surface_handles, 'colormap', 'hot'); + +view(-127, 33); +set(surface_handles, 'FaceAlpha', .85); +snapnow + +%% Explore on your own +% +% 1. Try to create your own custom brain surface to visualize one of the 3 +% statistic image maps we've been working on. Can you make a plot that +% really shows the important results in a clear way? +% +% 2. Try exploring movie_tools to create a movie where you rotate, zoom, +% and change the transparency of your surfaces in some way. +% You can write movie files to disk for use in presentations, too. +% +% 3. Try rendering another atlas, or subset of an atlas, with isosurface() +% Maybe you can pull out and render all the "ventral attention" network +% regions. +% +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m new file mode 100644 index 00000000..88c13409 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m @@ -0,0 +1,561 @@ +%% Regression Walkthrough + +%% +% This walkthrough works through multiple aspects of a basic 2nd-level (group) +% analysis. We are interested in activity related to emotion regulation +% for a contrast comparing reappraisal of negative images vs. looking at +% negative images without reappraisal, i.e., a [Reappraise Negative - Look Negative] +% contrast. This walkthrough will regress brain +% contrast values in each voxel (y) onto individual differences in +% "Reappraisal Success", defined as the drop in negative emotion ratings +% for the [Reappraise Negative - Look Negative] contrast. +% +% We are interested in both the brain correlates of individual differences +% and the group-average [Reappraise Negative - Look Negative] effect in the +% brain, which is the intercept in our regression model. + +%% +% This walkthrough includes multiple considerations you'll need when doing +% a practical data analysis. This goes far beyond just running a model! +% It's also important to understand the data and results. +% Understanding the data allows you to make informed a priori choices about +% the analysis and avoid "P-Hacking". Understanding the results helps +% ensure they are valid, and helps you interpret them. +% +% A) Understanding the data +% +% # Loading the dataset into an fmri_data object +% # Examining the brain data and checking +% +% outliers +% global signal +% which voxels are analyzed (coverage) +% homogeneity across subjects (images) +% +% # Examining the behavioral (individual differences) data and leverages +% # Use the info gained to make data preparation and analysis decisions +% +% Masking +% Transformation of the y variable +% Rescaling/normalizing image data +% Robust vs. OLS estimation +% +% B) Fitting the model +% +% # Running the regression and writing results images to disk +% +% C) Understanding the results +% +% # Explore variations in statistical threshold +% # Compare with other studies by comparing with a meta-analysis mask +% # Explore the impact of variations in data scaling/transformation +% # Localizing results, labeling regions, and making tables using a brain atlas +% # Extract and plot data from regions of interest +% # Explore the impact of plotting data in biased vs. unbiased regions of interest +% +% D) Multivariate prediction (overlaps with later walkthroughs) +% # Extra: Multivariate prediction within a set of unbiased ROIs + +%% Part 1: Load sample data +% Load a sample dataset, including a standard brain mask +% And a meta-analysis mask that provides a priori regions of interest. +% --------------------------------------------------------------- + +% Load sample data using load_image_set(), which produces an fmri_data +% object. Data loading exceeds the scope of this tutorial, but a more +% indepth demosntration may be provided by canlab_help_2_load_a_sample_dataset.m + +% These are [Reappraise - Look Neg] contrast images, one image per person + +[image_obj, networknames, imagenames] = load_image_set('emotionreg'); + +% Summarize and print a table with characteristics of the data: +desc = descriptives(image_obj); + +% Load behavioral data +% This is "Reappraisal success", one score per person, in our example +% If you do not have the file on your path, you will get an error. + +beh = importdata('Wager_2008_emotionreg_behavioral_data.txt') +success = beh.data(:, 2); % Reappraisal success + +% Load a mask that we would like to apply to analysis/results + +mask = which('gray_matter_mask.img') + +maskdat = fmri_data(mask, 'noverbose'); + +% Map of emotion-regulation related regions from neurosynth.org +metaimg = which('emotion regulation_pAgF_z_FDR_0.01_8_14_2015.nii') + +metamask = fmri_data(metaimg, 'noverbose'); + +%% Part 2: Examine the setup and coverage + +% Visualize the mask +% --------------------------------------------------------------- +% BLOCK 2: Check mask + +% This is an underlay brain with three separate montages so we can compare images: + +o2 = canlab_results_fmridisplay([], 'noverbose', 'multirow', 3); +drawnow, snapnow; + +% This is a basic gray-matter mask we will use for analysis: +% It can help reduce multiple comparisons relative to whole-image analysis +% but we should still look at what's happening in ventricles and +% out-of-brain space to check for artifacts. + +o2 = addblobs(o2, region(maskdat), 'wh_montages', 1:2); + +% Add a title to montage 2, the axial slice montage: +[o2, title_han] = title_montage(o2, 2, 'Mask image: Gray matter with border'); +set(title_han, 'FontSize', 18) + +% Visualize summary of brain coverage +% --------------------------------------------------------------- +% Check that we have valid data in most voxels for most subjects +% Always look at which voxels you are analyzing. The answer may surprise +% you! Some software programs mask out voxels in individual analyses, +% which may then be excluded in the group analysis. + +% Show summary of coverage - how many images have non-zero, non-NaN values in each voxel + +o2 = montage(region(desc.coverage_obj), o2, 'trans', 'transvalue', .3, 'wh_montages', 3:4); + +[o2, title_han] = title_montage(o2, 4, 'Coverage: Voxels included in images'); +set(title_han, 'FontSize', 18) + +% Visualize meta-analysis mask +% --------------------------------------------------------------- +o2 = montage(region(metamask), o2, 'colormap', 'wh_montages', 5:6); + +[o2, title_han] = title_montage(o2, 6, 'Meta-analysis mask'); +set(title_han, 'FontSize', 18) + +% There are many other alternatives to the montage - e.g.,: +% orthviews(desc.coverage_obj, 'continuous'); + +%% Part 3: Examine the brain data +% This is a really important step to understand what your images look like! +% Good data tend to produce good results. + +% The method plot() for an fmri_data object shows a number of plots that +% are helpful in diagnosing some problems. + +plot(image_obj); + +% Check histograms of individual subjects for global shifts in contrast values + +% The 'histogram' object method will allow you to examine a series of +% images in one panel each. See help fmri_data.histogram for more options, +% including breakdown by tissue type. + +hist_han = histogram(image_obj, 'byimage', 'color', 'b'); + +% This shows us that some of the images do not have the same mean as the +% others. This is fairly common, as individual subjects can often have +% global artifacts (e.g., task-correlated head motion or outliers) that +% influence the whole contrast image, even when baseline conditions are +% supposed to be "subtracted out". +% +% It suggests that we may want to do an outlier analysis and/or standardize +% the scale of the images. We'll return to this below. + +% If we want to visualize this separated by tissue class, try this: + +create_figure('histogram by tissue type') +hist_han = histogram(image_obj, 'byimage', 'by_tissue_type'); + +% Our previous peek at the data suggested that subjects 6 and 16 have +% whole-brain shifts towards acftivation and deactivation, respectively. + +% This plot shows us that there is a substantial shift in gray matter values +% across the brain for some individuals. +% Also, the means across gray and white matter are correlated. So are the standard deviations. +% We might consider a rescaling or nonparametric procedure (e.g., ranking) +% to standardize the images. + +% Some options are here: +help image_obj.rescale + + +%% Part 4: Examine the behavioral/task predictor(s) +% Examine predictor distribution and leverages +% Leverage is a measure of how much each point influences the regression +% line. The more extreme the predictor value, the higher the leverage. +% Outliers will have very high leverage. High-leverage behavioral observations +% can strongly influence, and sometimes invalidate, an analysis. + +X = scale(success, 1); X(:, end+1) = 1; % A simple design matrix, behavioral predictor + intercept +H = X * inv(X'* X) * X'; % The "Hat Matrix", which produces fits. Diagonals are leverage + +create_figure('levs', 2, 1); +plot(success, 'o', 'MarkerFaceColor', [0 .3 .7], 'LineWidth', 3); +set(gca, 'FontSize', 24); +xlabel('Subject number'); +ylabel('Reappraisal success'); + +subplot(2, 1, 2); +plot(diag(H), 'o', 'MarkerFaceColor', [0 .3 .7], 'LineWidth', 3); +set(gca, 'FontSize', 24); +xlabel('Subject number'); +ylabel('Leverage'); + +drawnow, snapnow + +% The distribution suggests that there are some high-leverage values to +% watch out for. They will influence the results disproportionately at every voxel in the brain! +% Ranking predictors and/or robust regression are good +% mitigation/prevention procedures in many cases. + +%% Part 5: Make data prep and analysis decisions +% Before we've "peeked" at the results, let's make some sensible decisions +% about what to do with the analysis. We'll compare the results to the +% "standard analysis" if we hadn't examined the data at all at end. + +% Choices: + +% 1. Masking: Apply a gray-matter mask. This will limit the area subjected to multiple-comparisons +% correction to a meaningful set. I considered {gray matter, no mask} + +image_obj = apply_mask(image_obj, maskdat); % Mask - analyze gray-matter + +% 2. Image scaling: I considered arbitrary rescaling of images, e.g., by +% the l2norm or z-scoring images. Z-scoring can cause artifactual +% de-activations when we are looking at group effects. L2-norming is +% sensible if there are bad scaling problems... but I'd rather minimize ad +% hoc procedures applied to the data. Ranking each voxel is also sensible, and +% if the predictor is ranked as well, this would implement Spearman's correlations. +% I considered {none, l2norm images, zscoreimages, rankvoxels}, and decided +% on ranking voxels. For the group average (intercept), this is not -- I'll use robust regression to help mitigate issues. + +image_obj = rescale(image_obj, 'rankvoxels'); % rescale images within this mask (relative activation) + +% 3. Predictor values: I considered dropping subject 16, which is a global +% outlier and also has high leverage, or windsorizing values to 3 sd. +% I also considered ranking predictor values, a very sensible choice. +% So I considered {do nothing, drop outlier, windsorize outliers, rank values} +% But robust regression (below) is partially redundant with this. I ultimately decided +% on ranking. + +% This is done below: +% image_obj.X = scale(rankdata(image_obj.X), 1); + +% 4. Regression method: I considered {standard OLS and robust IRLS +% regression}. I chose robust regression because it helps mitigate the +% effects of outliers. However, this is partly redundant with ranking, so +% it would have been sensible to choose ranking OR robust regression. + +%% Part 6: Run the regression analysis and get results +% The regress method takes predictors that are attached in the object's X +% attribute (X stands for design matrix) and regresses each voxel's +% activity (y) on the set of regressors. +% This is a group analysis, in this case correlating brain activity with +% reappraisal success at each voxel. + +% .X must have the same number of observations, n, in an n x k matrix. +% n images is the number of COLUMNS in image_obj.dat + +% mean-center success scores and attach them to image_obj in image_obj.X +image_obj.X = scale(success, 1); + +image_obj.X = scale(rankdata(image_obj.X), 1); + +% runs the regression at each voxel and returns statistic info and creates +% a visual image. regress = multiple regression. +% % Track the time it takes (about 45 sec for robust regression on Tor's 2018 laptop) +tic + +% out = regress(image_obj); % this is what we'd do for standard OLS regression +out = regress(image_obj, 'robust'); + +toc + +% out has statistic_image objects that have information about the betas +% (slopes) b, t-values and p-values (t), degrees of freedom (df), and sigma +% (error variance). The critical one is out.t. +% out = +% +% b: [1x1 statistic_image] +% t: [1x1 statistic_image] +% df: [1x1 fmri_data] +% sigma: [1x1 fmri_data] + +% This is a multiple regression, and there are two output t images, one for +% each regressor. We've only entered one regressor, why two images? The program always +% adds an intercept by default. The intercept is always the last column of the design matrix + +% Image 1 <--- brain contrast values correlated with "success" +% Positive effect: 0 voxels, min p-value: 0.00001192 +% Negative effect: 0 voxels, min p-value: 0.00170529 +% Image 2 FDR q < 0.050 threshold is 0.003193 +% +% Image 2 <--- intercept, because we have mean-centered, this is the +% average group effect (when success = average success). "Reapp - Look" +% contrast in the whole group. +% Positive effect: 3133 voxels, min p-value: 0.00000000 +% Negative effect: 51 voxels, min p-value: 0.00024068 + +% Now let's try thresholding the image at q < .05 fdr-corrected. +t = threshold(out.t, .05, 'fdr'); + +% Select the predictor image we care about: (the 2nd/last image is the intercept) +t = select_one_image(t, 1); + +% ...and display on a series of slices and surfaces +% There are many options for displaying results. +% montage and orthviews methods for statistic_image and region objects provide a good starting point. + +o2 = montage(t, 'trans', 'full'); +o2 = title_montage(o2, 2, 'Behavioral predictor: Reappraisal success (q < .05 FDR)'); +snapnow + +% Or: +% orthviews(t) + + +%% Write the thesholded t-image to disk + +% % t.fullpath = fullfile(pwd, 'Reapp_Success_005_FDR05.nii'); +% % write(t) +% % +% % disp(['Written to disk: ' t.fullpath]) + +%% Part 7. Explore: Display regions at a lower threshold, alongside meta-analysis regions +% Display the results on slices another way: +% multi_threshold lets us see the blobs with significant voxels at the +% highest (most stringent) threshold, and voxels that are touching +% (contiguous) down to the lowest threshold, in different colors. +% This shows results at p < .001 one-tailed (p < .002 two-tailed), with 10 contiguous voxels. +% Contiguous blobs are more liberal thresholds, down to .02 two-tailed = +% .01 one tailed, are shown in different colors. + +% Set up a display object with two separate brain montages: +o2 = canlab_results_fmridisplay([], 'multirow', 2); + +% Use multi-threshold to threshold and display the t-map: +% Pass the fmridisplay object with montages (o2) in as an argument to use +% that for display. Pass out the object with updated info. +% 'wh_montages', 1:2 tells it to plot only on the first two montages (one sagittal and one axial slice series) + +o2 = multi_threshold(t, o2, 'thresh', [.002 .01 .02], 'sizethresh', [10 1 1], 'wh_montages', 1:2); +o2 = title_montage(o2, 2, 'Behavioral predictor: Reappraisal success (p < .001 one-tailed and showing extent'); + +% Get regions from map from neurosynth.org +% Defined above: metaimg = which('emotion regulation_pAgF_z_FDR_0.01_8_14_2015.nii') + +r = region(metaimg); + +% Add these to the bottom montages +o2 = addblobs(o2, r, 'maxcolor', [1 0 0], 'mincolor', [.7 .2 .5], 'wh_montages', 3:4); + +o2 = title_montage(o2, 4, 'Neurosynth mask: Emotion regulation'); + +% Here, the regions predictive of success are largely different from those +% that respond to reappraisal demand in past literature. + +%% Part 8: Compare to a standard analysis with no scaling or variable transformation + +% Start over and re-load the images: +[image_obj_untouched] = load_image_set('emotionreg', 'noverbose'); + +% attach success scores to image_obj in image_obj.X +% Just mean-center them so we can interpret the group result +image_obj_untouched.X = scale(success, 1); + +% Do a standard OLS regression and view the results: +out = regress(image_obj_untouched, 'nodisplay'); + +% Treshold at q < .05 FDR and display +t_untouched = threshold(out.t, .05, 'fdr'); + +o2 = montage(t_untouched); +o2 = title_montage(o2, 2, 'No scaling: Behavioral predictor: Reappraisal success (q < .05 FDR)'); +o2 = title_montage(o2, 4, 'No scaling: Intercept: Group average activation'); + +% Since we have no results for the predictor, re-threshold at p < .005 +% uncorrected and display so we can see what's there: + +t_untouched = threshold(t_untouched, .005, 'unc'); + +o2 = montage(t_untouched); +o2 = title_montage(o2, 2, 'No scaling: Behavioral predictor: Reappraisal success (p < .005 uncorrected)'); +o2 = title_montage(o2, 4, 'No scaling: Intercept: Group average activation'); +snapnow + +% select just the image for reappraisal success: +t_untouched = select_one_image(t_untouched, 1); + +%% Make plots comparing the analysis with our choices and a naive analysis directly + +o2 = canlab_results_fmridisplay([], 'multirow', 4); + +t = threshold(t, .05, 'fdr'); +o2 = montage(region(t), o2, 'wh_montages', 1:2, 'colormap'); +o2 = title_montage(o2, 2, 'Reappraisal success with informed choices (q < .05 corrected)'); + +t_untouched = threshold(t_untouched, .05, 'fdr'); +o2 = montage(region(t_untouched), o2, 'wh_montages', 3:4, 'colormap'); +o2 = title_montage(o2, 4, 'Reappraisal success with no scaling (q < .05 corrected)'); + +% Here we'll use "multi-threshold" to display blobs at lower thresholds +% contiguous with those at higher thresholds. We'll use uncorrected p < +% .005 so we can see results in both maps: + +o2 = multi_threshold(t, o2, 'thresh', [.002 .01 .02], 'k', [10 1 1], 'wh_montages', 5:6); +o2 = title_montage(o2, 6, 'Reappraisal success with informed choices (p < .001 one-tailed unc, showing extent to .01 one-tailed)'); + +o2 = multi_threshold(t_untouched, o2, 'thresh', [.002 .01 .02], 'k', [10 1 1], 'wh_montages', 7:8); +o2 = title_montage(o2, 8, 'Reappraisal success with no scaling (p < .001 unc, showing extent to .01 one-tailed)'); + +% One-tailed results are the default in SPM and FSL - so this is what you'd +% get with 0.001 in SPM. + + + +%% Part 9: Localize results using a brain atlas +% The table() command is a simple way to label regions +% In addition, we can visualize each region in table to check the accuracy +% of the labels and the extent of the region. + +% Select the Success regressor map +r = region(t); + +% Autolabel regions and print a table +r = table(r); + +% Make a montage showing each significant region +montage(r, 'colormap', 'regioncenters'); + +%% Part 10: Extract and plot data from regions of interest +% Here, we'll explore data extracted from two kinds of regions: + +% 1 - biased regions based on the significant results. This is good for +% visualizing the distribution of the data in a region and checking +% assumptions, but examining the correlation in these regions is +% "circular", as the correlation was used to select voxels in the first +% place. The observed correlations in these regions will be upwardly +% biased as a result. In general, you cannot estimate the effect size +% (correlation or other metric) in a region or voxel after selecting it +% using criteria that are non-independent of the effect of interest. + +% 2 - unbiased regions from the meta-analysis. These are fair. + + + +%% Extract and plot from (biased) regions of interest based on results +% Here, we'll compare the implications of picking BIASED ROIs (which is +% often done in the literature) vs. picking UNBIASED ROIs. +% +% First, let's visualize the correlation scatterplots in the areas we've +% discovered as related to Success + +% Extract data from all regions +% r(i).dat has the averages for each subject across voxels for region i +r = extract_data(r, image_obj); + +% Select only regions with 3+ voxels +% wh = cat(1, r.numVox) < 3; +% r(wh) = []; + +% Set up the scatterplots +nrows = floor(sqrt(length(r))); +ncols = ceil(length(r) ./ nrows); +create_figure('scatterplot_region', nrows, ncols); + +% Make a loop and plot each region +for i = 1:length(r) + + subplot(nrows, ncols, i); + + % Use this line for non-robust correlations: + %plot_correlation_samefig(r(i).dat, datno16.X); + + % Use this line for robust correlations: + plot_correlation_samefig(r(i).dat, image_obj.X, [], 'k', 0, 1); + + set(gca, 'FontSize', 12); + xlabel('Reappraise - Look Neg brain response'); + ylabel('Reappraisal success'); + + % input('Press a key to continue'); + +end + + +%% Extract and plot data from unbiased regions of interest +% Let's visualize the correlation scatterplots in some meta-analysis +% derived ROIs + +% Select the Neurosynth meta-analysis map +r = region(metaimg); + +% Extract data from all regions +r = extract_data(r, image_obj); + +% Select only regions with 20+ voxels +wh = cat(1, r.numVox) < 20; +r(wh) = []; + +r = table(r); + +% Make a montage showing each significant region +montage(r, 'colormap', 'regioncenters'); + +% Set up the scatterplots +nrows = floor(sqrt(length(r))); +ncols = ceil(length(r) ./ nrows); +create_figure('scatterplot_region', nrows, ncols); + +% Make a loop and plot each region +for i = 1:length(r) + + subplot(nrows, ncols, i); + + % Use this line for non-robust correlations: + %plot_correlation_samefig(r(i).dat, datno16.X); + + % Use this line for robust correlations: + plot_correlation_samefig(r(i).dat, image_obj.X, [], 'k', 0, 1); + + set(gca, 'FontSize', 12); + xlabel('Brain'); + ylabel('Success'); + title(' '); + +end + +%% Part 11: (bonus) Multivariate prediction from unbiased ROI averages +% Predict reappraisal success using brain images +% with 5-fold cross-validation + +% Since we are predicting outcome data from brain images, we have to attach +% data in the .Y field of the fmri_data object. +% .Y is the outcome to be explained, in this case success scores +image_obj.Y = image_obj.X; + +% 5-fold cross validated prediction, stratified on outcome + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', 'cv_lassopcrmatlab', 'nfolds', 5); + +% Plot cross-validated predicted outcomes vs. actual + +% Critical fields in stats output structure: +% stats.Y = actual outcomes +% stats.yfit = cross-val predicted outcomes +% pred_outcome_r: Correlation between .yfit and .Y +% weight_obj: Weight map used for prediction (across all subjects). + +% Continuous outcomes: + +create_figure('scatterplot'); +plot(stats.yfit, stats.Y, 'o') +axis tight +refline +xlabel('Predicted reappraisal success'); +ylabel('Observed reappraisal success'); + +% Though many areas show some significant effects, these are not strong +% enough to obtain a meaningful out-of-sample prediction of Success + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m new file mode 100644 index 00000000..1f27baae --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m @@ -0,0 +1,391 @@ +%% Multivariate prediction of a continuous outcome + +%% +% This walkthrough uses the |predict()| method for fmri_data objects to +% predict a continuous outcome using cross-validated principal component +% regression (PCR / LASSO-PCR). +% +% The |predict()| method can run multiple algorithms with a range of options. +% The main ones used in the CANlab are SVM for two-choice classification and PCR for +% regression. + +%% About the pain dataset +% ------------------------------------------------------------------------- +% +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). Each image is the average over several +% (4-8) trials of heat delivered at a single stimulus intensity, ranging +% from 44.3 - 49.3 degrees C in one-degree increments. Each image is also +% paired with an average reported pain value for that set of trials, rated +% immmediately after heat experience. +% +% This dataset is interesting for mixed-effects and predictive analyses, as +% it has both within-person and between-person sources of variance. +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +%% Load dataset with images and print descriptives +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj +% This is an fmri_data object from the CANlab Core Tools repository for neuroimaging data analysis. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + +%% Collect some variables we need + +% subject_id is useful for cross-validation +% +% this used as a custom holdout set in fmri_data.predict() below will +% implement leave-one-subject-out cross-validation + +subject_id = image_obj.additional_info.subject_id; + +% ratings: reconstruct a subjects x temperatures matrix +% so we can plot it +% +% The command below does this only because we have exactly 6 conditions +% nested within 33 subjects and no missing data. + +ratings = reshape(image_obj.Y, 6, 33)'; + +% temperatures are 44 - 49 (actually 44.3 - 49.3) in order for each person. + +temperatures = image_obj.additional_info.temperatures; + +%% Plot the ratings +% ------------------------------------------------------------------------- + +create_figure('ratings'); +hold on; plot(ratings', '-', 'Color', [.7 .7 .7], 'LineWidth', .5); +lineplot_columns(ratings, 'color', [.7 .3 .3], 'markerfacecolor', [1 .5 0]); +xlabel('Temperature'); +ylabel('Rating'); +set(gca, 'XTickLabel', [44.3:49.3], 'FontSize', 18); + +%% Prediction options and choices +% There are many choices for algorithm and process/parameter tuning +% Here, we describe a constrained set of some of the most common choices +% +% Base model: Linear, whole brain prediction +% +% Outcome distribution: +% Continuous: Regression (LASSO-PCR) Categorical/2 choice: SVM, logistic regression +% +% Model Options: +% - Feature selection (part of the brain) +% - Feature integration (new features) +% +% Testing options: +% - Input data (what level) +% - Cross-validation +% - What testing metric (classification accuracy? RMSE?) and at what level +% (single-trial, condition averages, one-map-per-subject) +% +% Inference options: +% - Component maps, voxel-wise maps +% - Thresholding (Bootstrapping, Permutation) +% - Global null hypothesis testing / sanity checks (Permutation) + +%% Run the Base model +% +% *relevant functions:* +% predict method (fmri_data) +% predict_test_suite method (fmri_data) +% + +%% +% *Define the holdout set for cross-validation* +% We want to define custom holdout set. If we use subject_id, which is a vector of +% integers with a unique integer per subject, then we are doing +% leave-one-subject-out cross-validation. +% +% Let's build five-fold cross-validation set that leaves out ALL the images +% from a subject together. That way, we are always predicting out-of-sample +% (new individuals). If not, dependence across images from the same +% subjects may invalidate the estimate of predictive accuracy. + +holdout_set = zeros(size(subject_id)); % holdout set membership for each image +n_subjects = length(unique(subject_id)); +C = cvpartition(n_subjects, 'KFold', 5); +for i = 1:5 + teidx = test(C, i); % which subjects to leave out + imgidx = ismember(subject_id, find(teidx)); % all images for these subjects + holdout_set(imgidx) = i; +end + +%% +% *Run the prediction* + +algoname = 'cv_lassopcr'; % cross-validated penalized regression. Predict pain ratings + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +%% Plot cross-validated predicted vs. actual outcomes +% +% Critical fields in stats output structure: +% stats.Y = actual outcomes +% stats.yfit = cross-val predicted outcomes +% pred_outcome_r: Correlation between .yfit and .Y +% weight_obj: Weight map used for prediction (across all subjects). + +% Continuous outcomes: + +create_figure('scatterplots', 1, 2); +plot(stats.yfit, stats.Y, 'o') +axis tight +refline +xlabel('Predicted pain'); +ylabel('Observed pain'); + +% Consider that each subject has their own series of multiple observations. +% Separate subjects into multiple cells, which will be plotted as separate +% individual lines: + +n = max(subject_id); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +subplot(1, 2, 2); +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +% Yfit is a cell array, one cell per subject, with fitted values for that subject + +%% Summarize within-person classification accuracy +% +% Classification can turn a continuous outcome, such as the magnitude of +% activity in an integrated model/pattern, into a binary choice (Type A or +% B, Yes/No, Pain/No pain). Classification accuracy is an easy-to-interpret +% metric, which makes it appealing to report. +% +% But though it seems easy to understand, it is perhaps tricker to +% interpret than we might think. +% +% Classification accuracy is not one metric, but a family of them. The accuracy +% you obtain depends on what type of classification you're doing, and what the +% units of analysis are. For example, is your model based on a single person +% (individualized for that person), or based on other participants' data? +% Are you classifying images corresponding to a single trial, or an average across a group of trials? +% Do you have knowledge about which conditions/images belong to the same person, or not? +% Here, we will: +% (a) use a model based on other people's data (cross-validated across participants) +% (b) classify images summarizing a group of trials, and +% (c) assume we have 2 images belonging to the same person, and we classify +% which is which +% +% On this last point, there are two basic senses of classification, which assume different +% levels of background knowledge: +% - In _single-interval classification_, we do not know which images belong to which participants, +% and we are classifying a single image as Type A or B based on whether +% the response in the model is above or below a threshold. +% - In _forced-choice_ classification, we have two images that we know come +% from the same person, and we know one is Type A and one is Type B...or, +% equivalently, that one is "more A" than the other (e.g., more painful). +% The classifier tries to pick out which is the best answer. +% +% Forced-choice classification is an inherently easier problem, with higher +% resulting accuracy, because (a) each participant serves as their own control +% (many sources of noise are canceled out), and (b) you assume more +% background knowledge, conferring an advantage. +% +% In this case, we want to summarize each subject with two images: One high +% pain, and one low pain. We'll use stimulus intensity as the variable to +% classify here, and ask if we can predict which intensity is highest given +% pairs of images one degree apart for each participant. + +% If the stim intensity increases by 1 degree, how often do we correctly +% predict an increase? +diffs = cellfun(@diff, Yfit, 'UniformOutput', false); % successive increases in predicted values with increases in stim intensity + +% Subjects (rows) x comparisons (cols), each comparison 1 degree apart +diffs = cat(2, diffs{:})'; + +acc_per_subject = sum(diffs > 0, 2) ./ size(diffs, 2); +acc_per_comparison = (sum(diffs > 0) ./ size(diffs, 1))'; +cohens_d_per_comparison = (mean(diffs) ./ std(diffs))'; + +fprintf('Mean acc is : %3.2f%% across all successive 1 degree increments\n', 100*mean(acc_per_subject)); + +% Create and display a table: + +comparison = {'45-44' '46-45' '47-46' '48-47' '49-48'}'; +results_table = table(comparison, acc_per_comparison, cohens_d_per_comparison); +disp(results_table) + +% 45-44 degrees, 46-45, 47-46, etc. + +%% Visualize the classifier weight map +% +% The weight map is stored in |stats.weight_obj|, as a statistic_image +% class object. If we have requested bootstrapping, the image will be +% thresholded (0.05 uncorrected by default). Otherwise, it will be unthresholded, +% and we'll see weights everywhere within the analysis mask. + +w = stats.weight_obj; % This is an image object that we can view, etc. + +orthviews(w) + +create_figure('montage'); +o2 = montage(w); +o2 = title_montage(o2, 5, 'Predictive model weights'); + + +%% Bootstrap weights: Get most reliable weights and p-values for voxels +% +% Here is an exercise: +% Re-run the predictive model, adding a 'bootstrap' flag. +% For a final analysis, at least 5K bootstrap samples is a good idea. +% For now, try it with just 1,000 bootstrap samples (and be prepared to +% wait a bit). What does the resulting weight map look like? +% +% (This is not run in the example; see |help fmri_data.predict| for how to +% add the bootstrap option.) + +%% Try normalizing/scaling data +% +% One of the biggest sources of noise can be whole-image (or +% large-spatial-scale shifts) in image values across participants. These can +% be apparent even in contrast or "subtraction" images where various +% sources of noise and artifacts are supposed to have been "subtracted +% out". This may result from different chance correlations between task +% regressors and physiological noise or artifacts/outliers across +% individuals. +% +% Rescaling the data can remove whole-brain signal. This should be used +% with caution, as you're also removing signal from the images, and the +% resulting weight maps should be interpreted as effects *relative* to the mean +% signal across the image. i.e., a task may be associated with whole-brain +% increases, and *relative* decreases in a local area that appear only +% after the global increase has been accounted for. +% +% Another reason to rescale images is to ask whether the *pattern* of +% activity across the brain (or within a local area) is predictive of task +% state, after removing the overall intensity of activation. This relates +% to pattern information. +% +% Here, we'll rescale the images and ask if this improves our +% cross-validated predictions or not. If it does, the global signal we've +% removed is likely more noise than signal. If not, the global signal may +% contain information about pain. + +% z-score every image, removing the mean and dividing by the std. +% Remove some scaling noise, but also some signal... + image_obj = rescale(image_obj, 'zscoreimages'); + +% re-run prediction and plot + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +create_figure('scatterplots'); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +fprintf('Overall prediction-outcome correlation is %3.2f\n', stats.pred_outcome_r) + +figure; o2 = montage(stats.weight_obj); +o2 = title_montage(o2, 5, 'Predictive model weights'); + +%% Try selecting an a priori 'network of interest' + +% re-load +load('bmrk3_6levels_pain_dataset.mat', 'image_obj') + +% load an a priori pain-related mask and apply it + +mask = fmri_data(which('pain_2s_z_val_FDR_05.img.gz')); % in Neuroimaging_Pattern_Masks repository + +image_obj = apply_mask(image_obj, mask); + +% show mean data with mask +m = mean(image_obj); +figure; o2 = montage(m); +o2 = title_montage(o2, 5, 'Mean image data, masked'); + + +%% re-run prediction and plot + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +create_figure('scatterplots'); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +fprintf('Overall prediction-outcome correlation is %3.2f\n', stats.pred_outcome_r) + +figure; o2 = montage(stats.weight_obj); +o2 = title_montage(o2, 5, 'Predictive model weights'); + +%% Other ideas +% There are many ways to build from this basic analysis. + +% Data +% ------------ +% Examine/clean outliers +% Transformations of data: Scaling/centering to remove global mean, +% ventricle mean +% Identify unaccounted sources of variation/subgroups (experimenter sex, +% etc.) + +% Features +% ------------ +% A priori regions / meta-analysis +% Feature selection in cross-val loop: Univariate, recursive feature elim + From 1b0c93507f1c740bde06b323353560d24b6cce9c Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Fri, 19 Jun 2026 19:00:51 -0400 Subject: [PATCH 8/8] Walkthrough CI: handle missing Image Processing Toolbox First CI run of the hardened nightly surfaced 3 failures, all the same environment gap: read_nifti_volume falls back to niftiinfo (Image Processing Toolbox), which the runner did not provision (only Statistics + Signal were installed), so plot(obj) failed with Undefined function 'niftiinfo'. - Provision Image_Processing_Toolbox in tests-walkthroughs.yml so plot() and NIfTI I/O actually run on CI. - Safety net: canlab_classify_environment_error now buckets an UndefinedFunction error for a known optional-toolbox function (niftiinfo, niftiread, niftiwrite, cfg_getfile) as 'capability' -> skipped, not failed, so a missing toolbox can never spuriously redden the nightly. Genuine missing-function bugs still classify as 'genuine'. Co-Authored-By: Claude Opus 4.8 (1M context) --- .github/workflows/tests-walkthroughs.yml | 1 + .../canlab_classify_environment_error.m | 21 +++++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/.github/workflows/tests-walkthroughs.yml b/.github/workflows/tests-walkthroughs.yml index e1617ef2..f7664896 100644 --- a/.github/workflows/tests-walkthroughs.yml +++ b/.github/workflows/tests-walkthroughs.yml @@ -49,6 +49,7 @@ jobs: products: | Statistics_and_Machine_Learning_Toolbox Signal_Processing_Toolbox + Image_Processing_Toolbox cache: true - name: Run walkthrough tests diff --git a/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m index a6e3e587..cb750be2 100644 --- a/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m +++ b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m @@ -36,6 +36,9 @@ % (orthviews, surface rendering, etc.). % 'input' - code tried to prompt for interactive user input, % unavailable in batch/CI. +% 'capability' - an UndefinedFunction error for a known optional +% MATLAB toolbox function (e.g. niftiinfo) that is not +% provisioned on this runner. % 'data' - an optional data file (signature, atlas, feature set) % is not on the CI path. % 'cascade' - an undefined-variable/field error following an earlier @@ -104,6 +107,24 @@ return end +% --------------------------------------------------------------------- +% Missing optional MATLAB toolbox. An UndefinedFunction error for a known +% MathWorks toolbox function (e.g. niftiinfo/niftiread from the Image +% Processing Toolbox) means that toolbox is not provisioned on this runner, +% not that the code is broken. Keep this list to functions that are *only* +% ever provided by a toolbox, so we never mask a genuine missing-CanlabCore +% function bug. +% --------------------------------------------------------------------- +optional_toolbox_fns = {'niftiinfo', 'niftiread', 'niftiwrite', ... + 'cfg_getfile'}; +if strcmp(id, 'MATLAB:UndefinedFunction') + undef = regexp(ME.message, "Undefined function '([^']+)'", 'tokens', 'once'); + if ~isempty(undef) && any(strcmpi(undef{1}, optional_toolbox_fns)) + category = 'capability'; + return + end +end + % --------------------------------------------------------------------- % Missing optional data files (NPS+ signatures, Neurosynth feature set, % Bianciardi atlas sources, etc.). load_image_set / annotate_* abort with a