diff --git a/CLAUDE.md b/CLAUDE.md index 3bfb7dd0..5b5cfdcd 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. Its property names mirror the fields of the results structure `fmri_data.regress` builds (the variable `out` in that method): the design spec (wrapping an `fmri_glm_design_matrix` in `.design`, or a direct `.X`), the fitted result maps (`.betas`/`.t`/`.contrast_estimates`/`.contrast_t` as `statistic_image`; `.df`/`.sigma`/`.residuals` as `fmri_data`), and three nested option/diagnostic structs `.input_parameters`, `.input_image_metadata`, `.diagnostics` (VIF, contrast VIF, leverage, condition number, collinearity). The historical out-struct field names are available as read/write aliases (`.b`→`.betas`, `.con_t`→`.contrast_t`, `.contrast_images`→`.contrast_estimates`, `.resid`→`.residuals`, `.variable_names`→`.regressor_names`, `.C`→`.contrasts`). Two workflows: (1) call `fmri_data.regress` directly — it now **returns a glm_map** (not a struct) — and query the object; (2) build an estimator: `g = glm_map(...)`, `add_contrasts`, `run_diagnostics`, `g = fit(g, fmri_data_obj)`, then `threshold`/`table`/`montage`. (The results live in the `.diagnostics` property; the method that computes them is `run_diagnostics` to avoid a property/method name clash.) `glm_map(out_struct)` re-casts a regress-style struct into an object. Supports 1st-level event designs (`build_design`, `import_SPM`) and 2nd-level/group designs; `fmri_data.regress` is the compute engine. Walkthrough: `docs/workflows/regression_with_glm_map.m`. ### MATLAB `@class` directories diff --git a/CanlabCore/@fmri_data/regress.m b/CanlabCore/@fmri_data/regress.m index 9c7302ad..482c33eb 100644 --- a/CanlabCore/@fmri_data/regress.m +++ b/CanlabCore/@fmri_data/regress.m @@ -99,10 +99,13 @@ % :Outputs: % % **regression_results:** -% A structure containing stats_img and fmri_data objects. -% In addition to the main outputs below, the regression_results structure also has -% fields for input_parameters, the design matrix (X), variable -% names, and warnings. +% A glm_map object containing stats_img and fmri_data objects. +% (Historically this was a plain struct; it is now re-cast as a +% glm_map whose property names match the field names documented here, +% so out.b, out.t, out.contrast_images, out.con_t, etc. all still +% work. See help glm_map.) In addition to the main outputs below, it +% also carries input_parameters, input_image_metadata, the design +% matrix (X), variable names, diagnostics, and warnings. % % **regression_results.b:** % stats_img object of beta values estimated from regression @@ -121,7 +124,8 @@ % % **regression_results.diagnostics:*** % A structure containing VIFs and leverage values for the design matrix -% regression_results.diagnostics.Variance_inflation_factors = VIFs +% regression_results.diagnostics.Variance_inflation_factors = VIFs (from VIF.m) +% regression_results.diagnostics.Contrast_variance_inflation_factors = contrast VIFs (cVIF.m), if C provided % regression_results.diagnostics.Leverages = leverage values % % :Examples: @@ -563,8 +567,17 @@ % Variance inflation % --------------------------------------------------------------------- - -vifs = getvif(X); +% Use canonical VIF.m (getvif is deprecated). When a contrast matrix C is +% present, VIF also returns contrast VIFs (cVIF) for the contrasts. C is +% stored column-wise [regressors x contrasts]; VIF/cVIF expect one contrast +% per row, so pass C'. + +contrast_vifs = []; +if ~isempty(C) && size(C, 1) == size(X, 2) + [vifs, contrast_vifs] = VIF(X, C'); +else + vifs = VIF(X); +end if any(vifs > 4) @@ -572,6 +585,10 @@ end +if ~isempty(contrast_vifs) && any(contrast_vifs > 4) + mywarnings{end+1} = 'Warning!!! Some contrasts have contrast variance inflation factors (cVIF) > 4. Check regression_results.diagnostics.Contrast_variance_inflation_factors'; +end + % Leverages % --------------------------------------------------------------------- @@ -1067,7 +1084,8 @@ end -out.diagnostics = struct('Variance_inflation_factors', vifs, 'Leverages', leverages); +out.diagnostics = struct('Variance_inflation_factors', vifs, ... + 'Contrast_variance_inflation_factors', contrast_vifs, 'Leverages', leverages); out.warnings = mywarnings; % Create objects @@ -1199,6 +1217,25 @@ warning on end +% --------------------------------------------------------------------- +% Re-cast the results structure as a glm_map object +% --------------------------------------------------------------------- +% The fields built above (b, t, df, sigma, contrast_images, con_t, +% input_parameters, input_image_metadata, diagnostics, ...) map onto glm_map +% properties of the same names. The object preserves struct-style field +% access (out.b, out.t, out.con_t, ...) via Dependent aliases, so existing +% callers continue to work. Guard the conversion so that any non-standard +% output path (or a missing glm_map on the path) falls back to the struct. +if isstruct(out) && isfield(out, 'b') && isfield(out, 't') && exist('glm_map', 'class') == 8 + try + out = glm_map(out); + out = validate_object(out); % ensure all nested-struct fields are present + catch ME + warning('fmri_data:regress:glm_mapCastFailed', ... + 'Returning results as a struct; could not cast to glm_map: %s', ME.message); + end +end + % --------------------------------------------------------------------- % Subfunctions % --------------------------------------------------------------------- diff --git a/CanlabCore/@fmri_glm_design_matrix/add.m b/CanlabCore/@fmri_glm_design_matrix/add.m index c8c9eb8a..f5133bed 100644 --- a/CanlabCore/@fmri_glm_design_matrix/add.m +++ b/CanlabCore/@fmri_glm_design_matrix/add.m @@ -97,18 +97,31 @@ isothermethod(i + 1) = false; case 'condition_names' - - % special function for checking SPM-style format rules - [names, nsess, nconds] = check_req_and_length(obj, varargin{i+1}); - - indx = 1; + + % Condition names are assumed to be the same for all sessions. + % Accept either nconds names (shared across sessions) or + % nsess*nconds names (session-major). Use the actual number of + % conditions from the onset structures when available. + names = varargin{i + 1}; + if ~iscell(names), names = {names}; end + + nsess = length(obj.nscan); + nconds = length(obj.Sess(1).U); + if nconds == 0 + nconds = length(names) ./ nsess; % legacy fallback (no onsets yet) + end + + shared = (numel(names) == nconds); for ss = 1:nsess for condname = 1:nconds - obj.Sess(ss).U(condname).name = names{indx}; - indx = indx + 1; + if shared + obj.Sess(ss).U(condname).name = names{condname}; + else + obj.Sess(ss).U(condname).name = names{(ss - 1) * nconds + condname}; + end end end - + isothermethod(i + 1) = false; case 'SPM' diff --git a/CanlabCore/@fmri_glm_design_matrix/build.m b/CanlabCore/@fmri_glm_design_matrix/build.m index c78e3931..9e3646c7 100644 --- a/CanlabCore/@fmri_glm_design_matrix/build.m +++ b/CanlabCore/@fmri_glm_design_matrix/build.m @@ -128,7 +128,10 @@ end for i = 1:nconds - obj.xBF(i).name = sprintf('%s for Condition %3.0f', obj.xBF(i).name, i); + % Idempotent: strip any prior " for Condition N" suffix so repeated + % build() calls do not keep appending to the basis-set name. + basename = regexprep(obj.xBF(i).name, '\s*for Condition\s+\d+\s*$', ''); + obj.xBF(i).name = sprintf('%s for Condition %3.0f', basename, i); end end diff --git a/CanlabCore/@fmri_glm_design_matrix/build_single_trial.m b/CanlabCore/@fmri_glm_design_matrix/build_single_trial.m index e1724ff5..a4ec6fbf 100644 --- a/CanlabCore/@fmri_glm_design_matrix/build_single_trial.m +++ b/CanlabCore/@fmri_glm_design_matrix/build_single_trial.m @@ -45,6 +45,8 @@ [sess_delta, sess_conditions, sess_X, sess_C, sess_B] = deal(cell(1, nsess)); +res = 16; % high-resolution samples per second (TR-independent) + for s = 1:nsess [ons, name] = deal(cell(1, nconds)); @@ -52,50 +54,48 @@ [ons{:}] = deal(obj.Sess(s).U(:).ons); %[name{:}] = deal(obj.Sess(s).U(:).name); -% make sure onsets are in TRs, not secs -switch obj.xBF(1).UNITS - case 'secs' - % onsets are in sec, convert - for i = 1:nconds - ons{i} = ons{i} ./ TR; - end - - case {'tr', 'TR', 'trs'} - % ok, do nothing +% Convert onsets to seconds (we build at high resolution, downsample to TR) +to_sec = 1; +if any(strcmpi(obj.xBF(1).UNITS, {'scans', 'tr', 'trs'})), to_sec = TR; end +for i = 1:nconds + ons{i} = ons{i}(:) * to_sec; end % ---------------------------------------------- % Build predictors for each session % ---------------------------------------------- -delta = onsets2delta(ons, obj.nscan(s)); - -ns = size(delta, 1); +nscan = obj.nscan(s); +len_hires = round(nscan * TR * res); [trialdelta, cond_assignment] = deal(cell(1, nconds)); % ---------------------------------------------- -% For each condition, parse into separate columns -% for each onset (single-trial) +% For each condition, parse into separate columns for each onset +% (single-trial), at high resolution; convolve and downsample to TR. % ---------------------------------------------- for j = 1:nconds - - wh = find(delta(:, j)); - - trialdelta{j} = false(ns, length(wh)); - - for k = 1:length(wh) - trialdelta{j}(wh(k), k) = true; + + os = ons{j}; + ntrials = numel(os); + + trialdelta{j} = zeros(len_hires, ntrials); % one hi-res column per trial + for k = 1:ntrials + idx = round(os(k) * res) + 1; + if idx >= 1 && idx <= len_hires, trialdelta{j}(idx, k) = 1; end end - - cond_assignment{j} = ones(length(wh), 1); - + + cond_assignment{j} = ones(ntrials, 1); + + % custom per-condition HRF, sampled at high resolution (res samples/s) bf = obj.xBF(j); xvals = 0:TR:bf.length - TR; - hrf = interp1(xvals, inputhrf{j}(1:length(xvals)), 1:bf.dt:bf.length, 'spline', 'extrap'); - - condX{j} = getPredictors(trialdelta{j}, hrf, 16); - + hrf = interp1(xvals, inputhrf{j}(1:length(xvals)), 0:1/res:bf.length, 'spline', 'extrap'); + hrf = hrf(:); + + Xi = getPredictors(trialdelta{j}, hrf, 'dsrate', res * TR, 'force_delta'); + condX{j} = Xi(1:nscan, :); + end % conditions sess_delta{s} = cat(2, trialdelta{:}); % not needed? @@ -116,13 +116,15 @@ % Put the pieces together % ---------------------------------------------- -obj.xX.X = [blkdiag(sess_X{:}) blkdiag(sess_C{:}) blkdiag(sess_B{:})]; +% Index xX(1) explicitly: obj.xX starts as an empty (0x0) struct, and a +% dot-assignment to an empty struct array is illegal. +obj.xX(1).X = [blkdiag(sess_X{:}) blkdiag(sess_C{:}) blkdiag(sess_B{:})]; -obj.xX.cond_assignments = cat(1, sess_conditions{:}); +obj.xX(1).cond_assignments = cat(1, sess_conditions{:}); -obj.xX.iH = find(any(obj.xX.cond_assignments, 2)); -obj.xX.iC = find(~any(obj.xX.cond_assignments, 2)); -obj.xX.iB = size(obj.xX.X, 2) - nsess : size(obj.xX.X, 2); +obj.xX(1).iH = find(any(obj.xX(1).cond_assignments, 2)); +obj.xX(1).iC = find(~any(obj.xX(1).cond_assignments, 2)); +obj.xX(1).iB = size(obj.xX(1).X, 2) - nsess + 1 : size(obj.xX(1).X, 2); end % function @@ -174,7 +176,10 @@ end for i = 1:nconds - obj.xBF(i).name = sprintf('%s for Condition %3.0f', obj.xBF(i).name, i); + % Idempotent: strip any prior " for Condition N" suffix so repeated + % build calls do not keep appending to the basis-set name. + basename = regexprep(obj.xBF(i).name, '\s*for Condition\s+\d+\s*$', ''); + obj.xBF(i).name = sprintf('%s for Condition %3.0f', basename, i); end end diff --git a/CanlabCore/@fmri_glm_design_matrix/create_orthogonal_contrast_set.m b/CanlabCore/@fmri_glm_design_matrix/create_orthogonal_contrast_set.m new file mode 100644 index 00000000..662ed759 --- /dev/null +++ b/CanlabCore/@fmri_glm_design_matrix/create_orthogonal_contrast_set.m @@ -0,0 +1,89 @@ +function obj = create_orthogonal_contrast_set(obj, varargin) +% Assign an orthogonal set of contrasts spanning the regressors of interest. +% +% Uses create_orthogonal_contrast_set.m to build a set of mutually orthogonal, +% zero-sum contrasts spanning the space of differences among the regressors of +% interest (the H partition, obj.xX.iH), and stores them in obj.xX.contrasts +% (one row per contrast, one column per design regressor) with placeholder +% names in obj.xX.contrast_names. Nuisance covariates and session constants +% get zero weight. The model must be built first (so obj.xX.iH is defined). +% +% :Usage: +% :: +% +% obj = create_orthogonal_contrast_set(obj) +% obj = create_orthogonal_contrast_set(obj, 'names', {'c1','c2',...}) +% +% :Inputs: +% +% **obj:** +% A built fmri_glm_design_matrix object (obj.xX.iH non-empty). +% +% :Optional Inputs: +% +% **'names', {cellstr}:** contrast names (default 'OrthC1', ...). +% +% :Outputs: +% +% **obj:** +% with obj.xX.contrasts [n_contrasts x n_regressors] and +% obj.xX.contrast_names set. +% +% :See also: +% - create_orthogonal_contrast_set (function), build, glm_map.create_orthogonal_contrast_set +% +% .. +% 2026 - Initial implementation. +% .. + +names = {}; +wh = find(strcmpi(varargin, 'names')); +if ~isempty(wh), names = varargin{wh(1) + 1}; end + +% Regressors of interest = H partition. Requires a built design. +whI = []; +if isstruct(obj.xX) && isscalar(obj.xX) && isfield(obj.xX, 'iH') && ~isempty(obj.xX.X) + whI = obj.xX.iH(:)'; +end + +if isempty(whI) + error('fmri_glm_design_matrix:NoInterestRegressors', '%s', local_interest_guidance()); +end +if numel(whI) < 2 + error('fmri_glm_design_matrix:TooFewInterest', ... + 'Orthogonal contrasts require at least 2 regressors of interest (obj.xX.iH); found %d.', numel(whI)); +end + +nreg = size(obj.xX.X, 2); +Cint = create_orthogonal_contrast_set(numel(whI)); % [(n-1) x n], rows = contrasts +ncon = size(Cint, 1); + +C = zeros(ncon, nreg); +C(:, whI) = Cint; + +if isempty(names) + names = arrayfun(@(i) sprintf('OrthC%d', i), 1:ncon, 'UniformOutput', false); +elseif numel(names) ~= ncon + error('fmri_glm_design_matrix:ContrastNames', 'Number of names (%d) must match number of contrasts (%d).', numel(names), ncon); +end + +obj.xX(1).contrasts = C; +obj.xX(1).contrast_names = names(:)'; + +if ~iscell(obj.history), obj.history = {}; end +obj.history{end + 1} = sprintf('create_orthogonal_contrast_set: %d contrast(s) over %d interest regressors', ncon, numel(whI)); + +end % create_orthogonal_contrast_set + + +% ========================================================================= +function s = local_interest_guidance() +s = sprintf([ ... + 'No regressors of interest are defined (obj.xX.iH is empty).\n' ... + 'Regressors of interest are the event/task regressors (the H partition).\n' ... + 'Add event onsets and build the design first, e.g.:\n' ... + ' d = fmri_glm_design_matrix(TR, ''nscan'', nscan, ''units'', ''secs'', ...\n' ... + ' ''onsets'', onsets, ''condition_names'', names);\n' ... + ' d = build(d);\n' ... + 'After build, the event regressors are placed in obj.xX.iH.']); +end diff --git a/CanlabCore/@fmri_glm_design_matrix/fmri_glm_design_matrix.m b/CanlabCore/@fmri_glm_design_matrix/fmri_glm_design_matrix.m index dea58d65..b7ecd8da 100644 --- a/CanlabCore/@fmri_glm_design_matrix/fmri_glm_design_matrix.m +++ b/CanlabCore/@fmri_glm_design_matrix/fmri_glm_design_matrix.m @@ -24,7 +24,9 @@ % % my_model = fmri_glm_design_matrix(2); % creates an empty fmri_glm_design_matrix object with a TR of 2 (which is used to -% create a default b-spline basis set for the HRF. +% create a default basis set: the canonical SPM HRF, one basis function per +% condition). To use a multi-basis set (b-splines, HRF + derivatives, ...), +% pass a custom xBF or call replace_basis_set after construction. % % my_model = fmri_glm_design_matrix(2, 'nscan', [198 198 198]); % creates an empty structure but assigns data to the field nscans, number @@ -167,11 +169,22 @@ obj.Sess = []; % : [1xs struct] - Session structure array obj.xX = []; - % Default basis set + % Default basis set: the canonical SPM HRF (a single basis + % function per condition), so each condition becomes one + % HRF-convolved regressor. This matches SPM's default and is the + % most common choice. To use a multi-basis set (e.g. b-splines or + % HRF + derivatives) instead, pass a custom xBF (e.g. from + % fmri_spline_basis or spm_get_bf) or call replace_basis_set after + % construction. % ---------------------------------------------------------------------- obj.xY.RT = obj.TR; % : - repetition time {seconds) - - obj.xBF = fmri_spline_basis(obj.TR, 0); + + % Sample the basis set at high resolution (16 samples/second, + % TR-independent). get_session_X builds a 16 Hz onset/epoch delta, + % convolves with this basis, and downsamples to TR -- the standard + % microtime pipeline, robust to fractional TRs. + obj.xBF = struct('dt', 1/16, 'name', 'hrf', 'length', 32, 'order', 1); + obj.xBF = spm_get_bf(obj.xBF); % adds obj.xBF.bf (canonical HRF, single column) % obj.xBF.name: - name of basis set % obj.xBF.length: - support of basis set {seconds} % obj.xBF.order: - order of basis set diff --git a/CanlabCore/@fmri_glm_design_matrix/get_session_X.m b/CanlabCore/@fmri_glm_design_matrix/get_session_X.m index 8108f9e7..474a8c5e 100644 --- a/CanlabCore/@fmri_glm_design_matrix/get_session_X.m +++ b/CanlabCore/@fmri_glm_design_matrix/get_session_X.m @@ -1,164 +1,225 @@ function [Xs, delta, C, B, names] = get_session_X(obj, s) -% Get design matrix (predictors) for one session of fmri_model object, using -% basis functions defined in the object and onsets for one session (s). +% Get design matrix (predictors) for one session of an fmri_glm_design_matrix +% object, using the basis functions defined in the object and the onsets for +% one session (s). % % :Usage: % :: % % [Xs, delta, C, B, names] = get_session_X(obj, session number) % +% Method: a high-resolution (microtime) pipeline. Onsets and event durations +% are sampled into a high-resolution indicator at RES samples per second +% (RES = 16, TR-independent), convolved with the per-condition basis set +% (resampled to the same resolution), and downsampled to the TR with +% getPredictors. This matches onsets2fmridesign and is robust to fractional +% TRs. Each condition may have its own basis set (one or more basis +% functions), and parametric modulators are supported. +% +% :Outputs: +% **Xs:** [nscan x k] of-interest predictors (conditions x basis fns, then PM regressors) +% **delta:** [nscan x nconds] TR-resolution onset indicator +% **C:** user-specified covariates for this session ([] if none) +% **B:** [nscan x 1] baseline (intercept) for this session +% **names:** 1 x k cell of regressor names +% +% :See also: getPredictors, onsets2fmridesign, fmri_glm_design_matrix.build % .. -% Define sessions and number of conditions +% High-resolution microtime sampling rate (samples per second). TR- +% independent, matching onsets2fmridesign (res = 16). % .. +res = 16; nsess = length(obj.Sess); - if s > nsess, error('Session %3.0f does not exist', s); end - TR = obj.xY.RT; +nscan = obj.nscan(s); -nconds = length(obj.Sess(s).U); - -[ons, name] = deal(cell(1, nconds)); +if isempty(nscan) || nscan < 1 + error('fmri_glm_design_matrix:get_session_X', ... + 'obj.nscan(%d) is empty or invalid; set the number of scans for each session.', s); +end -[ons{:}] = deal(obj.Sess(s).U(:).ons); -[name{:}] = deal(obj.Sess(s).U(:).name); +nconds = length(obj.Sess(s).U); -% make sure onsets are in TRs, not secs -switch obj.xBF(1).UNITS - case 'secs' - % onsets are in sec, convert - for i = 1:nconds - ons{i} = ons{i} ./ TR; - end - - case {'tr', 'TR', 'trs'} - % ok, do nothing +% Units of onsets/durations: 'secs' (default) or 'scans'/'TR' +units = ''; +if ~isempty(obj.xBF) && isfield(obj.xBF(1), 'UNITS') && ~isempty(obj.xBF(1).UNITS) + units = obj.xBF(1).UNITS; end +if any(strcmpi(units, {'scans', 'tr', 'trs'})), to_sec = TR; else, to_sec = 1; end + +% Session length in high-resolution samples. nscan*TR*res is an integer +% multiple of res*TR, so getPredictors downsamples cleanly to nscan rows. +len_hires = round(nscan * TR * res); % ---------------------------------------------- -% Predictors +% Of-interest predictors: one block per condition (HRF-convolved) % ---------------------------------------------- +Xs = cell(1, nconds); +delta = zeros(nscan, nconds); +condname = cell(1, nconds); -delta = onsets2delta(ons, obj.nscan(s)); +for i = 1:nconds -% Of-interest part of design matrix -% time res is defined as TR / 16, so build and downsample by 16 to TR -% Allow for different basis sets for each condition. + U = obj.Sess(s).U(i); -Xs = cell(1, nconds); + condname{i} = local_getfield(U, 'name', sprintf('R%d', i)); + if iscell(condname{i}) + if isempty(condname{i}), condname{i} = sprintf('R%d', i); else, condname{i} = condname{i}{1}; end + end -for i = 1:nconds - - bf = obj.xBF(i).bf; - - Xs{i} = getPredictors(delta(:, i), bf, 16); -end + ons_sec = U.ons(:) * to_sec; -Xs = cat(2, Xs{:}); + % Durations (sec): empty/scalar -> applied to all events; else per-event + dur = local_getfield(U, 'dur', []); + if isempty(dur), dur = 0; end + dur = dur(:); + if isscalar(dur), dur = repmat(dur, numel(ons_sec), 1); end + dur_sec = dur * to_sec; -% covariate part of design matrix -C = []; -if ~isempty(obj.Sess(s).C) && isfield(obj.Sess(s).C, 'C') - C = obj.Sess(s).C.C; -end + % High-resolution onset/epoch indicator for this condition + dhr = local_hires_indicator(ons_sec, dur_sec, res, len_hires); -% baseline -% ---------------------------------------------- + % Basis set for this condition, resampled to res samples/second + bf = local_bf_at_res(obj.xBF(min(i, numel(obj.xBF))), res); + + Xi = getPredictors(dhr, bf, 'dsrate', res * TR, 'force_delta'); + Xs{i} = local_fit_rows(Xi, nscan); -B = ones(obj.nscan(s), 1); + % TR-resolution indicator (for reference/output) + tr_idx = round(ons_sec / TR) + 1; + tr_idx = tr_idx(tr_idx >= 1 & tr_idx <= nscan); + delta(tr_idx, i) = 1; +end + +Xs = cat(2, Xs{:}); % ---------------------------------------------- -% modulators +% Parametric modulators (appended after the of-interest block) % ---------------------------------------------- - Xs_pm = cell(1, nconds); for i = 1:nconds - is_pm = isfield(obj.Sess(s).U(i), 'P') && ~isempty(obj.Sess(s).U(i).P) && isfield(obj.Sess(s).U(i).P, 'P') && ~isempty(obj.Sess(s).U(i).P.P); - - if is_pm - - pm_vals = {obj.Sess(s).U(i).P.P}; - - model = onsets2parametric_mod_X(ons(i), pm_vals, obj.nscan(s), obj.xBF(i).bf, 16); - - Xs_pm{i} = model; - end + U = obj.Sess(s).U(i); + if ~local_has_pm(U), continue, end + + ons_sec = U.ons(:) * to_sec; + dur = local_getfield(U, 'dur', []); if isempty(dur), dur = 0; end + dur = dur(:); if isscalar(dur), dur = repmat(dur, numel(ons_sec), 1); end + dhr = local_hires_indicator(ons_sec, dur * to_sec, res, len_hires); + + bf = local_bf_at_res(obj.xBF(min(i, numel(obj.xBF))), res); + pm_vals = U.P.P(:); + Xpm = getPredictors(dhr, bf, 'dsrate', res * TR, 'force_delta', ... + 'parametric_singleregressor', {pm_vals}); + Xs_pm{i} = local_fit_rows(Xpm, nscan); end -Xs_pm = cat(2, Xs_pm{:}); - -% add to Xs -Xs = [Xs Xs_pm]; +Xs = [Xs cat(2, Xs_pm{:})]; % ---------------------------------------------- -% get names for each BF -% Allow variable number of basis fcns for each condition +% Covariates and baseline % ---------------------------------------------- +C = []; +if ~isempty(obj.Sess(s).C) && isfield(obj.Sess(s).C, 'C') + C = obj.Sess(s).C.C; +end -names = {}; +B = ones(nscan, 1); +% ---------------------------------------------- +% Names: one per basis function per condition, then PM names +% ---------------------------------------------- +names = {}; for i = 1:nconds - - nbf = size(obj.xBF(i).bf, 2); - + nbf = size(local_bf_at_res(obj.xBF(min(i, numel(obj.xBF))), res), 2); for j = 1:nbf - myname = [name{i} ' BF' num2str(j)]; - myname = replaceblanks(myname); - - names{end+1} = myname; + names{end + 1} = replaceblanks([condname{i} ' BF' num2str(j)]); %#ok end end -all_pmnames = {}; - for i = 1:nconds - is_pm = isfield(obj.Sess(s).U(i), 'P') && ~isempty(obj.Sess(s).U(i).P) && isfield(obj.Sess(s).U(i).P, 'P') && ~isempty(obj.Sess(s).U(i).P.P); - - if is_pm - % Get param mod names - % ------------------------------------------------ - pmnames = cell(1, nbf); - - for j = 1:nbf - myname = [obj.Sess(s).U(i).P.name ' BF' num2str(j)]; - myname = replaceblanks(myname); - - pmnames{1, j} = myname; - end - - all_pmnames = cat(2, all_pmnames, pmnames{:}); - + U = obj.Sess(s).U(i); + if ~local_has_pm(U), continue, end + nbf = size(local_bf_at_res(obj.xBF(min(i, numel(obj.xBF))), res), 2); + for j = 1:nbf + names{end + 1} = replaceblanks([U.P.name ' BF' num2str(j)]); %#ok end end names = names(:)'; -names = [names all_pmnames]; +end % get_session_X + + +% ========================================================================= +% Local helpers +% ========================================================================= +function dhr = local_hires_indicator(ons_sec, dur_sec, res, len_hires) +% Build a high-resolution onset/epoch indicator (len_hires x 1). Events with +% positive duration become epochs; zero-duration events are single impulses. +dhr = zeros(len_hires, 1); +for j = 1:numel(ons_sec) + k0 = round(ons_sec(j) * res) + 1; % first scan is time 0 -> element 1 + if k0 < 1 || k0 > len_hires, continue, end + if numel(dur_sec) >= j && dur_sec(j) > 0 + k1 = min(len_hires, round(k0 + dur_sec(j) * res)); + dhr(k0:k1) = 1; + else + dhr(k0) = 1; + end +end +end + +function bf = local_bf_at_res(xBF, res) +% Return the basis-set matrix resampled to res samples/second, using the +% stored sampling interval xBF.dt (seconds/sample). No-op when already at res. +bf = xBF.bf; +dt0 = []; +if isfield(xBF, 'dt') && ~isempty(xBF.dt), dt0 = xBF.dt; end +if isempty(dt0) || dt0 <= 0, dt0 = 1 / res; end +if abs(dt0 - 1 / res) < 1e-9, return, end + +n0 = size(bf, 1); +t0 = (0:n0 - 1)' * dt0; % times of stored samples (sec) +tnew = (0:1 / res:t0(end))'; % resample at res Hz +bf = interp1(t0, bf, tnew, 'linear', 'extrap'); end -% ---------------------------------------------- -% ---------------------------------------------- +function X = local_fit_rows(X, nscan) +% Trim or zero-pad X to exactly nscan rows. +if size(X, 1) >= nscan + X = X(1:nscan, :); +else + X(end + 1:nscan, :) = 0; +end +end -% ---------------------------------------------- -% ---------------------------------------------- +function tf = local_has_pm(U) +tf = isfield(U, 'P') && ~isempty(U.P) && isfield(U.P, 'P') && ~isempty(U.P.P); +end -function myname = replaceblanks(myname) -myname(myname == ' ') = '-'; -% this works in 2010b but not a... -%tmp = diff(num2str(myname)) == 0; +function v = local_getfield(s, f, default) +% Return s.(f) if the field exists and is non-empty, else default. +if isfield(s, f) && ~isempty(s.(f)) + v = s.(f); +else + v = default; +end +end + +function myname = replaceblanks(myname) +myname(myname == ' ') = '-'; tmp = diff(double(myname)) == 0; myname([false tmp] & myname == '-') = []; - end - diff --git a/CanlabCore/@fmri_glm_design_matrix/import_onsets.m b/CanlabCore/@fmri_glm_design_matrix/import_onsets.m index 17384787..595645ec 100644 --- a/CanlabCore/@fmri_glm_design_matrix/import_onsets.m +++ b/CanlabCore/@fmri_glm_design_matrix/import_onsets.m @@ -1,10 +1,224 @@ -function data_obj = import_onsets(data_obj, design_file) -% Import onsets from a csv or Excel file and attach them to the fmri_glm_design_matrix object - - % DesginmatrixTable INPUT - % TRY readtable() -- may be more flexible than xlsread. - xlsread(DesginmatrixTable) +function obj = import_onsets(obj, source, varargin) +% Import event onsets/durations/condition names (and parametric modulators) +% into an fmri_glm_design_matrix object. Two input styles are supported. +% +% (1) Tabular / FSL-style -- a file (.csv/.txt/.xlsx), or a MATLAB table, with +% one row per event and columns for onset, (optional) duration, and a +% condition label. The label column may be a string name OR an integer +% event-type code (condition and event type mean the same thing). Events +% are grouped by label into conditions (one Sess.U entry each). +% +% (2) SPM-style -- cell arrays with one cell per condition: a cell of onset +% vectors, and (optionally) matching cells of durations and parametric +% modulator values. +% +% :Usage: +% :: +% +% % Tabular / FSL-style +% obj = import_onsets(obj, 'events.csv') +% obj = import_onsets(obj, events_table, 'onset_col','Onset', 'name_col','EV') +% +% % SPM-style cell arrays (one cell per condition / event type) +% obj = import_onsets(obj, onsets_cell) +% obj = import_onsets(obj, onsets_cell, durations_cell) +% obj = import_onsets(obj, onsets_cell, durations_cell, pmods_cell, ... +% 'names', {'A','B'}, 'pm_names', {'rt',''}) +% +% :Inputs: +% +% **obj:** an fmri_glm_design_matrix object. +% +% **source:** +% A filename, a MATLAB table (tabular style), OR a cell array of onset +% vectors, one cell per condition (SPM style). +% +% :Optional Inputs (tabular style): +% +% **'onset_col' / 'dur_col' / 'name_col', columnname:** +% Override the column used for onsets / durations / condition labels. +% Defaults are matched case-insensitively against common names +% (onset/duration/name/condition/trial_type/event_type/value/code). +% +% :Optional Inputs (SPM style): +% +% **second/third positional cells:** +% durations (cell, one per condition) and parametric-modulator values +% (cell, one per condition); pass [] to skip either. +% +% **'names', {...}:** condition names, one per condition. +% **'pm_names', {...}:** parametric-modulator names, one per condition. +% +% :Outputs: +% +% **obj:** with obj.Sess(1).U populated (one entry per condition). +% +% :Examples: +% :: +% +% d = fmri_glm_design_matrix(2, 'nscan', 200, 'units', 'secs'); +% d = import_onsets(d, 'events.csv'); % FSL/tabular +% d = import_onsets(d, {[10 40]' [25 55]'}, {4 4}); % SPM-style +% d = build(d); +% +% :See also: fmri_glm_design_matrix.add, Add_Event_Info, readtable +% +% .. +% 2026 - Reimplemented (previous version was a non-functional stub) and +% extended to FSL/tabular files with event-type codes and SPM-style cells. +% .. +if iscell(source) + obj = local_import_spm_cells(obj, source, varargin{:}); +elseif istable(source) || ischar(source) || isstring(source) + obj = local_import_table(obj, source, varargin{:}); +else + error('fmri_glm_design_matrix:import_onsets', ... + 'source must be a filename, a table, or a cell array of onset vectors.'); +end + +end % import_onsets + + +% ========================================================================= +% SPM-style cell arrays +% ========================================================================= +function obj = local_import_spm_cells(obj, onsets, varargin) + +% Leading positional cell/[] args are durations then pmods +durations = {}; pmods = {}; +k = 0; +while ~isempty(varargin) && (iscell(varargin{1}) || isempty(varargin{1})) && ~ischar(varargin{1}) + k = k + 1; + if k == 1, durations = varargin{1}; elseif k == 2, pmods = varargin{1}; else, break, end + varargin(1) = []; +end + +p = inputParser; +p.addParameter('names', {}, @iscell); +p.addParameter('pm_names', {}, @iscell); +p.parse(varargin{:}); +names = p.Results.names; pm_names = p.Results.pm_names; + +nconds = numel(onsets); + +for i = 1:nconds + + obj.Sess(1).U(i).ons = onsets{i}(:); + + % name + if i <= numel(names) && ~isempty(names{i}), obj.Sess(1).U(i).name = names{i}; + else, obj.Sess(1).U(i).name = sprintf('Cond%d', i); + end + + % duration + if i <= numel(durations) && ~isempty(durations{i}) + obj.Sess(1).U(i).dur = durations{i}(:); + else + obj.Sess(1).U(i).dur = zeros(numel(onsets{i}), 1); + end + + % parametric modulator (single per condition) + if i <= numel(pmods) && ~isempty(pmods{i}) + P = struct('name', '', 'P', pmods{i}(:), 'h', 1, 'dur', [], 'i', []); + if i <= numel(pm_names) && ~isempty(pm_names{i}), P.name = pm_names{i}; + else, P.name = sprintf('pm%d', i); + end + obj.Sess(1).U(i).P = P; + end +end + +obj.history{end + 1} = sprintf('import_onsets: %d conditions from SPM-style cell arrays', nconds); + +end + + +% ========================================================================= +% Tabular / FSL-style +% ========================================================================= +function obj = local_import_table(obj, source, varargin) + +p = inputParser; +p.addParameter('onset_col', '', @(x) ischar(x) || isstring(x)); +p.addParameter('dur_col', '', @(x) ischar(x) || isstring(x)); +p.addParameter('name_col', '', @(x) ischar(x) || isstring(x)); +p.parse(varargin{:}); +opt = p.Results; + +if istable(source) + T = source; +else + if exist(char(source), 'file') ~= 2 + error('fmri_glm_design_matrix:import_onsets', 'File not found: %s', char(source)); + end + T = readtable(char(source)); +end + +vn = T.Properties.VariableNames; + +onset_col = local_pick_col(vn, opt.onset_col, {'onset', 'onsets', 'ons', 'time'}); +name_col = local_pick_col(vn, opt.name_col, {'name', 'condition', 'cond', 'trial_type', ... + 'event_type', 'eventtype', 'type', 'value', 'code', 'ev'}); +dur_col = local_pick_col(vn, opt.dur_col, {'duration', 'durations', 'dur'}); + +if isempty(onset_col) + error('fmri_glm_design_matrix:import_onsets', ... + 'Could not find an onset column. Columns: %s. Use the ''onset_col'' option.', strjoin(vn, ', ')); +end +if isempty(name_col) + error('fmri_glm_design_matrix:import_onsets', ... + 'Could not find a condition/event-type column. Columns: %s. Use the ''name_col'' option.', strjoin(vn, ', ')); +end + +onsets = T.(onset_col); +labels = T.(name_col); + +% Labels may be string names or integer event-type codes; normalize to a +% cellstr (codes become 'Cond'). +if isnumeric(labels) + labels = arrayfun(@(v) sprintf('Cond%g', v), labels, 'UniformOutput', false); +elseif iscategorical(labels) || isstring(labels) + labels = cellstr(labels); +elseif ~iscell(labels) + labels = cellstr(string(labels)); +end + +if ~isempty(dur_col), durs = T.(dur_col); else, durs = zeros(size(onsets)); end + +% Group by condition label (first-appearance order) +[uconds, ~, grp] = unique(labels, 'stable'); +for c = 1:numel(uconds) + wh = (grp == c); + obj.Sess(1).U(c).name = uconds{c}; + obj.Sess(1).U(c).ons = onsets(wh); + obj.Sess(1).U(c).dur = durs(wh); end +obj.history{end + 1} = sprintf('import_onsets: %d conditions, %d events from %s', ... + numel(uconds), numel(onsets), local_src_name(source)); + +end + + +% ========================================================================= +function col = local_pick_col(varnames, override, candidates) +col = ''; +if ~isempty(override) + wh = find(strcmpi(varnames, override), 1); + if isempty(wh) + error('fmri_glm_design_matrix:import_onsets', 'Column ''%s'' not found.', char(override)); + end + col = varnames{wh}; + return +end +for k = 1:numel(candidates) + wh = find(strcmpi(varnames, candidates{k}), 1); + if ~isempty(wh), col = varnames{wh}; return, end +end +end + + +function s = local_src_name(source) +if istable(source), s = '(table)'; else, s = char(source); end +end diff --git a/CanlabCore/@fmri_glm_design_matrix/plot.m b/CanlabCore/@fmri_glm_design_matrix/plot.m index 45dcf2ca..d1f6a852 100644 --- a/CanlabCore/@fmri_glm_design_matrix/plot.m +++ b/CanlabCore/@fmri_glm_design_matrix/plot.m @@ -12,8 +12,6 @@ function plot(obj) TR = obj.xY.RT; -bf = obj.xBF.bf; - % Define sessions and number of conditions nsess = length(obj.Sess); diff --git a/CanlabCore/@fmri_glm_design_matrix/replace_basis_set.m b/CanlabCore/@fmri_glm_design_matrix/replace_basis_set.m index 8891908e..2b81c9de 100644 --- a/CanlabCore/@fmri_glm_design_matrix/replace_basis_set.m +++ b/CanlabCore/@fmri_glm_design_matrix/replace_basis_set.m @@ -25,15 +25,52 @@ % %save this to get info that is not typically in basis set until after % %model is built. -oldBF = obj.xBF(1); +oldBF = obj.xBF(1); -% SPM adds these things here, so we will too, for consistency -xBF_hires.T = oldBF.T; -xBF_hires.T0 = oldBF.T0; -xBF_hires.UNITS = oldBF.UNITS; -xBF_hires.Volterra = oldBF.Volterra; +% Expand the basis-set array to one entry per condition (filled with the +% current default), so that replacing one condition's basis leaves the others +% unchanged regardless of whether the model has been built yet. (build's +% check_model only replicates a length-1 xBF, which would otherwise overwrite +% all conditions with the replacement.) +nconds = condition_num; +if ~isempty(obj.Sess) && isfield(obj.Sess, 'U') && ~isempty(obj.Sess(1).U) + nconds = max(condition_num, length(obj.Sess(1).U)); +end +for c = length(obj.xBF) + 1 : nconds + obj.xBF(c) = oldBF; +end + +% SPM adds these things here, so we will too, for consistency. Inherit them +% from the existing basis when the replacement does not provide them. +for f = {'T', 'T0', 'UNITS', 'Volterra'} + if ~isfield(xBF_hires, f{1}) && isfield(oldBF, f{1}) + xBF_hires.(f{1}) = oldBF.(f{1}); + end +end + +% Assigning into the struct array obj.xBF(condition_num) requires identical +% field sets and order. Different basis-set builders (spm_get_bf, +% fmri_spline_basis, ...) return different fields, so harmonize first: add any +% missing fields (filled with []) to both sides, then match the ordering. +allflds = union(fieldnames(obj.xBF), fieldnames(xBF_hires), 'stable'); +obj.xBF = local_ensure_fields(obj.xBF, allflds); +xBF_hires = local_ensure_fields(xBF_hires, allflds); +xBF_hires = orderfields(xBF_hires, obj.xBF(1)); obj.xBF(condition_num) = xBF_hires; end + +% ========================================================================= +function s = local_ensure_fields(s, flds) +% Add any missing fields (default []) to every element of struct array s, +% then reorder fields to match flds. +for k = 1:numel(flds) + if ~isfield(s, flds{k}) + [s.(flds{k})] = deal([]); + end +end +s = orderfields(s, flds); +end + 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..9520e7d0 --- /dev/null +++ b/CanlabCore/@glm_map/build_design.m @@ -0,0 +1,84 @@ +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'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: +% :: +% +% obj = build_design(obj, varargin) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a non-empty .design (fmri_glm_design_matrix) +% that has conditions/onsets assigned. +% +% :Optional Inputs: +% +% **'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: +% 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 + +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/check_properties.m b/CanlabCore/@glm_map/check_properties.m new file mode 100644 index 00000000..e8de1c77 --- /dev/null +++ b/CanlabCore/@glm_map/check_properties.m @@ -0,0 +1,78 @@ +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 + +% Nested struct fields (mirror fmri_data.regress out) +if ~isstruct(obj.input_parameters), obj.input_parameters = struct(); end +if ~isstruct(obj.input_image_metadata), obj.input_image_metadata = struct(); end +if ~isstruct(obj.diagnostics), obj.diagnostics = struct(); 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/create_orthogonal_contrast_set.m b/CanlabCore/@glm_map/create_orthogonal_contrast_set.m new file mode 100644 index 00000000..f9cdb927 --- /dev/null +++ b/CanlabCore/@glm_map/create_orthogonal_contrast_set.m @@ -0,0 +1,100 @@ +function obj = create_orthogonal_contrast_set(obj, varargin) +% Assign an orthogonal set of contrasts spanning the regressors of interest. +% +% Uses create_orthogonal_contrast_set.m to build a set of mutually orthogonal, +% zero-sum contrasts that span the space of differences among the regressors +% of interest (obj.wh_interest), and stores them in obj.contrasts with +% placeholder names. Nuisance covariates and the intercept get zero weight. +% Any previously assigned contrasts are replaced. +% +% :Usage: +% :: +% +% obj = create_orthogonal_contrast_set(obj) +% obj = create_orthogonal_contrast_set(obj, 'names', {'c1','c2',...}) +% +% :Inputs: +% +% **obj:** +% A glm_map object whose regressors of interest are defined +% (obj.wh_interest has at least 2 true entries). +% +% :Optional Inputs: +% +% **'names', {cellstr}:** +% Names for the contrasts (default 'OrthC1', 'OrthC2', ...). +% +% :Outputs: +% +% **obj:** +% glm_map with obj.contrasts ([regressors x contrasts]) and +% obj.contrast_names set to the orthogonal set. +% +% :Examples: +% :: +% +% d = fmri_glm_design_matrix(2,'nscan',120,'units','secs', ... +% 'onsets',{[10 40]' [25 55]' [12 42]'}, 'condition_names',{'A','B','C'}); +% g = glm_map(d); g = build_design(g); +% g = create_orthogonal_contrast_set(g); % 2 orthogonal contrasts over A,B,C +% g = run_diagnostics(g); +% +% :See also: +% - create_orthogonal_contrast_set (function), add_contrasts, diagnostics, calcEfficiency +% +% .. +% 2026 - Initial implementation. +% .. + +names = {}; +wh = find(strcmpi(varargin, 'names')); +if ~isempty(wh), names = varargin{wh(1) + 1}; end + +whI = find(obj.wh_interest); + +if isempty(whI) + error('glm_map:NoInterestRegressors', '%s', local_interest_guidance()); +end +if numel(whI) < 2 + error('glm_map:TooFewInterest', ... + ['Orthogonal contrasts require at least 2 regressors of interest; found %d. ' ... + 'Add more event types of interest, or check obj.wh_interest.'], numel(whI)); +end + +% Orthogonal contrast set over the of-interest regressors (rows = contrasts) +Cint = create_orthogonal_contrast_set(numel(whI)); % [(n-1) x n], dispatches to the function + +% Embed into the full design width (zeros on nuisance / intercept columns) +ncon = size(Cint, 1); +Cfull = zeros(ncon, obj.num_regressors); +Cfull(:, whI) = Cint; + +if isempty(names) + names = arrayfun(@(i) sprintf('OrthC%d', i), 1:ncon, 'UniformOutput', false); +elseif numel(names) ~= ncon + error('glm_map:ContrastNames', 'Number of names (%d) must match number of contrasts (%d).', numel(names), ncon); +end + +% Store (replace any existing contrasts). glm_map stores [regressors x contrasts]. +obj.contrasts = Cfull'; +obj.contrast_names = names(:)'; + +obj.history{end + 1} = sprintf(['create_orthogonal_contrast_set: %d orthogonal contrast(s) ' ... + 'spanning %d regressors of interest'], ncon, numel(whI)); + +end % create_orthogonal_contrast_set + + +% ========================================================================= +function s = local_interest_guidance() +s = sprintf([ ... + 'No regressors of interest are defined for this glm_map.\n' ... + 'Regressors of interest are the task-event regressors.\n' ... + ' - Event / 1st-level designs: add event onsets and build the design\n' ... + ' (g = build_design(g), or import_onsets / import_SPM); event\n' ... + ' regressors are then flagged as of interest automatically.\n' ... + ' - Direct / group designs: by default all non-intercept columns are of\n' ... + ' interest. Mark covariates of no interest with\n' ... + ' g.nuisance_columns = [column indices]; at least one column must\n' ... + ' remain of interest.']); +end diff --git a/CanlabCore/@glm_map/fit.m b/CanlabCore/@glm_map/fit.m new file mode 100644 index 00000000..eab05b44 --- /dev/null +++ b/CanlabCore/@glm_map/fit.m @@ -0,0 +1,251 @@ +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 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: +% :: +% +% 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), and optionally contrasts. +% +% **data:** +% 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: +% +% **'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. +% +% **'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). +% +% :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, ... +% '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, threshold +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. Thin orchestration layer over +% fmri_data.regress: glm_map owns the design/diagnostics/result container, +% regress remains the compute engine. +% .. + +% ------------------------------------------------------------------------- +% 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 +% Always request residuals so diagnostics can compute Cook's distance; they +% are dropped again below unless the caller asked to keep them (do_resid). +regress_args{end + 1} = 'residual'; + +% 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. fmri_data.regress now returns a glm_map object whose +% properties mirror the historical out-struct fields. +% ------------------------------------------------------------------------- +out = regress(data, regress_args{:}); + +% ------------------------------------------------------------------------- +% Unpack result maps and the nested input structs +% ------------------------------------------------------------------------- +obj.betas = out.betas; +obj.t = out.t; +obj.sigma = out.sigma; +obj.df = out.df; +obj.input_parameters = out.input_parameters; +obj.input_image_metadata = out.input_image_metadata; +obj.contrast_summary_table = out.contrast_summary_table; + +% Error degrees of freedom (scalar summary; per-voxel df lives in obj.df) +if ~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 ~isempty(out.contrast_estimates) + obj.contrast_estimates = out.contrast_estimates; + obj.contrast_t = out.contrast_t; + obj.contrast_t.dfe = obj.dfe; +end + +% Stash residuals so run_diagnostics() can compute Cook's distance. If the caller +% did not ask to keep them, they are cleared after diagnostics (below). +if ~isempty(out.residuals) + obj.residuals = out.residuals; +end + +% Sync names back from regress (it may have generated/added them) +if isempty(obj.design) + obj.regressor_names_direct = out.regressor_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, Cook's distance, condition +% number, collinearity report; uses canonical VIF/cVIF rather than getvif) +obj = run_diagnostics(obj, 'noverbose'); + +% Drop the (potentially large) residual maps unless the caller kept them +if ~do_resid + obj.residuals = []; +end + +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/glm_map.m b/CanlabCore/@glm_map/glm_map.m new file mode 100644 index 00000000..4fa3d9fd --- /dev/null +++ b/CanlabCore/@glm_map/glm_map.m @@ -0,0 +1,740 @@ +% 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 +% g = glm_map(regress_out_struct) % re-cast a regress() output struct +% +% The object property names mirror the fields of the results structure +% returned by fmri_data.regress (the variable "out" in that method). Related +% options are grouped into nested structs (.input_parameters, +% .input_image_metadata, .diagnostics), and per-voxel df / sigma are kept as +% fmri_data objects in .df and .sigma. For backward compatibility the +% historical out-struct field names are available as aliases that read/write +% the canonical properties: +% +% .b -> .betas .contrast_images -> .contrast_estimates .con_t -> .contrast_t +% .resid -> .residuals .variable_names -> .regressor_names .C -> .contrasts +% +% :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 +% run_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 + + % --- Provenance / metadata -------------------------------------- + analysis_name = ''; % Short descriptive name for this analysis (mirrors regress out.analysis_name) + + % --- 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 (read/written as out.C via the .C alias) + contrast_names = {}; % Cell array of contrast names, one per column of C (mirrors out.contrast_names) + contrast_summary_table = table(); % Table summarizing contrasts over conditions (mirrors out.contrast_summary_table) + + % --- Nested input/option structs (mirror fmri_data.regress out) -- + input_parameters = struct(); % Struct: options used by the fit (brain_is_predictor, do_robust, grandmeanscale, do_intercept, do_resid, initial_statistical_threshold, ...). Mirrors out.input_parameters. + input_image_metadata = struct(); % Struct: provenance of the fitted images (source_notes, history, image_names, fullpath). Mirrors out.input_image_metadata. + + % --- Fitted result maps (populated by fit / ingested from regress) - + betas % statistic_image (type 'Beta'), [voxels x regressors]. Alias: .b + t % statistic_image (type 'T') for regressors, [voxels x regressors] + contrast_estimates % statistic_image (type 'Contrast'), [voxels x contrasts]. Alias: .contrast_images + contrast_t % statistic_image (type 'T') for contrasts, [voxels x contrasts]. Alias: .con_t + df % fmri_data object: per-voxel error degrees of freedom (mirrors out.df) + sigma % fmri_data object: residual standard deviation per voxel (mirrors out.sigma) + residuals % fmri_data object: residuals (optional; only if requested). Alias: .resid + dfe % Scalar error degrees of freedom for the fit (median of df, convenience summary) + + % --- Design diagnostics (nested struct; populated by fit/run_diagnostics) - + % Fields: Variance_inflation_factors, Leverages (same names as + % fmri_data.regress out.diagnostics), plus + % Contrast_variance_inflation_factors, Cooks_distance (per-observation + % influence; populated when residuals are available), condition_number, + % rank_deficient, collinearity_report, vif_threshold. + diagnostics = struct(); + + warnings = {}; % Cell array of warning messages accumulated during build/fit (mirrors out.warnings) + + % --- Provenance / metadata -------------------------------------- + fit_parameters = struct(); % Struct recording options used in fit (robust, AR order, threshold, ...) + 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 + + % Column indices of X that are nuisance covariates (covariates of no + % interest). Used mainly in direct/group mode, where there is no + % event/covariate partition to read; in event mode the partition is + % taken from design.xX (iH = of interest, iC/iG = nuisance, iB = + % intercept). See the Dependent .wh_interest / .wh_nuisance accessors. + nuisance_columns = []; + + property_descriptions = { ... + 'analysis_name: short descriptive name (out.analysis_name)' ... + '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 (.C): [regressors x contrasts] contrast matrix' ... + 'contrast_names: cell array of contrast names' ... + 'contrast_summary_table: table of contrast weights over conditions (out.contrast_summary_table)' ... + 'input_parameters: struct of fit options (out.input_parameters)' ... + 'input_image_metadata: struct of input-image provenance (out.input_image_metadata)' ... + 'betas (.b) / t: statistic_image maps of regression coefficients and their t-statistics, [voxels x regressors]' ... + 'contrast_estimates (.contrast_images) / contrast_t (.con_t): statistic_image maps for linear contrasts, [voxels x contrasts]' ... + 'df: fmri_data object with per-voxel error degrees of freedom (out.df)' ... + 'sigma: fmri_data object with residual standard deviation per voxel (out.sigma)' ... + 'residuals (.resid): fmri_data object with residuals (optional)' ... + 'dfe: scalar error degrees of freedom (median of df)' ... + 'diagnostics: struct with VIF/cVIF (full + interest-only), leverages, Cooks_distance, condition_number, rank_deficient, collinearity_report' ... + 'nuisance_columns: column indices of X that are nuisance covariates (direct mode). See wh_interest/wh_nuisance/wh_intercept.' ... + '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 + + % --- Regressor-role indicators (logical, [1 x num_regressors]) ---- + wh_interest % true for regressors of interest (task events). Event mode: design.xX.iH. Direct mode: non-intercept, non-nuisance columns. + wh_nuisance % true for nuisance covariates (of no interest). Event mode: design.xX.iC/iG. Direct mode: obj.nuisance_columns. + wh_intercept % true for intercept/baseline columns. Event mode: design.xX.iB. Direct mode: constant columns. + + % --- Aliases to the fmri_data.regress out-struct field names ----- + % These read/write the canonical properties above, so that an object + % returned by fmri_data.regress supports the historical struct-style + % field access (out.b, out.con_t, ...) unchanged. + b % Alias for betas (out.b) + contrast_images % Alias for contrast_estimates (out.contrast_images) + con_t % Alias for contrast_t (out.con_t) + resid % Alias for residuals (out.resid) + variable_names % Alias for regressor_names (out.variable_names) + C % Alias for contrasts (out.C) + + end % dependent properties + + + methods + + % ================================================================= + % Constructor + % ================================================================= + function obj = glm_map(varargin) + + if nargin > 0 + + % If the first argument is a struct, treat it as a regression + % results structure (the historical fmri_data.regress output) + % and re-cast it as a glm_map object. Any remaining 'field', + % value pairs are then applied as overrides. + if ~isempty(varargin) && isstruct(varargin{1}) + obj = local_from_regress_struct(obj, varargin{1}); + varargin(1) = []; + + % If first argument is an fmri_glm_design_matrix, wrap it as a + % 1st-level (event) design and consume that argument. + elseif ~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'); % includes public Dependent props + + % Names of Dependent properties that have setters + settable_dependent = {'TR', 'X', 'regressor_names', ... + 'b', 'contrast_images', 'con_t', 'resid', 'variable_names', 'C'}; + + 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 % nargin > 0 + + % Populate every nested option/diagnostic struct with its full set + % of valid fields (missing fields are filled with []), so a freshly + % constructed object always exposes the complete schema. + obj = validate_object(obj); + + 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) && isscalar(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) && 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) + 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 + + function val = get.wh_interest(obj) + [val, ~, ~] = local_partition(obj); + end + + function val = get.wh_nuisance(obj) + [~, val, ~] = local_partition(obj); + end + + function val = get.wh_intercept(obj) + [~, ~, val] = local_partition(obj); + 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) && 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 ' ... + '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 + + + % ================================================================= + % Aliases to the fmri_data.regress out-struct field names + % (read/write the canonical properties) + % ================================================================= + function val = get.b(obj), val = obj.betas; end + function obj = set.b(obj, val), obj.betas = val; end + + function val = get.contrast_images(obj), val = obj.contrast_estimates; end + function obj = set.contrast_images(obj, val), obj.contrast_estimates = val; end + + function val = get.con_t(obj), val = obj.contrast_t; end + function obj = set.con_t(obj, val), obj.contrast_t = val; end + + function val = get.resid(obj), val = obj.residuals; end + function obj = set.resid(obj, val), obj.residuals = val; end + + function val = get.variable_names(obj), val = obj.regressor_names; end + function obj = set.variable_names(obj, val), obj.regressor_names = val; end + + function val = get.C(obj), val = obj.contrasts; end + function obj = set.C(obj, val), obj.contrasts = val; end + + + % ================================================================= + % disp: object summary with full property listing + % ================================================================= + function disp(obj) + + line = repmat('-', 1, 64); + fprintf(' glm_map object\n %s\n', line); + + % -------- One-line status header -------- + 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 | fitted: %s\n', levelstr, local_tf(obj.is_fitted, 'YES', 'no')); + + % Design matrix (obj.X reads through to design.xX.X in event mode) + nI = sum(obj.wh_interest); nN = sum(obj.wh_nuisance); nB = sum(obj.wh_intercept); + fprintf(' design X: %d observations x %d regressors (%d of interest, %d nuisance, %d intercept)\n', ... + obj.num_images, obj.num_regressors, nI, nN, nB); + fprintf(' contrasts: %d\n', obj.num_contrasts); + fprintf(' %s\n', line); + + % -------- Full property listing -------- + % Curated order; nested structs are expanded one level so their + % fields are visible at a glance. (See properties(obj) for the + % full set, including the out-struct aliases b/con_t/.../C.) + props = { ... + 'analysis_name', 'design', 'level', 'is_timeseries', ... + 'contrasts', 'contrast_names', 'contrast_summary_table', ... + 'input_parameters', 'input_image_metadata', ... + 'betas', 't', 'contrast_estimates', 'contrast_t', ... + 'df', 'sigma', 'residuals', 'dfe', ... + 'diagnostics', 'warnings', ... + 'fit_parameters', 'notes', 'history'}; + + for i = 1:numel(props) + p = props{i}; + fprintf(' %-22s : %s\n', p, local_summarize(obj.(p))); + + % Expand the three nested option/diagnostic structs one level + if any(strcmp(p, {'input_parameters', 'input_image_metadata', 'diagnostics'})) ... + && isstruct(obj.(p)) && ~isempty(fieldnames(obj.(p))) + fn = fieldnames(obj.(p)); + for j = 1:numel(fn) + fprintf(' %-22s .%-18s %s\n', '', fn{j}, ... + local_summarize(obj.(p).(fn{j}))); + end + end + end + + fprintf(' %s\n', line); + fprintf(' methods(glm_map) for operations; properties(glm_map) for all fields.\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 [interest, nuisance, intercept] = local_partition(obj) +% Partition the columns of obj.X into regressors of interest (task events), +% nuisance covariates, and intercept/baseline columns. Returns three logical +% row vectors of length num_regressors. + +k = obj.num_regressors; +[interest, nuisance, intercept] = deal(false(1, k)); +if k == 0, return, end + +% --- Event mode: read the partition from the built design (xX.iH/iC/iG/iB) - +if ~isempty(obj.design) && isstruct(obj.design.xX) && isscalar(obj.design.xX) ... + && isfield(obj.design.xX, 'X') && ~isempty(obj.design.xX.X) ... + && size(obj.design.xX.X, 2) == k && isfield(obj.design.xX, 'iH') + + xX = obj.design.xX; + interest(local_idx(xX, 'iH')) = true; + nuisance(local_idx(xX, 'iC')) = true; + nuisance(local_idx(xX, 'iG')) = true; + intercept(local_idx(xX, 'iB')) = true; + + % Any column not assigned by the partition defaults to "of interest" + unassigned = ~(interest | nuisance | intercept); + interest(unassigned) = true; + return +end + +% --- Direct mode: detect intercept (constant columns); use stored nuisance - +X = obj.X; +isconst = all(X == X(1, :), 1) & any(X ~= 0, 1); % constant, nonzero columns +intercept(isconst) = true; + +nz = obj.nuisance_columns; +if ~isempty(nz) + if islogical(nz) + nuisance(1:min(k, numel(nz))) = nz(1:min(k, numel(nz))); + else + nz = nz(nz >= 1 & nz <= k); + nuisance(nz) = true; + end +end +nuisance = nuisance & ~intercept; +interest = ~(intercept | nuisance); + +end % local_partition + + +function idx = local_idx(xX, f) +% Safely return field xX.(f) as a row of indices (empty if absent). +idx = []; +if isfield(xX, f) && ~isempty(xX.(f)), idx = xX.(f)(:)'; end +end + + +function s = local_tf(tf, a, b) +% Inline if: return a if tf is true, else b ('' if b omitted). +if nargin < 3, b = ''; end +if tf, s = a; else, s = b; end +end % local_tf + + +function s = local_summarize(v) +% Compact one-line description of any property value, for disp(). + +if isempty(v) + s = '[]'; + return +end + +if ischar(v) + if size(v, 1) > 1 + s = sprintf('[%dx%d char]', size(v, 1), size(v, 2)); % char matrix + else + s = v; + if numel(s) > 52, s = [s(1:49) '...']; end + end + return +end + +if islogical(v) && isscalar(v) + s = local_tf(v, 'true', 'false'); + return +end + +if isnumeric(v) && isscalar(v) + s = num2str(v); + return +end + +if isa(v, 'statistic_image') || isa(v, 'fmri_data') + try + s = sprintf('%s [%d voxels x %d images]', class(v), size(v.dat, 1), size(v.dat, 2)); + catch + s = class(v); + end + return +end + +if isa(v, 'fmri_glm_design_matrix') + try + s = sprintf('fmri_glm_design_matrix (TR = %g)', v.TR); + catch + s = 'fmri_glm_design_matrix'; + end + return +end + +if istable(v) + s = sprintf('table [%d x %d]', size(v, 1), size(v, 2)); + return +end + +if isstruct(v) + fn = fieldnames(v); + if isempty(fn) + s = 'struct (empty)'; + else + s = sprintf('struct with %d field(s)', numel(fn)); + end + return +end + +if iscell(v) + % Show contents if it is a short cellstr + if numel(v) <= 6 && all(cellfun(@(x) ischar(x) || isstring(x), v)) + s = ['{' strjoin(cellfun(@char, v(:)', 'UniformOutput', false), ', ') '}']; + if numel(s) > 52, s = sprintf('{%dx%d cell}', size(v, 1), size(v, 2)); end + else + s = sprintf('{%dx%d cell}', size(v, 1), size(v, 2)); + end + return +end + +if isnumeric(v) + s = sprintf('[%dx%d %s]', size(v, 1), size(v, 2), class(v)); + return +end + +s = sprintf('[%s]', class(v)); + +end % local_summarize + + +function obj = local_from_regress_struct(obj, S) +% Re-cast a regression results structure (the historical fmri_data.regress +% output) into a glm_map object. Maps the out-struct field names onto the +% canonical glm_map properties; tolerates missing fields. + +obj.level = 2; % regress output carries no level; default to group + +% src field in struct -> dst glm_map property (or settable alias) +map = { ... + 'analysis_name', 'analysis_name'; ... + 'input_parameters', 'input_parameters'; ... + 'input_image_metadata', 'input_image_metadata'; ... + 'X', 'X'; ... % set.X -> Xdirect + 'variable_names', 'regressor_names'; ... % set.regressor_names -> regressor_names_direct + 'C', 'contrasts'; ... + 'contrast_names', 'contrast_names'; ... + 'contrast_summary_table', 'contrast_summary_table'; ... + 'warnings', 'warnings'; ... + 'b', 'betas'; ... + 't', 't'; ... + 'df', 'df'; ... + 'sigma', 'sigma'; ... + 'resid', 'residuals'; ... + 'contrast_images', 'contrast_estimates'; ... + 'con_t', 'contrast_t'}; + +for i = 1:size(map, 1) + src = map{i, 1}; + if ~isfield(S, src), continue, end + v = S.(src); + % Skip genuinely-empty primitives, but always ingest result-map objects + % (statistic_image / fmri_data may overload isempty based on volInfo). + if (isnumeric(v) || ischar(v) || iscell(v) || isstruct(v)) && isempty(v) + continue + end + obj.(map{i, 2}) = v; +end + +% Diagnostics: ingest as-is (regress already names the fields +% Variance_inflation_factors / Leverages, which glm_map reuses) +if isfield(S, 'diagnostics') && ~isempty(S.diagnostics) + obj.diagnostics = S.diagnostics; +end + +% Scalar error-df summary from the per-voxel df image +if isfield(S, 'df') && ~isempty(S.df) && isprop(S.df, 'dat') && ~isempty(S.df.dat) + obj.dfe = double(median(S.df.dat(:))); +end + +obj.history{end + 1} = 'constructed from fmri_data.regress() output structure'; + +end % local_from_regress_struct + + diff --git a/CanlabCore/@glm_map/import_SPM.m b/CanlabCore/@glm_map/import_SPM.m new file mode 100644 index 00000000..7d17753d --- /dev/null +++ b/CanlabCore/@glm_map/import_SPM.m @@ -0,0 +1,248 @@ +function obj = import_SPM(obj, SPM, varargin) +% 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. +% 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') % path to SPM.mat +% obj = import_SPM(glm_map, '/path/to/dir') % directory containing SPM.mat +% +% :Inputs: +% +% **obj:** +% A glm_map object (typically empty: glm_map). +% +% **SPM:** +% 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/.img images into obj.betas +% (looked up in SPM.swd), labeled by SPM.xX.name. +% +% **'noverbose':** +% Suppress console output. +% +% :Outputs: +% +% **obj:** +% 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, fit +% +% .. +% Programmers' notes: +% 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, ...). +% .. + +% ------------------------------------------------------------------------- +% 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) + +% Derive the of-interest / nuisance / intercept partition from SPM's event +% semantics. In SPM first-level models the task-event regressors are named +% '... *bf()' (and parametric mods '... ^

*bf()'), session +% constants are named '... constant', and movement / multiple_regressors +% columns are everything else. SPM's own xX.iH is typically empty for +% event-related first-level designs, so we set the partition here (preserving +% SPM's original indices under spm_*). glm_map's wh_interest reads xX.iH. +design.xX = local_spm_partition(design.xX); + +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 xX = local_spm_partition(xX) +% Set xX.iH (of interest), xX.iC (nuisance), xX.iB (intercept/constant), +% xX.iG (empty) from SPM regressor-name semantics. Preserve SPM's original +% partition indices under spm_iH / spm_iC / spm_iB / spm_iG. + +xX.spm_iH = local_getfield(xX, 'iH'); +xX.spm_iC = local_getfield(xX, 'iC'); +xX.spm_iB = local_getfield(xX, 'iB'); +xX.spm_iG = local_getfield(xX, 'iG'); + +X = xX.X; +ncol = size(X, 2); +nm = {}; +if isfield(xX, 'name'), nm = xX.name; end + +is_event = false(1, ncol); +is_const = false(1, ncol); +for j = 1:ncol + if j <= numel(nm) && (ischar(nm{j}) || isstring(nm{j})) + s = char(nm{j}); + if contains(s, '*bf('), is_event(j) = true; end % convolved task event + if contains(lower(s), 'constant'), is_const(j) = true; end % session constant + end +end + +% Numeric fallback for constants (all-equal nonzero columns) +const_cols = all(X == X(1, :), 1) & any(X ~= 0, 1); +is_const = is_const | const_cols; +is_event = is_event & ~is_const; + +xX.iH = find(is_event); % of interest +xX.iB = find(is_const); % intercept / session constants +xX.iG = []; +xX.iC = find(~is_event & ~is_const); % nuisance (movement, covariates, ...) + +end % local_spm_partition + + +function v = local_getfield(s, f) +if isfield(s, f), v = s.(f); else, v = []; end +end + + +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 diff --git a/CanlabCore/@glm_map/import_onsets.m b/CanlabCore/@glm_map/import_onsets.m new file mode 100644 index 00000000..26eb01ef --- /dev/null +++ b/CanlabCore/@glm_map/import_onsets.m @@ -0,0 +1,150 @@ +function obj = import_onsets(obj, source, varargin) +% Import event onsets into a glm_map's wrapped fmri_glm_design_matrix design. +% +% Passes its inputs through to fmri_glm_design_matrix.import_onsets (which +% accepts FSL/tabular files or tables, and SPM-style cell arrays), then builds +% the design matrix if enough information is available. Event regressors are +% flagged as of interest automatically (they form the design's H partition); +% any covariates of no interest become nuisance regressors. If the glm_map has +% no design yet, you can bootstrap one by passing 'TR' and 'nscan'; otherwise +% the method explains what is needed and how to attach it. +% +% :Usage: +% :: +% +% g = import_onsets(g, source, ...) % design already attached +% g = import_onsets(g, source, 'TR', 2, 'nscan', 200) % bootstrap the design +% +% :Inputs: +% +% **obj:** a glm_map object. +% **source:** a filename / table (FSL/tabular) or a cell of onset vectors +% (SPM-style); see fmri_glm_design_matrix.import_onsets. +% +% :Optional Inputs: +% +% **'TR', value / 'nscan', value / 'units', 'secs'|'scans':** +% Used to create (or update) the wrapped design when bootstrapping. +% Remaining options are passed through to the design's import_onsets. +% +% :Outputs: +% +% **obj:** glm_map with the design imported and (if possible) built; level +% set to 1. obj.wh_interest then flags the event regressors. +% +% :Examples: +% :: +% +% g = import_onsets(glm_map, 'events.csv', 'TR', 2, 'nscan', 200); +% g = run_diagnostics(g); +% +% :See also: +% - fmri_glm_design_matrix.import_onsets, build_design, import_SPM +% +% .. +% 2026 - Initial implementation. +% .. + +% ------------------------------------------------------------------------- +% Pull out design-bootstrap options (TR / nscan / units); pass the rest on +% ------------------------------------------------------------------------- +[TR, nscan, units, rest] = local_extract_design_opts(varargin); + +% ------------------------------------------------------------------------- +% Ensure a design object exists +% ------------------------------------------------------------------------- +if isempty(obj.design) || ~isa(obj.design, 'fmri_glm_design_matrix') + if ~isempty(TR) && ~isempty(nscan) + d = fmri_glm_design_matrix(TR, 'nscan', nscan); + if ~isempty(units), d = add(d, 'units', units); end + obj.design = d; + else + error('glm_map:NoDesignForImport', '%s', local_design_guidance()); + end +else + if ~isempty(TR), obj.design.TR = TR; obj.design.xY.RT = TR; end + if ~isempty(nscan), obj.design.nscan = nscan; end + if ~isempty(units), obj.design = add(obj.design, 'units', units); end +end + +% ------------------------------------------------------------------------- +% Delegate the import to the wrapped design object +% ------------------------------------------------------------------------- +obj.design = import_onsets(obj.design, source, rest{:}); +obj.level = 1; + +% ------------------------------------------------------------------------- +% Build if we have enough information; otherwise direct the user +% ------------------------------------------------------------------------- +if local_can_build(obj.design) + obj.design = build(obj.design); + obj = validate_object(obj); + obj.history{end + 1} = sprintf('import_onsets: imported and built design [%d x %d]', ... + size(obj.X, 1), size(obj.X, 2)); +else + obj.history{end + 1} = 'import_onsets: imported onsets; design not built (missing info)'; + fprintf('%s\n', local_missing_info(obj.design)); +end + +end % import_onsets + + +% ========================================================================= +% Local helpers +% ========================================================================= +function [TR, nscan, units, rest] = local_extract_design_opts(args) +[TR, nscan, units] = deal([]); +rest = args; +for key = {'TR', 'nscan', 'units'} + wh = find(strcmpi(rest, key{1}), 1); + if ~isempty(wh) && wh < numel(rest) + switch lower(key{1}) + case 'tr', TR = rest{wh + 1}; + case 'nscan', nscan = rest{wh + 1}; + case 'units', units = rest{wh + 1}; + end + rest(wh:wh + 1) = []; + end +end +end + + +function tf = local_can_build(d) +% Buildable if TR is set, nscan is set, and at least one condition has onsets. +tf = false; +if isempty(d) || ~isa(d, 'fmri_glm_design_matrix'), return, end +if isempty(d.nscan) || any(isnan(d.TR)), return, end +if isempty(d.Sess) || ~isfield(d.Sess(1), 'U') || isempty(d.Sess(1).U), return, end +if ~isfield(d.Sess(1).U(1), 'ons') || isempty(d.Sess(1).U(1).ons), return, end +tf = true; +end + + +function s = local_missing_info(d) +miss = {}; +if isempty(d.nscan) || any(isnan(d.TR)) + if any(isnan(d.TR)), miss{end + 1} = 'TR (repetition time)'; end + if isempty(d.nscan), miss{end + 1} = 'nscan (scans per session)'; end +end +if isempty(d.Sess) || ~isfield(d.Sess(1), 'U') || isempty(d.Sess(1).U) ... + || ~isfield(d.Sess(1).U(1), 'ons') || isempty(d.Sess(1).U(1).ons) + miss{end + 1} = 'event onsets'; +end +s = sprintf([ ... + ' import_onsets: design imported but not built; still missing: %s.\n' ... + ' Attach the missing info and build, e.g.:\n' ... + ' g.design.nscan = nscan; %% number of scans per session\n' ... + ' g.design.TR = TR; %% repetition time (s)\n' ... + ' g = build_design(g);'], strjoin(miss, ', ')); +end + + +function s = local_design_guidance() +s = sprintf([ ... + 'glm_map has no design to import onsets into.\n' ... + 'Either bootstrap one by passing TR and nscan:\n' ... + ' g = import_onsets(g, source, ''TR'', 2, ''nscan'', 200, ''units'', ''secs'');\n' ... + 'or attach an fmri_glm_design_matrix first:\n' ... + ' g.design = fmri_glm_design_matrix(2, ''nscan'', 200, ''units'', ''secs'');\n' ... + ' g = import_onsets(g, source);']); +end diff --git a/CanlabCore/@glm_map/montage.m b/CanlabCore/@glm_map/montage.m new file mode 100644 index 00000000..e7b1e593 --- /dev/null +++ b/CanlabCore/@glm_map/montage.m @@ -0,0 +1,62 @@ +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, 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 + 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..95a3c6d6 --- /dev/null +++ b/CanlabCore/@glm_map/plot_design.m @@ -0,0 +1,306 @@ +function plot_design(obj, varargin) +% Plot the design of a glm_map object. +% +% If the model is an event/1st-level design with a small number of event +% types of interest (<= 12), each event type is drawn in its own panel showing +% the object's *actual* basis-convolved regressor(s) -- one line per basis +% function, so the real basis set is respected -- with boxes marking each event +% and its duration. A heat map of the full design matrix (including any +% nuisance covariates and the intercept) is shown alongside. Otherwise +% (direct/group designs, or many event types) only the full-design heat map is +% shown. Heat maps use 'YDir','reverse' (observation 0 at the top) and axis +% tight. +% +% :Usage: +% :: +% +% plot_design(obj) +% plot_design(obj, 'max_event_types', 12) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design available (obj.X non-empty). +% +% :Optional Inputs: +% +% **'max_event_types', [k]:** +% Use the per-event-type line-plot format only when the number of event +% types of interest is <= k (default 12). +% +% :See also: +% - diagnostics, fmri_glm_design_matrix.plot, drawbox +% +% .. +% 2026 - Renders the object's own basis-convolved regressors per event type +% (one line per basis function) instead of re-convolving with a canonical +% HRF, so the actual basis set is respected. +% .. + +X = obj.X; +if isempty(X) + error('glm_map:NoDesign', 'No design matrix available (obj.X is empty).'); +end + +max_event_types = 12; +wh = find(strcmpi(varargin, 'max_event_types')); +if ~isempty(wh), max_event_types = varargin{wh(1) + 1}; end + +% Per-event-type info for session 1 (onsets/durations/names + the X columns +% that hold each condition's basis-function regressors). +ev = local_event_info(obj); +do_panels = ~isempty(ev) && numel(ev) <= max_event_types; + +if do_panels + nconds = numel(ev); + colors = local_colors(nconds); + + create_figure('glm_map design'); + + % Left column: one panel per event type. Right column: full-design heat map. + for i = 1:nconds + subplot(nconds, 2, 2 * i - 1); + local_event_panel(ev(i), colors{i}, i == nconds); + end + + subplot(nconds, 2, 2:2:2 * nconds); + local_heatmap(obj, X); +else + create_figure('glm_map design'); + local_heatmap(obj, X); +end + +% If diagnostics have been computed, add a VIF / cVIF figure. +local_vif_figure(obj); + +end % plot_design + + +% ========================================================================= +% Local helpers +% ========================================================================= +function ev = local_event_info(obj) +% Build a struct array (one per event type of interest, session 1) with: +% .name, .time (sec), .Xcols ([nscan1 x nbf] regressor block), .ons, .dur +% Returns [] if this is not a built event design. +ev = []; + +if isempty(obj.design) || ~isa(obj.design, 'fmri_glm_design_matrix'), return, end +if isempty(obj.design.Sess) || ~isfield(obj.design.Sess(1), 'U') || isempty(obj.design.Sess(1).U), return, end + +X = obj.X; +if isempty(X), return, end + +U = obj.design.Sess(1).U; +TR = obj.TR; + +ns1 = obj.design.nscan(1); +if isempty(ns1) || ns1 < 1 || ns1 > size(X, 1), ns1 = size(X, 1); end +time = (0:ns1 - 1)' * TR; + +% onset/duration units -> seconds +to_sec = 1; +if ~isempty(obj.design.xBF) && isfield(obj.design.xBF(1), 'UNITS') ... + && any(strcmpi(obj.design.xBF(1).UNITS, {'scans', 'tr', 'trs'})) + to_sec = TR; +end + +% Interest basis-function columns are ordered condition-by-condition at the +% front of X (session 1 block), nbf columns per condition. +col = 0; +nconds = numel(U); +ev = repmat(struct('name', '', 'time', time, 'Xcols', [], 'ons', [], 'dur', []), 1, nconds); + +for i = 1:nconds + nbf = size(obj.design.xBF(min(i, numel(obj.design.xBF))).bf, 2); + cols = col + (1:nbf); + col = col + nbf; + if any(cols > size(X, 2)), cols = cols(cols <= size(X, 2)); end + + nm = U(i).name; if iscell(nm), if isempty(nm), nm = sprintf('Cond%d', i); else, nm = nm{1}; end, end + + o = U(i).ons(:) * to_sec; + dd = []; if isfield(U(i), 'dur'), dd = U(i).dur; end + if isempty(dd), dd = 0; end + dd = dd(:); if isscalar(dd), dd = repmat(dd, numel(o), 1); end + dd = dd * to_sec; + + ev(i).name = nm; + ev(i).Xcols = X(1:ns1, cols); + ev(i).ons = o; + ev(i).dur = dd; +end +end + + +function local_event_panel(e, color, is_last) +% One panel: the event type's actual regressor line(s) (one per basis +% function) plus boxes marking each event and its duration. + +Xi = e.Xcols; +nbf = size(Xi, 2); + +% Regressor lines: solid color, slightly varying line style per basis function +hold on; +styles = {'-', '--', ':', '-.'}; +for j = 1:nbf + plot(e.time, Xi(:, j), 'Color', color, 'LineStyle', styles{mod(j - 1, numel(styles)) + 1}, 'LineWidth', 1.5); +end + +% Event boxes along the bottom of the panel +yl = [min([Xi(:); 0]) max([Xi(:); eps])]; +if diff(yl) == 0, yl = yl + [-1 1]; end +boxh = 0.12 * diff(yl); +ybot = yl(1) - 1.6 * boxh; +minw = max(0.4 * (e.time(2) - e.time(1)), 0.001); % give impulses a visible width +for k = 1:numel(e.ons) + w = max(e.dur(k), minw); + local_drawbox(e.ons(k), w, ybot, boxh, color); +end + +axis tight; +ylim([ybot - 0.4 * boxh, yl(2) + 0.05 * diff(yl)]); +ylabel(e.name, 'Interpreter', 'none', 'Rotation', 0, 'HorizontalAlignment', 'right'); +set(gca, 'YTick', []); +% Time axis is always in seconds (onsets/durations converted from the design's +% units, and the regressor time base is (0:nscan-1)*TR). +if is_last, xlabel('Time (seconds)'); else, set(gca, 'XTickLabel', []); end +if nbf > 1 + legend(arrayfun(@(j) sprintf('BF%d', j), 1:nbf, 'UniformOutput', false), ... + 'Location', 'northeast', 'Box', 'off'); +end +hold off; +end + + +function local_drawbox(t, dur, ystart, yheight, color) +% Filled rectangle for one event (borrowed from drawbox.m). +x = [0 1 1 0] * dur + t; +y = [0 0 1 1] * yheight + ystart; +fill(x, y, color, 'FaceAlpha', 0.4, 'EdgeColor', 'none'); +end + + +function c = local_colors(n) +% n distinct colors as a cell array. +base = get(groot, 'DefaultAxesColorOrder'); +if isempty(base), base = lines(7); end +while size(base, 1) < n, base = [base; base]; end %#ok +c = mat2cell(base(1:n, :), ones(n, 1), 3)'; +end + + +function local_heatmap(obj, X) +% Heat map of the full design matrix (incl. nuisance covariates / intercept). +% Each regressor (column) is scaled to unit L2 norm for display, so that +% regressors with very different magnitudes (e.g. tiny event regressors and a +% constant intercept) are shown on a comparable scale instead of the +% large-norm columns washing the others out to black. +nrm = vecnorm(X, 2, 1); +nrm(nrm == 0) = 1; +Xdisp = X ./ nrm; + +imagesc(Xdisp); +colormap(gray); +colorbar; +set(gca, 'YDir', 'reverse'); % observation 0 at the top +axis tight; +xlabel('Regressor'); +ylabel('Image / observation'); +title(sprintf('Full design matrix, columns scaled to unit norm (%d interest, %d nuisance, %d intercept)', ... + sum(obj.wh_interest), sum(obj.wh_nuisance), sum(obj.wh_intercept))); + +rn = obj.regressor_names; +if ~isempty(rn) && numel(rn) == size(X, 2) + labels = rn(:)'; + for j = 1:numel(labels) + if obj.wh_intercept(j), labels{j} = [labels{j} ' (icpt)']; + elseif obj.wh_nuisance(j), labels{j} = [labels{j} ' (nuis)']; + end + end + set(gca, 'XTick', 1:size(X, 2), 'XTickLabel', labels, 'XTickLabelRotation', 45, 'TickLabelInterpreter', 'none'); +end +end + + +function local_vif_figure(obj) +% When diagnostics have been computed, show variance inflation factors for the +% regressors of interest -- in the full design (with nuisance covariates) and, +% if nuisance indicators exist, without them -- plus contrast VIFs if present. +% Style and severity lines borrowed from scn_spm_design_check. +dg = obj.diagnostics; +if ~isstruct(dg) || ~isfield(dg, 'Variance_inflation_factors') || isempty(dg.Variance_inflation_factors) + return +end + +whI = obj.wh_interest; +rn = obj.regressor_names; +names_int = {}; +if ~isempty(rn) && numel(rn) == numel(whI), names_int = rn(whI); end + +vif_full = dg.Variance_inflation_factors(whI); % interest regressors, full design + +% Interest-only VIFs (nuisance removed), restricted to the interest regressors +vif_io = []; +has_io = any(obj.wh_nuisance) && isfield(dg, 'Variance_inflation_factors_interest_only') ... + && ~isempty(dg.Variance_inflation_factors_interest_only); +if has_io + io_cols = dg.wh_interest_only_columns; + vif_io = dg.Variance_inflation_factors_interest_only(obj.wh_interest(io_cols)); +end + +% Contrast VIFs +cvif = []; +if isfield(dg, 'Contrast_variance_inflation_factors') && ~isempty(dg.Contrast_variance_inflation_factors) + cvif = dg.Contrast_variance_inflation_factors; +end + +npanels = 1 + has_io + ~isempty(cvif); +create_figure('glm_map VIFs', 1, npanels); + +p = 1; +subplot(1, npanels, p); p = p + 1; +local_vif_panel(vif_full, names_int, 'VIF: regressors of interest (full design)'); + +if has_io + subplot(1, npanels, p); p = p + 1; + local_vif_panel(vif_io, names_int, 'VIF: of interest, nuisance removed'); +end + +if ~isempty(cvif) + subplot(1, npanels, p); + local_vif_panel(cvif, obj.contrast_names, 'Contrast VIFs (cVIF)'); +end +end + + +function local_vif_panel(vals, names, ttl) +% One VIF panel: orange markers + severity reference lines at VIF = 1/2/4/8 +% (1 = best/orthogonal; each line is a doubling of variance inflation). +vals = vals(:)'; +n = numel(vals); + +% Set x-limits first so the horizontal reference lines span the axis +plot(1:n, vals, 'ko', 'MarkerFaceColor', [1 .5 0], 'MarkerSize', 7); hold on; +xlim([0.5, n + 0.5]); + +plot_horizontal_line(1, 'k'); % minimum possible (best) +plot_horizontal_line(2, 'b--'); % doublings of variance inflation +plot_horizontal_line(4, 'r--'); +plot_horizontal_line(8, 'r-'); + +finite_vals = vals(isfinite(vals)); +ymax = 9; if ~isempty(finite_vals), ymax = max(9, 1.15 * max(finite_vals)); end +ylim([0, ymax]); + +ylabel('Variance inflation factor'); +xlabel('Predictor'); +title(ttl); + +if ~isempty(names) && numel(names) == n && n <= 20 + set(gca, 'XTick', 1:n, 'XTickLabel', names, 'XTickLabelRotation', 45, 'TickLabelInterpreter', 'none'); +else + set(gca, 'XTick', 1:n); +end +hold off; +end diff --git a/CanlabCore/@glm_map/private/select_map.m b/CanlabCore/@glm_map/private/select_map.m new file mode 100644 index 00000000..e61f95f8 --- /dev/null +++ b/CanlabCore/@glm_map/private/select_map.m @@ -0,0 +1,54 @@ +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 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 '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; + 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/replace_basis_set.m b/CanlabCore/@glm_map/replace_basis_set.m new file mode 100644 index 00000000..376ba66f --- /dev/null +++ b/CanlabCore/@glm_map/replace_basis_set.m @@ -0,0 +1,140 @@ +function obj = replace_basis_set(obj, condition_num, xBF_hires, varargin) +% Replace the basis set for one condition of an event/1st-level glm_map. +% +% Delegates to fmri_glm_design_matrix.replace_basis_set, rebuilds the design +% matrix, and clears any fitted results / contrasts / diagnostics that were +% computed with the previous basis set (they no longer apply once the number +% and meaning of the design columns change). If a data object is supplied (or +% has been retained), the model is re-fit; otherwise design diagnostics are +% recomputed on the new design. +% +% :Usage: +% :: +% +% obj = replace_basis_set(obj, condition_num, xBF_hires) +% obj = replace_basis_set(obj, condition_num, xBF_hires, 'data', fmri_data_obj, ...) +% +% :Inputs: +% +% **obj:** +% A glm_map object built around an fmri_glm_design_matrix (event / +% 1st-level mode; obj.design non-empty). +% +% **condition_num:** +% Index of the condition whose basis set to replace. +% +% **xBF_hires:** +% A basis-set struct (e.g. from fmri_spline_basis or spm_get_bf) with a +% .bf matrix and .dt sampling interval (seconds/sample). +% +% :Optional Inputs: +% +% **'data', fmri_data_obj:** +% If provided, re-fit the model on this data after rebuilding. Any +% remaining optional arguments are passed through to fit(). +% +% :Outputs: +% +% **obj:** +% The glm_map with the new basis set, a rebuilt design, stale fields +% cleared, and either re-fit (if data given) or with refreshed design +% diagnostics. +% +% :Examples: +% :: +% +% d = fmri_glm_design_matrix(2, 'nscan', 200, 'units', 'secs', ... +% 'onsets', {[10 40 70]' [25 55 85]'}, 'condition_names', {'A' 'B'}); +% g = glm_map(d); g = build_design(g); +% [xBF_hires, ~] = fmri_spline_basis(2, 'length', 20, 'nbasis', 4, 'order', 3); +% g = replace_basis_set(g, 1, xBF_hires); % condition A -> spline +% g = replace_basis_set(g, 1, xBF_hires, 'data', bold); % and re-fit +% +% :See also: +% - fmri_glm_design_matrix.replace_basis_set, build_design, fit, diagnostics +% +% .. +% 2026 - Initial implementation. +% .. + +% ------------------------------------------------------------------------- +% Validate +% ------------------------------------------------------------------------- +if isempty(obj.design) || ~isa(obj.design, 'fmri_glm_design_matrix') + error('glm_map:NoEventDesign', ... + ['replace_basis_set requires an event/1st-level glm_map with an ' ... + 'fmri_glm_design_matrix in obj.design.']); +end + +% Optional data for re-fitting (and pass-through fit options) +data = []; +wh = find(strcmpi(varargin, 'data')); +if ~isempty(wh) + data = varargin{wh(1) + 1}; + varargin(wh(1):wh(1) + 1) = []; +end + +had_contrasts = ~isempty(obj.contrasts); + +% ------------------------------------------------------------------------- +% Replace the basis set in the wrapped design and rebuild X +% ------------------------------------------------------------------------- +obj.design = replace_basis_set(obj.design, condition_num, xBF_hires); +obj.design = build(obj.design); +obj.level = 1; + +% ------------------------------------------------------------------------- +% Clear fields left over from the previous basis set that no longer apply. +% Changing the basis set changes the number and meaning of design columns, +% so fitted maps, contrasts, the fit's input structs, and diagnostics are +% all invalidated. +% ------------------------------------------------------------------------- +obj.betas = []; +obj.t = []; +obj.contrast_estimates = []; +obj.contrast_t = []; +obj.df = []; +obj.sigma = []; +obj.residuals = []; +obj.dfe = []; + +obj.contrasts = []; % defined over the old regressors +obj.contrast_names = {}; +obj.contrast_summary_table = table(); + +obj.input_parameters = struct(); +obj.input_image_metadata = struct(); +obj.fit_parameters = struct(); +obj.diagnostics = struct(); + +% Restore the full (empty) schema for the nested structs +obj = validate_object(obj); + +if had_contrasts + obj.warnings{end + 1} = ['replace_basis_set: contrasts were cleared because the ' ... + 'design columns changed with the new basis set. Re-add contrasts with add_contrasts.']; +end + +obj.history{end + 1} = sprintf('replace_basis_set: condition %d -> %s [%d x %d design]', ... + condition_num, local_bf_name(xBF_hires), size(obj.X, 1), size(obj.X, 2)); + +% ------------------------------------------------------------------------- +% Re-fit if data supplied; otherwise refresh design diagnostics +% ------------------------------------------------------------------------- +if ~isempty(data) + obj = fit(obj, data, varargin{:}); +elseif ~isempty(obj.X) + obj = run_diagnostics(obj, 'noverbose'); +end + +end % replace_basis_set + + +% ========================================================================= +function nm = local_bf_name(xBF) +if isstruct(xBF) && isfield(xBF, 'name') && ~isempty(xBF.name) + nm = xBF.name; +else + nm = 'custom basis set'; +end +end diff --git a/CanlabCore/@glm_map/run_diagnostics.m b/CanlabCore/@glm_map/run_diagnostics.m new file mode 100644 index 00000000..896e5c6c --- /dev/null +++ b/CanlabCore/@glm_map/run_diagnostics.m @@ -0,0 +1,575 @@ +function obj = run_diagnostics(obj, varargin) +% Compute and report design diagnostics for a glm_map object. +% +% Note: this method is named run_diagnostics (a verb) to avoid clashing with +% the .diagnostics property, which stores the results it computes. Call it +% with function syntax: g = run_diagnostics(g, ...). The computed results are +% stored in g.diagnostics. +% +% Evaluates the conditioning of the design matrix X (and contrasts C): +% variance inflation factors (VIF) per regressor, contrast VIFs (cVIF), +% per-observation leverage and Cook's distance, condition number, rank +% deficiency, and a redundant/near-collinear column report. VIFs/cVIFs and the +% condition number are computed both for the full design and for the +% regressors of interest only (excluding nuisance covariates), so that the +% influence of nuisance covariates on the of-interest estimates is visible. +% By default a narrative report is printed with interpretation and warnings. +% +% :Usage: +% :: +% +% obj = run_diagnostics(obj, varargin) +% +% :Inputs: +% +% **obj:** +% A glm_map object with a design matrix available (obj.X non-empty). +% Regressor roles (of interest vs nuisance vs intercept) are read from +% obj.wh_interest / obj.wh_nuisance / obj.wh_intercept. +% +% :Optional Inputs: +% +% **'doverbose' / 'noverbose':** +% Print the narrative report (default true). +% +% **'vif_thresh', [value]:** +% VIF/cVIF warning threshold (default 4). +% +% **'cond_thresh', [value]:** +% Condition-number warning threshold (default 30; > 100 = severe). +% +% :Outputs: +% +% **obj:** +% glm_map with obj.diagnostics populated. Fields include +% Variance_inflation_factors, Contrast_variance_inflation_factors, +% Leverages, Cooks_distance, condition_number, rank_deficient, +% collinearity_report (all for the FULL design), plus the +% *_interest_only counterparts, and warnings appended. +% +% :Examples: +% :: +% +% g = glm_map('X', [ones(30,1) zscore((1:30)')], 'level', 2); +% g = run_diagnostics(g); +% +% :See also: +% - VIF, cVIF, fmri_data.regress +% +% .. +% 2026 - Operates on the design only (no fit required), so it can be used to +% screen a design before fitting. Cook's distance additionally needs +% residuals and is computed only when those are present. +% .. + +% ------------------------------------------------------------------------- +% Parse inputs +% ------------------------------------------------------------------------- +doverbose = ~any(strcmpi(varargin, 'noverbose')); + +vif_thresh = 4; +wh = find(strcmpi(varargin, 'vif_thresh')); +if ~isempty(wh), vif_thresh = varargin{wh(1) + 1}; end + +cond_thresh = 30; +wh = find(strcmpi(varargin, 'cond_thresh')); +if ~isempty(wh), cond_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 = {}; + +d = obj.diagnostics; +if ~isstruct(d), d = struct(); end +d.vif_threshold = vif_thresh; +d.cond_threshold = cond_thresh; + +% Regressor roles +whI = obj.wh_interest(:)'; +whN = obj.wh_nuisance(:)'; +whB = obj.wh_intercept(:)'; +wh_io = whI | whB; % interest + intercept (model without nuisance) +has_nuisance = any(whN); + +% ------------------------------------------------------------------------- +% Variance inflation factors -- full design +% ------------------------------------------------------------------------- +d.Variance_inflation_factors = VIF(X); + +if any(d.Variance_inflation_factors > vif_thresh) + mywarnings{end + 1} = sprintf(['Design multicollinearity: %d regressor(s) have VIF > %g. ' ... + 'Check obj.diagnostics.Variance_inflation_factors and obj.regressor_names.'], ... + sum(d.Variance_inflation_factors > vif_thresh), vif_thresh); +end + +% Contrast VIFs -- full design +d.Contrast_variance_inflation_factors = []; +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 + d.Contrast_variance_inflation_factors = cVIF(X, obj.contrasts'); + if any(d.Contrast_variance_inflation_factors > vif_thresh) + mywarnings{end + 1} = sprintf('%d contrast(s) have cVIF > %g (full design).', ... + sum(d.Contrast_variance_inflation_factors > vif_thresh), vif_thresh); + end + end +end + +% ------------------------------------------------------------------------- +% Variance inflation factors -- regressors of interest only (no nuisance) +% Stored over the interest+intercept columns (see d.wh_interest_only_columns). +% ------------------------------------------------------------------------- +d.Variance_inflation_factors_interest_only = []; +d.Contrast_variance_inflation_factors_interest_only = []; +d.wh_interest_only_columns = find(wh_io); + +if has_nuisance && any(wh_io) + Xio = X(:, wh_io); + d.Variance_inflation_factors_interest_only = VIF(Xio); + + if ~isempty(obj.contrasts) && size(obj.contrasts, 1) == size(X, 2) + Cio = obj.contrasts(wh_io, :); + keepcon = any(Cio ~= 0, 1); % contrasts that survive restriction + if any(keepcon) + d.Contrast_variance_inflation_factors_interest_only = cVIF(Xio, Cio(:, keepcon)'); + end + end +end + +% ------------------------------------------------------------------------- +% Leverage (per observation) +% ------------------------------------------------------------------------- +H = X * pinv(X); +d.Leverages = diag(H)'; + +if any(abs(zscore(d.Leverages)) >= 3) + mywarnings{end + 1} = ['Some observations have high leverage (abs(z(leverage)) >= 3); ' ... + 'the fit may be unstable. Check obj.diagnostics.Leverages.']; +end + +% ------------------------------------------------------------------------- +% Cook's distance (per observation; requires residuals + sigma) +% D_i = mean_v [ r_iv^2 / (k * s_v^2) ] * h_i / (1 - h_i)^2 +% ------------------------------------------------------------------------- +d.Cooks_distance = []; +have_resid = ~isempty(obj.residuals) && isa(obj.residuals, 'fmri_data') && ~isempty(obj.residuals.dat); +have_sigma = ~isempty(obj.sigma) && isa(obj.sigma, 'fmri_data') && ~isempty(obj.sigma.dat); + +if have_resid && have_sigma + h = d.Leverages(:); + k = size(X, 2); + r = obj.residuals.dat; + s2 = obj.sigma.dat(:) .^ 2; + if size(r, 2) == numel(h) + good = isfinite(s2) & s2 > 0; + if any(good) + msr = mean( (r(good, :) .^ 2) ./ s2(good), 1 ); + hfac = (h ./ (1 - h) .^ 2)'; + d.Cooks_distance = (msr / k) .* hfac; + if any(d.Cooks_distance > 1) + mywarnings{end + 1} = sprintf(['%d observation(s) have Cook''s distance > 1 ' ... + '(highly influential). Check obj.diagnostics.Cooks_distance.'], sum(d.Cooks_distance > 1)); + end + end + end +end + +% ------------------------------------------------------------------------- +% Conditioning / rank (full and interest-only) +% +% The condition number is computed after scaling each column to unit L2 norm +% (the Belsley scaled condition index). This makes it scale-invariant -- like +% VIF -- so it reflects true collinearity rather than differences in regressor +% magnitude. Without scaling, a constant intercept (large norm) alongside +% small-amplitude event regressors produces a large condition number that is a +% pure scaling artifact, even when the regressors are nearly orthogonal. +% ------------------------------------------------------------------------- +d.condition_number = local_scaled_cond(X); +d.condition_number_interest_only = []; +if has_nuisance && any(wh_io), d.condition_number_interest_only = local_scaled_cond(X(:, wh_io)); end +d.rank_deficient = rank(X) < size(X, 2); + +if d.rank_deficient + mywarnings{end + 1} = 'Design matrix X is rank deficient (rank(X) < number of regressors).'; +end +if d.condition_number > 100 + mywarnings{end + 1} = sprintf('Design is severely ill-conditioned (condition number %.0f > 100).', d.condition_number); +elseif d.condition_number > cond_thresh + mywarnings{end + 1} = sprintf('Design is moderately ill-conditioned (condition number %.0f > %g).', d.condition_number, cond_thresh); +end + +% Does adding nuisance covariates substantially inflate of-interest VIFs? +infl = local_inflation(d, whI, wh_io); +if ~isempty(infl) && infl.max_ratio >= 1.5 + mywarnings{end + 1} = sprintf(['Nuisance covariates inflate of-interest VIFs by up to %.1fx ' ... + '(max VIF %.2f vs %.2f without nuisance); events of interest may correlate with nuisance variables.'], ... + infl.max_ratio, infl.max_full, infl.max_io); +end + +% ------------------------------------------------------------------------- +% Design efficiency for contrasts (calcEfficiency). Uses the entered +% contrasts, or an orthogonal set spanning the regressors of interest when +% none are entered. Higher efficiency = lower contrast-estimate variance. +% ------------------------------------------------------------------------- +[d.efficiency, d.efficiency_per_contrast, d.efficiency_contrast_names, ... + d.efficiency_contrast_source, eff_note] = local_efficiency(obj, X, whI); + +% ------------------------------------------------------------------------- +% Temporal frequency content and high-pass filter recommendation. Only +% meaningful for within-run timeseries designs (obj.is_timeseries). For each +% regressor of interest (and each contrast), compute the cumulative power by +% frequency and the high-pass filter cutoff that removes < 5% of its variance; +% recommend the cutoff that keeps the MAX loss across regressors (and across +% contrasts) below 5%. Borrows the approach from scn_spm_choose_hpfilter. +% ------------------------------------------------------------------------- +d.hpfilter = []; +if obj.is_timeseries + d.hpfilter = local_hpfilter(obj, X, whI, whB); +end + +% ------------------------------------------------------------------------- +% Redundant / near-collinear column report +% ------------------------------------------------------------------------- +report = struct(); +report.vif_threshold = vif_thresh; +report.high_vif_columns = find(d.Variance_inflation_factors > vif_thresh); + +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]; end %#ok + 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.diagnostics.collinearity_report).', size(dup_pairs, 1)); +end + +R = corrcoef(X); +R(logical(eye(ncol))) = 0; +[ia, ib] = find(triu(abs(R) > 0.95, 1)); +report.high_correlation_pairs = [ia ib]; + +d.collinearity_report = report; + +% ------------------------------------------------------------------------- +% Store +% ------------------------------------------------------------------------- +obj.diagnostics = d; +obj.warnings = [obj.warnings(:); mywarnings(:)]'; +obj.history{end + 1} = 'diagnostics: VIF/cVIF (full + interest-only), leverage, Cook''s D, conditioning, collinearity'; + +% ------------------------------------------------------------------------- +% Report +% ------------------------------------------------------------------------- +if doverbose + local_report(obj, d, infl, vif_thresh, cond_thresh, mywarnings, eff_note); +end + +end % diagnostics + + +% ========================================================================= +% Local helpers +% ========================================================================= +function [eff, eff_vec, cnames, src, note] = local_efficiency(obj, X, whI) +% Design efficiency for contrasts via calcEfficiency. Uses entered contrasts, +% else an orthogonal set spanning the regressors of interest. Returns [] for +% eff and an informative note when efficiency cannot be computed. +[eff, eff_vec, cnames, src, note] = deal([], [], {}, '', ''); + +Crows = []; +if ~isempty(obj.contrasts) && size(obj.contrasts, 1) == size(X, 2) + Crows = obj.contrasts'; % [n_con x n_reg] + cnames = obj.contrast_names; + src = 'entered contrasts'; + +elseif any(whI) + whIidx = find(whI); + if numel(whIidx) >= 2 + Cint = create_orthogonal_contrast_set(numel(whIidx)); % dispatches to the function + Crows = zeros(size(Cint, 1), size(X, 2)); + Crows(:, whIidx) = Cint; + cnames = arrayfun(@(i) sprintf('OrthC%d', i), 1:size(Crows, 1), 'UniformOutput', false); + src = 'auto orthogonal set spanning the regressors of interest'; + else + note = 'Efficiency not computed: need >= 2 regressors of interest, or enter contrasts.'; + return + end +else + note = ['Efficiency not computed: no regressors of interest are defined and no contrasts ' ... + 'entered. Define regressors of interest (build an event design, or set ' ... + 'obj.nuisance_columns for direct designs) or add contrasts.']; + return +end + +try + [eff, eff_vec] = calcEfficiency(ones(1, size(Crows, 1)), Crows, pinv(X), []); + eff_vec = eff_vec(:)'; +catch ME + [eff, eff_vec, cnames, src] = deal([], [], {}, ''); + note = sprintf('Efficiency not computed (calcEfficiency error): %s', ME.message); +end +end + + +function c = local_scaled_cond(X) +% Condition number after scaling each column to unit L2 norm (scale-invariant +% Belsley condition index). +nrm = vecnorm(X, 2, 1); +nrm(nrm == 0) = 1; +c = cond(X ./ nrm); +end + + +function hp = local_hpfilter(obj, X, whI, whB) +% Cumulative power by frequency and recommended high-pass filter cutoff for a +% timeseries design. Returns [] if no interest regressors or no valid TR. +hp = []; + +TR = obj.TR; +if isempty(TR) || ~isscalar(TR) || isnan(TR) || TR <= 0, return, end +if ~any(whI), return, end + +% Regressors of interest, with the intercept/baseline projected out +Xint = X(:, whI); +Xb = X(:, whB); +if ~isempty(Xb), Xint = Xint - Xb * pinv(Xb) * Xint; end + +[cp_reg, fr_reg, hplen_reg, maxloss_reg] = local_cumulative_power(Xint, TR); + +% Contrast time courses (also intercept-projected) +hplen_con = []; maxloss_con = []; cp_con = []; fr_con = []; +if ~isempty(obj.contrasts) && size(obj.contrasts, 1) == size(X, 2) + Xcon = X * obj.contrasts; + if ~isempty(Xb), Xcon = Xcon - Xb * pinv(Xb) * Xcon; end + [cp_con, fr_con, hplen_con, maxloss_con] = local_cumulative_power(Xcon, TR); +end + +hp = struct( ... + 'TR', TR, ... + 'recommended_hpfilter_sec_regressors', hplen_reg, ... + 'recommended_hpfilter_sec_contrasts', hplen_con, ... + 'max_variance_loss_regressors', maxloss_reg, ... % fraction lost at the recommended cutoff + 'max_variance_loss_contrasts', maxloss_con, ... + 'cumulative_power_regressors', cp_reg, ... + 'frequencies_regressors', fr_reg, ... + 'cumulative_power_contrasts', cp_con, ... + 'frequencies_contrasts', fr_con); +end + + +function [cumpower, freqs, hplen, maxloss] = local_cumulative_power(X, TR) +% For each column of X, cumulative (low-to-high frequency) power fraction. +% Returns the high-pass filter length (sec) whose cutoff removes < 5% of the +% variance of every column, and the actual max fraction removed at that cutoff. +% (Approach borrowed from scn_spm_choose_hpfilter / cumulative_power.) +ncol = size(X, 2); +timepts = floor(size(X, 1) / 2); +cumpower = zeros(timepts, ncol); +freqs = zeros(timepts, ncol); +loss_freq = zeros(1, ncol); +loss_at_cut = zeros(1, ncol); + +for i = 1:ncol + [mag, freq] = fft_calc(X(:, i), TR, 'noplot'); % normalized power (sums to 1), freq in Hz + cp = cumsum(mag(:)); + cumpower(:, i) = cp; + freqs(:, i) = freq(:); + + % Highest frequency at which cumulative (removed) power is still < 5% + wh = find(cp < 0.05); + if isempty(wh), wh = 1; end % even the lowest bin loses >= 5% + wh = wh(end); + loss_freq(i) = freq(wh); + loss_at_cut(i) = cp(wh); +end + +% Use the lowest cutoff frequency (longest period) so that all columns lose +% < 5%. Guard the DC bin (freq 0). +posfreq = loss_freq(loss_freq > 0); +if isempty(posfreq) + hplen = Inf; % cannot filter without > 5% loss +else + hplen = 1 / min(posfreq); +end +maxloss = max(loss_at_cut); +end + + +function infl = local_inflation(d, whI, wh_io) +% Compare full-design VIFs of the of-interest regressors against their +% interest-only VIFs. Returns a struct with the max inflation ratio, or []. +infl = []; +if isempty(d.Variance_inflation_factors_interest_only), return, end + +vfull = d.Variance_inflation_factors; +vio = d.Variance_inflation_factors_interest_only; % over wh_io columns +io_cols = find(wh_io); + +% map each interest column to its position within the interest-only design +isI_in_io = whI(io_cols); % which io columns are of-interest +full_on_I = vfull(io_cols(isI_in_io)); +io_on_I = vio(isI_in_io); + +valid = io_on_I > 0 & isfinite(io_on_I) & isfinite(full_on_I); +if ~any(valid), return, end + +ratios = full_on_I(valid) ./ io_on_I(valid); +[infl.max_ratio, w] = max(ratios); +fv = full_on_I(valid); iv = io_on_I(valid); +infl.max_full = fv(w); +infl.max_io = iv(w); +end + + +function local_report(obj, d, infl, vif_thresh, cond_thresh, mywarnings, eff_note) + +line = repmat('-', 1, 70); +fprintf('\n glm_map design diagnostics\n %s\n', line); +if ~isempty(obj.analysis_name), fprintf(' Analysis: %s\n', obj.analysis_name); end +fprintf(' %d images x %d regressors (%d of interest, %d nuisance, %d intercept)\n', ... + obj.num_images, obj.num_regressors, sum(obj.wh_interest), sum(obj.wh_nuisance), sum(obj.wh_intercept)); + +% ---- VIFs (full design) ---- +rn = obj.regressor_names; +fprintf('\n Variance inflation factors (VIF), full design\n'); +fprintf(' VIF >= 1 (1 = orthogonal). ~1-2 low, 2-4 moderate, >%g high collinearity, Inf = exact dependence.\n', vif_thresh); +fprintf(' %-28s %8s role\n', 'Regressor', 'VIF'); +for i = 1:numel(d.Variance_inflation_factors) + nm = local_name(rn, i, 'R'); + flag = local_flag(d.Variance_inflation_factors(i) > vif_thresh, ' <-- high'); + fprintf(' %-28s %8.2f %s%s\n', nm, d.Variance_inflation_factors(i), local_role(obj, i), flag); +end + +% ---- Interest-only comparison ---- +if ~isempty(d.Variance_inflation_factors_interest_only) + fprintf('\n VIF, regressors of interest only (nuisance covariates removed)\n'); + io_cols = d.wh_interest_only_columns; + vio = d.Variance_inflation_factors_interest_only; + for j = 1:numel(io_cols) + if obj.wh_intercept(io_cols(j)), continue, end % skip intercept in this listing + nm = local_name(rn, io_cols(j), 'R'); + fprintf(' %-28s %8.2f (full: %.2f)\n', nm, vio(j), d.Variance_inflation_factors(io_cols(j))); + end + if ~isempty(infl) + if infl.max_ratio >= 1.5 + fprintf(' => Nuisance covariates inflate of-interest VIFs up to %.1fx (%.2f vs %.2f).\n', ... + infl.max_ratio, infl.max_full, infl.max_io); + fprintf(' This suggests events of interest are correlated with nuisance variables.\n'); + else + fprintf(' => Nuisance covariates have little effect on of-interest VIFs (max %.1fx).\n', infl.max_ratio); + end + end +end + +% ---- Contrast VIFs ---- +if ~isempty(d.Contrast_variance_inflation_factors) + fprintf('\n Contrast VIFs (cVIF), full design [same scale as VIF]\n'); + for i = 1:numel(d.Contrast_variance_inflation_factors) + nm = local_name(obj.contrast_names, i, 'Con'); + flag = local_flag(d.Contrast_variance_inflation_factors(i) > vif_thresh, ' <-- high'); + fprintf(' %-28s %8.2f%s\n', nm, d.Contrast_variance_inflation_factors(i), flag); + end + if ~isempty(d.Contrast_variance_inflation_factors_interest_only) + fprintf(' Interest-only cVIF: %s\n', mat2str(round(d.Contrast_variance_inflation_factors_interest_only, 2))); + end +end + +% ---- Conditioning ---- +fprintf('\n Conditioning\n'); +fprintf(' Scaled condition number (columns at unit norm; scale-invariant, like VIF).\n'); +fprintf(' >= 1; <%g OK, %g-100 moderate, >100 severe multicollinearity.\n', cond_thresh, cond_thresh); +fprintf(' condition number (full) : %.1f%s\n', d.condition_number, ... + local_flag(d.condition_number > cond_thresh, ' <-- elevated')); +if ~isempty(d.condition_number_interest_only) + fprintf(' condition number (interest only): %.1f\n', d.condition_number_interest_only); +end +fprintf(' rank deficient : %s\n', local_yn(d.rank_deficient)); + +% ---- Cook's distance ---- +if ~isempty(d.Cooks_distance) + [maxcd, wcd] = max(d.Cooks_distance); + fprintf('\n Influence: Cook''s distance per observation (image)\n'); + fprintf(' > 1 = highly influential (rule of thumb; 4/n = %.3f is a stricter cutoff).\n', 4 / obj.num_images); + fprintf(' max Cook''s distance: %.3f (observation %d)%s\n', maxcd, wcd, local_flag(maxcd > 1, ' <-- influential')); +end + +% ---- Design efficiency ---- +fprintf('\n Design efficiency for contrasts (calcEfficiency)\n'); +fprintf(' Relative measure (no absolute scale); higher = lower contrast-estimate\n'); +fprintf(' variance = more efficient. Compare designs/contrasts on the same data.\n'); +if ~isempty(d.efficiency) + fprintf(' Contrasts: %s\n', d.efficiency_contrast_source); + for i = 1:numel(d.efficiency_per_contrast) + nm = local_name(d.efficiency_contrast_names, i, 'Con'); + fprintf(' %-28s %12.4g\n', nm, d.efficiency_per_contrast(i)); + end + fprintf(' %-28s %12.4g\n', 'overall (mean)', d.efficiency); +else + fprintf(' %s\n', eff_note); +end + +% ---- High-pass filter (timeseries designs only) ---- +if obj.is_timeseries && isstruct(d.hpfilter) && ~isempty(fieldnames(d.hpfilter)) + hp = d.hpfilter; + fprintf('\n Temporal frequency / high-pass filter (timeseries design)\n'); + fprintf(' Recommended HP filter length keeps low-frequency variance loss < 5%% for every\n'); + fprintf(' regressor (and contrast). A *longer* length filters less; use the larger value\n'); + fprintf(' to protect both, then set it as the SPM/first-level high-pass cutoff.\n'); + fprintf(' regressors of interest : %s\n', ... + local_hpline(hp.recommended_hpfilter_sec_regressors, hp.max_variance_loss_regressors)); + if ~isempty(hp.recommended_hpfilter_sec_contrasts) + fprintf(' contrasts : %s\n', ... + local_hpline(hp.recommended_hpfilter_sec_contrasts, hp.max_variance_loss_contrasts)); + end + rec = max([hp.recommended_hpfilter_sec_regressors, hp.recommended_hpfilter_sec_contrasts]); + if isfinite(rec) + fprintf(' => suggested cutoff : %.0f s (the more conservative of the two)\n', rec); + else + fprintf(' => a high-pass filter would remove > 5%% of variance; consider not high-pass filtering.\n'); + end +end + +% ---- Warnings ---- +if ~isempty(mywarnings) + fprintf('\n %d warning(s):\n', numel(mywarnings)); + for i = 1:numel(mywarnings), fprintf(' - %s\n', mywarnings{i}); end +else + fprintf('\n No diagnostic warnings.\n'); +end +fprintf(' %s\n\n', line); + +end + + +function r = local_role(obj, i) +if obj.wh_intercept(i), r = 'intercept'; +elseif obj.wh_nuisance(i), r = 'nuisance'; +else, r = 'of interest'; +end +end + +function nm = local_name(names, idx, prefix) +if iscell(names) && idx <= numel(names) && ~isempty(names{idx}), nm = names{idx}; +else, nm = sprintf('%s%d', prefix, idx); +end +end + +function s = local_flag(tf, str), if tf, s = str; else, s = ''; end, end +function s = local_yn(tf), if tf, s = 'yes'; else, s = 'no'; end, end +function s = local_hpline(len, maxloss) +% Format one high-pass-filter recommendation line. +if isempty(len) || ~isfinite(len) + s = 'no high-pass filter recommended (even long filters would lose > 5%)'; +else + if isempty(maxloss), maxloss = 0; end + s = sprintf('%.0f s (loses %.1f%% at this cutoff)', len, 100 * maxloss); +end +end diff --git a/CanlabCore/@glm_map/summary.m b/CanlabCore/@glm_map/summary.m new file mode 100644 index 00000000..c49e612a --- /dev/null +++ b/CanlabCore/@glm_map/summary.m @@ -0,0 +1,258 @@ +function summary(obj) +% Print a narrative summary report for a glm_map object. +% +% Unlike disp/display (which list the object's properties), summary prints a +% human-readable report: the analysis name, a description of the model +% (level and input variables), the design diagnostics (VIF, conditioning, +% leverage, Cook's distance, collinearity), and -- if the model has been +% fitted -- a summary of the results, including how many result maps there +% are, the statistical threshold, and the number of significant voxels at +% that threshold for each predictor and contrast. +% +% :Usage: +% :: +% +% summary(obj) +% +% :Inputs: +% +% **obj:** +% A glm_map object (fitted or unfitted). +% +% :See also: +% - glm_map, glm_map.diagnostics, glm_map.disp +% +% .. +% Programmers' notes: +% 2026 - Rewritten as a narrative report (was a thin wrapper over disp). +% .. + +line = repmat('=', 1, 68); +fprintf('\n%s\n glm_map analysis summary\n%s\n', line, line); + +% --------------------------------------------------------------------- +% Analysis name +% --------------------------------------------------------------------- +if ~isempty(obj.analysis_name) + fprintf(' Analysis: %s\n\n', obj.analysis_name); +end + +% --------------------------------------------------------------------- +% Model description (level and input variables) +% --------------------------------------------------------------------- +switch obj.level + case 1, levelstr = 'first-level (within-run / single-subject)'; + case 2, levelstr = 'second-level (group)'; + otherwise, levelstr = num2str(obj.level); +end + +if ~isempty(obj.design) + modestr = 'event/1st-level (HRF-convolved design)'; +else + modestr = 'direct design matrix'; +end + +fprintf(' Model\n %s\n', repmat('-', 1, 64)); +fprintf(' Level : %s\n', levelstr); +fprintf(' Design : %s', modestr); +if obj.is_timeseries, fprintf(' [timeseries; AR errors allowed]'); end +fprintf('\n'); +fprintf(' Observations : %d images\n', obj.num_images); +fprintf(' Regressors : %d (%d of interest, %d nuisance, %d intercept)\n', ... + obj.num_regressors, sum(obj.wh_interest), sum(obj.wh_nuisance), sum(obj.wh_intercept)); + +rn = obj.regressor_names; +if ~isempty(rn) + for i = 1:numel(rn) + nm = rn{i}; if isempty(nm), nm = sprintf('R%d', i); end + fprintf(' %2d. %-26s [%s]\n', i, nm, local_role(obj, i)); + end +end + +if obj.num_contrasts > 0 + fprintf(' Contrasts : %d\n', obj.num_contrasts); + cn = obj.contrast_names; + for i = 1:obj.num_contrasts + if i <= numel(cn) && ~isempty(cn{i}), nm = cn{i}; else, nm = sprintf('Con%d', i); end + fprintf(' %2d. %s\n', i, nm); + end +end +fprintf('\n'); + +% --------------------------------------------------------------------- +% Design diagnostics +% --------------------------------------------------------------------- +d = obj.diagnostics; +if isstruct(d) && ~isempty(fieldnames(d)) + fprintf(' Design diagnostics\n %s\n', repmat('-', 1, 64)); + + if isfield(d, 'Variance_inflation_factors') && ~isempty(d.Variance_inflation_factors) + [maxvif, wv] = max(d.Variance_inflation_factors); + nm = local_name(rn, wv, 'R'); + fprintf(' Max VIF : %.2f (%s)%s\n', maxvif, nm, local_flag(maxvif > 4, ' high')); + end + + if isfield(d, 'Contrast_variance_inflation_factors') && ~isempty(d.Contrast_variance_inflation_factors) + [maxcv, wc] = max(d.Contrast_variance_inflation_factors); + nm = local_name(obj.contrast_names, wc, 'Con'); + fprintf(' Max contrast VIF: %.2f (%s)%s\n', maxcv, nm, local_flag(maxcv > 4, ' high')); + end + + if isfield(d, 'condition_number') && ~isempty(d.condition_number) + fprintf(' Condition number: %.2f\n', d.condition_number); + end + + if isfield(d, 'rank_deficient') && ~isempty(d.rank_deficient) + fprintf(' Rank deficient : %s\n', local_yn(d.rank_deficient)); + end + + if isfield(d, 'Leverages') && ~isempty(d.Leverages) + nhi = sum(abs(zscore(d.Leverages(:))) >= 3); + fprintf(' High leverage : %d observation(s) with abs(z) >= 3\n', nhi); + end + + if isfield(d, 'Cooks_distance') && ~isempty(d.Cooks_distance) + [maxcd, wcd] = max(d.Cooks_distance); + fprintf(' Max Cook''s dist : %.3f (observation %d)%s\n', maxcd, wcd, local_flag(maxcd > 1, ' influential')); + end + + if isfield(d, 'collinearity_report') && isstruct(d.collinearity_report) + cr = d.collinearity_report; + ndup = 0; if isfield(cr, 'duplicate_column_pairs'), ndup = size(cr.duplicate_column_pairs, 1); end + ncor = 0; if isfield(cr, 'high_correlation_pairs'), ncor = size(cr.high_correlation_pairs, 1); end + if ndup > 0 || ncor > 0 + fprintf(' Collinearity : %d duplicate column pair(s), %d high-correlation pair(s)\n', ndup, ncor); + end + end + fprintf('\n'); +end + +% --------------------------------------------------------------------- +% Results (if fitted) +% --------------------------------------------------------------------- +if obj.is_fitted + fprintf(' Results\n %s\n', repmat('-', 1, 64)); + + % How many maps + fprintf(' Maps : betas [%d], t [%d]', local_nimg(obj.betas), local_nimg(obj.t)); + if ~isempty(obj.contrast_estimates) + fprintf(', contrast estimates [%d], contrast t [%d]', ... + local_nimg(obj.contrast_estimates), local_nimg(obj.contrast_t)); + end + fprintf('\n'); + fprintf(' Error df : %g\n', obj.dfe); + + % Statistical threshold (from the thresholded t map) + fprintf(' Threshold : %s\n', local_threshold_str(obj.t, obj.fit_parameters)); + + % Significant voxels per predictor (t map) + fprintf(' Significant voxels at threshold (per regressor):\n'); + local_print_sig(obj.t, rn, 'R'); + + % Significant voxels per contrast (contrast t map) + if ~isempty(obj.contrast_t) + fprintf(' Significant voxels at threshold (per contrast):\n'); + local_print_sig(obj.contrast_t, obj.contrast_names, 'Con'); + end + fprintf('\n'); +else + fprintf(' Results : not fitted (run fit(obj, fmri_data_obj))\n\n'); +end + +% --------------------------------------------------------------------- +% Warnings +% --------------------------------------------------------------------- +if ~isempty(obj.warnings) + fprintf(' Warnings (%d)\n %s\n', numel(obj.warnings), repmat('-', 1, 64)); + for i = 1:numel(obj.warnings) + fprintf(' - %s\n', obj.warnings{i}); + end + fprintf('\n'); +end + +fprintf('%s\n\n', line); + +end % summary + + +% ===================================================================== +% Local helpers +% ===================================================================== +function r = local_role(obj, i) +if obj.wh_intercept(i), r = 'intercept'; +elseif obj.wh_nuisance(i), r = 'nuisance'; +else, r = 'of interest'; +end +end + +function nm = local_name(names, idx, prefix) +if iscell(names) && idx <= numel(names) && ~isempty(names{idx}) + nm = names{idx}; +else + nm = sprintf('%s%d', prefix, idx); +end +end + + +function s = local_flag(tf, str) +if tf, s = str; else, s = ''; end +end + + +function s = local_yn(tf) +if tf, s = 'yes'; else, s = 'no'; end +end + + +function n = local_nimg(map) +if isempty(map) || ~isprop(map, 'dat') || isempty(map.dat), n = 0; else, n = size(map.dat, 2); end +end + + +function s = local_threshold_str(map, fit_parameters) +% Build a readable threshold description. Prefer the glm_map fit_parameters +% (clean p-value/type), then fall back to the statistic_image's own record. +s = '(unthresholded)'; + +if isstruct(fit_parameters) && isfield(fit_parameters, 'pthresh') && ~isempty(fit_parameters.pthresh) + ttype = ''; if isfield(fit_parameters, 'thresh_type'), ttype = char(fit_parameters.thresh_type); end + s = strtrim(sprintf('p < %g %s', fit_parameters.pthresh, ttype)); + +elseif ~isempty(map) && isprop(map, 'threshold') && ~isempty(map.threshold) + thr = unique(map.threshold(:)'); % collapse duplicated pos/neg thresholds + ttype = ''; + if isprop(map, 'thr_type') && ~isempty(map.thr_type) + if iscell(map.thr_type), ttype = strjoin(cellfun(@char, map.thr_type(:)', 'UniformOutput', false), '/'); + else, ttype = char(map.thr_type); end + end + if isscalar(thr) + s = strtrim(sprintf('p < %g %s', thr, ttype)); + else + s = strtrim(sprintf('%s %s', mat2str(thr, 4), ttype)); + end +end +end + + +function local_print_sig(map, names, prefix) +% Print the number of significant voxels per image of a thresholded +% statistic_image map. +if isempty(map) || ~isprop(map, 'dat') || isempty(map.dat) + fprintf(' (no map)\n'); + return +end + +nimg = size(map.dat, 2); + +% Prefer the .sig logical mask; fall back to nonzero thresholded values +if isprop(map, 'sig') && ~isempty(map.sig) && size(map.sig, 2) == nimg + counts = sum(map.sig ~= 0, 1); +else + counts = sum(map.dat ~= 0 & ~isnan(map.dat), 1); +end + +for i = 1:nimg + nm = local_name(names, i, prefix); + fprintf(' %-26s %d voxels\n', nm, counts(i)); +end +end diff --git a/CanlabCore/@glm_map/table.m b/CanlabCore/@glm_map/table.m new file mode 100644 index 00000000..73ecd549 --- /dev/null +++ b/CanlabCore/@glm_map/table.m @@ -0,0 +1,67 @@ +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, 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 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 diff --git a/CanlabCore/@glm_map/validate_object.m b/CanlabCore/@glm_map/validate_object.m new file mode 100644 index 00000000..bb5347bc --- /dev/null +++ b/CanlabCore/@glm_map/validate_object.m @@ -0,0 +1,151 @@ +function obj = validate_object(obj) +% Ensure a glm_map's nested structs expose their full set of valid fields. +% +% glm_map groups related outputs into nested structs (.input_parameters, +% .input_image_metadata, .diagnostics, and .fit_parameters). Depending on how +% an object was built (empty constructor, re-cast from an fmri_data.regress +% results struct, or partially populated), some of those structs may be empty +% or missing fields. validate_object fills in any missing fields with empty +% defaults so that every struct always exposes the complete schema, and +% normalizes the nested diagnostics.collinearity_report the same way. +% +% Existing field values are never overwritten; only absent fields are added. +% A struct property that is not a struct (e.g. left as []) is replaced by the +% empty-field template. The field templates are the single source of truth for +% the valid field names produced by fmri_data.regress and glm_map.diagnostics. +% +% :Usage: +% :: +% +% obj = validate_object(obj) +% +% :Inputs: +% +% **obj:** +% A glm_map object. +% +% :Outputs: +% +% **obj:** +% The glm_map object with .input_parameters, .input_image_metadata, +% .diagnostics (and .diagnostics.collinearity_report), and +% .fit_parameters guaranteed to contain every valid field. +% +% :Examples: +% :: +% +% g = validate_object(glm_map); % empty object, full schema +% fieldnames(g.diagnostics) % all diagnostic fields present +% +% :See also: +% - glm_map, glm_map.check_properties, fmri_data.regress +% +% .. +% Programmers' notes: +% 2026 - Initial implementation. Templates below define the valid field +% set for each nested struct; keep them in sync with the fields produced by +% fmri_data.regress (out.input_parameters / out.input_image_metadata / +% out.diagnostics) and glm_map.diagnostics. +% .. + +obj.input_parameters = local_fill(obj.input_parameters, local_tmpl_input_parameters()); +obj.input_image_metadata = local_fill(obj.input_image_metadata, local_tmpl_input_image_metadata()); +obj.diagnostics = local_fill(obj.diagnostics, local_tmpl_diagnostics()); +obj.fit_parameters = local_fill(obj.fit_parameters, local_tmpl_fit_parameters()); + +% Nested struct inside .diagnostics +obj.diagnostics.collinearity_report = local_fill(obj.diagnostics.collinearity_report, ... + local_tmpl_collinearity_report()); + +end % validate_object + + +% ===================================================================== +% Field-filling helper +% ===================================================================== +function s = local_fill(s, tmpl) +% Add every field of tmpl that is missing from s (default = template value). +% Existing fields are preserved. A non-struct (or non-scalar struct) s is +% replaced by the template. +if ~isstruct(s) || ~isscalar(s) + s = tmpl; + return +end + +fn = fieldnames(tmpl); +for i = 1:numel(fn) + if ~isfield(s, fn{i}) + s.(fn{i}) = tmpl.(fn{i}); + end +end +end + + +% ===================================================================== +% Field templates (valid field sets). Defaults are [] = "not set / not +% computed". Keep in sync with fmri_data.regress and glm_map.diagnostics. +% ===================================================================== +function t = local_tmpl_input_parameters() +t = struct( ... + 'brain_is_predictor', [], ... + 'do_robust', [], ... + 'grandmeanscale', [], ... + 'do_intercept', [], ... + 'do_resid', [], ... + 'doverbose', [], ... + 'do_display', [], ... + 'covdat', [], ... + 'initial_statistical_threshold', []); +end + + +function t = local_tmpl_input_image_metadata() +t = struct( ... + 'source_notes', [], ... + 'history', [], ... + 'image_names', [], ... + 'fullpath', []); +end + + +function t = local_tmpl_diagnostics() +t = struct( ... + 'Variance_inflation_factors', [], ... + 'Contrast_variance_inflation_factors', [], ... + 'Variance_inflation_factors_interest_only', [], ... + 'Contrast_variance_inflation_factors_interest_only',[], ... + 'wh_interest_only_columns', [], ... + 'Leverages', [], ... + 'Cooks_distance', [], ... + 'condition_number', [], ... + 'condition_number_interest_only', [], ... + 'rank_deficient', [], ... + 'collinearity_report', [], ... + 'efficiency', [], ... + 'efficiency_per_contrast', [], ... + 'efficiency_contrast_names', [], ... + 'efficiency_contrast_source', [], ... + 'hpfilter', [], ... + 'vif_threshold', [], ... + 'cond_threshold', []); +end + + +function t = local_tmpl_collinearity_report() +t = struct( ... + 'vif_threshold', [], ... + 'high_vif_columns', [], ... + 'duplicate_column_pairs', [], ... + 'high_correlation_pairs', []); +end + + +function t = local_tmpl_fit_parameters() +t = struct( ... + 'robust', [], ... + 'ar_order', [], ... + 'is_timeseries', [], ... + 'pthresh', [], ... + 'thresh_type', [], ... + 'do_resid', []); +end 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: diff --git a/CanlabCore/Data_processing_tools/smooth_timeseries.m b/CanlabCore/Data_processing_tools/smooth_timeseries.m index 3bb442d7..eb5cf883 100644 --- a/CanlabCore/Data_processing_tools/smooth_timeseries.m +++ b/CanlabCore/Data_processing_tools/smooth_timeseries.m @@ -23,7 +23,7 @@ xc = xc - min(xc); %xc = xc ./ sum(xc); irrelevant -V = toeplitz(pad(xc',size(x,1)-length(xc))); +V = toeplitz(canlab_pad(xc',size(x,1)-length(xc))); % canlab_pad: renamed from pad.m (avoids builtin pad) V = V ./ repmat(sum(V),size(V,1),1); V = V'; diff --git a/CanlabCore/Model_building_tools/onsets2fmridesign.m b/CanlabCore/Model_building_tools/onsets2fmridesign.m index dd0a230d..dffcc35b 100644 --- a/CanlabCore/Model_building_tools/onsets2fmridesign.m +++ b/CanlabCore/Model_building_tools/onsets2fmridesign.m @@ -533,12 +533,14 @@ if ~isempty(varargin) for i = 1:length(ons) - try - % 2019 internal matlab function - delta{i} = padarray(delta{i}, ceil(len./res./TR - length(delta{i}))); - catch - % This is in CanlabCore/Misc_utilities - delta{i} = pad(delta{i}, ceil(len./res./TR - length(delta{i}))); + % Zero-pad the TR-resolution onset indicator to the run length. + % (Was: padarray, falling back to CanlabCore's pad.m. pad.m was + % renamed to canlab_pad.m and 'pad' now resolves to MATLAB's builtin + % string pad -- which errors on numeric input in R2024b+ -- so append + % zeros explicitly, matching canlab_pad's one-sided behavior.) + npad = ceil(len ./ res ./ TR - length(delta{i})); + if npad > 0 + delta{i} = [delta{i}; zeros(npad, 1)]; end end diff --git a/CanlabCore/Statistics_tools/nonlin_parammod_predfun.m b/CanlabCore/Statistics_tools/nonlin_parammod_predfun.m index 2e6ef758..32714962 100644 --- a/CanlabCore/Statistics_tools/nonlin_parammod_predfun.m +++ b/CanlabCore/Statistics_tools/nonlin_parammod_predfun.m @@ -183,7 +183,9 @@ yhi = p(5) + fast_conv_fft(hrf, xbox); if length(yhi) < runDuration - yhi = pad(yhi, runDuration); + % pad yhi out to runDuration (canlab_pad takes a count of zeros to append; + % 'pad' was renamed to canlab_pad and now collides with the builtin pad) + yhi = canlab_pad(yhi, runDuration - length(yhi)); elseif length(yhi) > runDuration yhi = yhi(1:runDuration); end @@ -254,7 +256,7 @@ hrf = hrf/max(hrf); end - hrf = pad(hrf,x); + hrf = canlab_pad(hrf, x); % pad hrf to length(x) (renamed from pad.m) % Check if length(pm_vals) ~= length(ons) 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..9d0d8266 --- /dev/null +++ b/CanlabCore/Unit_tests/glm_map/canlab_test_glm_map.m @@ -0,0 +1,500 @@ +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 = run_diagnostics(g, 'noverbose'); + +tc.verifyNumElements(g.diagnostics.Variance_inflation_factors, 3); +tc.verifyNumElements(g.diagnostics.Contrast_variance_inflation_factors, 1); +tc.verifyGreaterThan(g.diagnostics.condition_number, 0); +tc.verifyFalse(g.diagnostics.rank_deficient); +tc.verifyTrue(isstruct(g.diagnostics.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 = run_diagnostics(g, 'noverbose'); + +tc.verifyTrue(g.diagnostics.rank_deficient); +tc.verifySize(g.diagnostics.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.diagnostics.Variance_inflation_factors); +tc.verifyClass(g.df, 'fmri_data'); % per-voxel df kept as fmri_data +tc.verifyClass(g.sigma, 'fmri_data'); % per-voxel sigma kept as fmri_data +tc.verifyTrue(isstruct(g.input_parameters)); % nested option struct populated +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 + + +% ===================================================================== +% fmri_data.regress now returns a glm_map; constructor re-casts a struct +% ===================================================================== +function test_regress_returns_glm_map(tc) +dat = canlab_get_sample_fmri_data(); +n = size(dat.dat, 2); +dat.X = [zscore((1:n)') ones(n, 1)]; +out = regress(dat, 0.05, 'unc', 'noverbose', 'nodisplay'); + +tc.verifyClass(out, 'glm_map'); +% Historical struct-style field access still works via Dependent aliases +tc.verifyClass(out.b, 'statistic_image'); % alias for betas +tc.verifyClass(out.t, 'statistic_image'); +tc.verifyClass(out.df, 'fmri_data'); +tc.verifyClass(out.sigma, 'fmri_data'); +tc.verifyTrue(isstruct(out.input_parameters)); +tc.verifyTrue(isstruct(out.diagnostics)); +tc.verifyNotEmpty(out.diagnostics.Variance_inflation_factors); +end + + +function test_construct_from_struct_and_aliases(tc) +% A struct with regress-style field names re-casts into a glm_map, and the +% out-struct aliases read the canonical properties. +S = struct(); +S.b = statistic_image; S.b.dat = randn(50, 2); +S.t = statistic_image; S.t.dat = randn(50, 2); +S.variable_names = {'slope', 'intercept'}; +S.C = [1 0]'; +S.contrast_names = {'slope'}; +S.analysis_name = 'from_struct_test'; + +g = glm_map(S); +tc.verifyClass(g, 'glm_map'); +tc.verifyEqual(g.analysis_name, 'from_struct_test'); +tc.verifyEqual(g.betas.dat, S.b.dat); % b -> betas +tc.verifyEqual(g.regressor_names, {'slope', 'intercept'}); % variable_names -> regressor_names +tc.verifyEqual(g.contrasts, [1 0]'); % C -> contrasts +tc.verifyEqual(g.contrast_names, {'slope'}); +% Round-trip: alias getters return the canonical values +tc.verifyEqual(g.b.dat, g.betas.dat); +tc.verifyEqual(g.variable_names, g.regressor_names); +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 +% Default canonical HRF: one column per condition + one intercept = 3 +tc.verifyEqual(size(g.X, 2), 3); +tc.verifyNotEmpty(g.regressor_names); +tc.verifyNumElements(g.onsets, 2); % onsets read-through to design + +% Timing: condition A's regressor should peak ~3 TRs after its first onset +% (canonical HRF peak ~6 s = 3 TRs); first onset 10 s -> TR 6 -> peak ~TR 9. +[~, pk] = max(g.X(:, 1)); +tc.verifyLessThanOrEqual(abs(pk - 9), 1); +end + + +function test_event_design_matches_onsets2fmridesign(tc) +% The high-resolution pipeline should match onsets2fmridesign across TRs, +% including fractional TRs. +ons = {[10 40 70 100 150]', [25 55 85 130 200]'}; +for TR = [2 1.3 2.5] + nscan = round(260 / TR); + d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', ons, 'condition_names', {'A', 'B'}); + w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok + g = glm_map(d); g = build_design(g); + Xgt = onsets2fmridesign(ons, TR, nscan * TR, 'hrf'); % intercept last col + tc.verifyEqual(size(g.X, 1), size(Xgt, 1)); + tc.verifyGreaterThan(corr(g.X(:, 1), Xgt(:, 1)), 0.999); + tc.verifyGreaterThan(corr(g.X(:, 2), Xgt(:, 2)), 0.999); +end +end + + +function test_replace_basis_set(tc) +TR = 2; nscan = 120; +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', {[10 40 70 100]' [25 55 85 115]'}, 'condition_names', {'A', 'B'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok + +g = glm_map(d); g.is_timeseries = true; g = build_design(g); +rng(0); sim = fmri_data; sim.dat = randn(50, nscan); +g = add_contrasts(g, [1 -1 0], {'AmB'}); +g = fit(g, sim, 'noverbose'); +tc.verifyEqual(g.num_regressors, 3); +tc.verifyTrue(g.is_fitted); + +[xBF_hires, ~] = fmri_spline_basis(TR, 'length', 20, 'nbasis', 4, 'order', 3); + +% Without data: rebuild, clear stale fit/contrasts, refresh diagnostics +g2 = replace_basis_set(g, 1, xBF_hires); +tc.verifyEqual(g2.num_regressors, 6); % 4 spline (A) + 1 hrf (B) + intercept +tc.verifyFalse(g2.is_fitted); % fitted maps cleared +tc.verifyEqual(g2.num_contrasts, 0); % contrasts cleared (columns changed) +tc.verifyNotEmpty(g2.diagnostics.Variance_inflation_factors); % diagnostics refreshed + +% With data: re-fit on the new design +g3 = replace_basis_set(g, 1, xBF_hires, 'data', sim, 'noverbose'); +tc.verifyTrue(g3.is_fitted); +tc.verifyEqual(size(g3.betas.dat, 2), 6); +end + + +% ===================================================================== +% Regressor-role indicators (of interest vs nuisance vs intercept) +% ===================================================================== +function test_interest_nuisance_indicators(tc) +% Event mode: conditions are of interest, baseline is the intercept +d = fmri_glm_design_matrix(2, 'nscan', 120, 'units', 'secs', ... + 'onsets', {[10 40 70]' [25 55 85]'}, 'condition_names', {'A' 'B'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = glm_map(d); g = build_design(g); +tc.verifyEqual(g.wh_interest, [true true false]); +tc.verifyEqual(g.wh_nuisance, [false false false]); +tc.verifyEqual(g.wh_intercept, [false false true]); + +% Direct mode: mark column 3 as a nuisance covariate; col 4 is the intercept +n = 30; +X = [zscore((1:n)') zscore(randn(n, 1)) zscore((n:-1:1)') ones(n, 1)]; +g2 = glm_map('X', X, 'level', 2, 'nuisance_columns', 3); +tc.verifyEqual(g2.wh_interest, [true true false false]); +tc.verifyEqual(g2.wh_nuisance, [false false true false]); +tc.verifyEqual(g2.wh_intercept, [false false false true]); +end + + +function test_diagnostics_interest_only(tc) +% A nuisance covariate correlated with a task regressor should inflate the +% full-design VIF relative to the interest-only VIF. +rng(1); n = 40; +task = zscore((1:n)'); +nuis = zscore(task + 0.4 * randn(n, 1)); +X = [task zscore(randn(n, 1)) nuis ones(n, 1)]; +g = glm_map('X', X, 'level', 2, 'nuisance_columns', 3); +g = run_diagnostics(g, 'noverbose'); + +tc.verifyNotEmpty(g.diagnostics.Variance_inflation_factors_interest_only); +% Full-design VIF for the task regressor exceeds its interest-only VIF +full_task = g.diagnostics.Variance_inflation_factors(1); +io_task = g.diagnostics.Variance_inflation_factors_interest_only(1); +tc.verifyGreaterThan(full_task, 2 * io_task); +tc.verifyNotEmpty(g.diagnostics.condition_number_interest_only); +end + + +% ===================================================================== +% import_onsets: tabular (FSL) and SPM-style cell arrays +% ===================================================================== +function test_import_onsets_variants(tc) +TR = 2; nscan = 120; + +% FSL/tabular with string condition names +T = table([5; 35; 65; 20; 50], [2; 2; 2; 3; 3], {'A'; 'A'; 'A'; 'B'; 'B'}, ... + 'VariableNames', {'onset', 'duration', 'name'}); +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs'); +d = import_onsets(d, T); +tc.verifyEqual(numel(d.Sess(1).U), 2); +tc.verifyEqual(d.Sess(1).U(1).ons(:)', [5 35 65]); + +% FSL/tabular with integer event-type codes +T2 = table([5; 35; 20; 50], [0; 0; 0; 0], [1; 1; 2; 2], ... + 'VariableNames', {'onset', 'duration', 'trial_type'}); +d2 = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs'); +d2 = import_onsets(d2, T2); +tc.verifyEqual(numel(d2.Sess(1).U), 2); + +% SPM-style cell arrays with durations + parametric modulators +d3 = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs'); +d3 = import_onsets(d3, {[10 40 70]' [25 55 85]'}, {4 4}, {[1 2 3]' []}, ... + 'names', {'cue' 'pain'}, 'pm_names', {'intensity' ''}); +tc.verifyEqual(d3.Sess(1).U(1).name, 'cue'); +tc.verifyEqual(d3.Sess(1).U(1).P.P(:)', [1 2 3]); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +d3 = build(d3); +tc.verifyEqual(size(d3.xX.X, 1), nscan); +end + + +% ===================================================================== +% Orthogonal contrast set + design efficiency +% ===================================================================== +function test_create_orthogonal_contrast_set(tc) +d = fmri_glm_design_matrix(2, 'nscan', 150, 'units', 'secs', ... + 'onsets', {[10 40 70]' [25 55 85]' [12 42 72]'}, 'condition_names', {'A' 'B' 'C'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = glm_map(d); g = build_design(g); + +g = create_orthogonal_contrast_set(g); +tc.verifyEqual(g.num_contrasts, 2); % 3 conditions -> 2 contrasts +% Orthogonal and zero on the intercept column +tc.verifyEqual(g.contrasts(end, :), [0 0]); % intercept row all zeros +offdiag = g.contrasts' * g.contrasts; offdiag(logical(eye(2))) = 0; +tc.verifyLessThan(max(abs(offdiag(:))), 1e-10); + +% Graceful error when there are no regressors of interest +g2 = glm_map('X', ones(20, 1), 'level', 2); +tc.verifyError(@() create_orthogonal_contrast_set(g2), 'glm_map:NoInterestRegressors'); +end + + +function test_diagnostics_efficiency(tc) +d = fmri_glm_design_matrix(2, 'nscan', 150, 'units', 'secs', ... + 'onsets', {[10 40 70]' [25 55 85]' [12 42 72]'}, 'condition_names', {'A' 'B' 'C'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = glm_map(d); g = build_design(g); + +% No contrasts entered: efficiency uses an auto orthogonal set over interest +g = run_diagnostics(g, 'noverbose'); +tc.verifyNotEmpty(g.diagnostics.efficiency); +tc.verifyEqual(numel(g.diagnostics.efficiency_per_contrast), 2); +tc.verifyTrue(contains(g.diagnostics.efficiency_contrast_source, 'orthogonal')); + +% No regressors of interest -> efficiency skipped gracefully (not an error) +g2 = glm_map('X', ones(20, 1), 'level', 2); +g2 = run_diagnostics(g2, 'noverbose'); +tc.verifyEmpty(g2.diagnostics.efficiency); +end + + +function test_diagnostics_scaled_condition_number(tc) +% A well-separated 2-condition event design with an intercept is nearly +% orthogonal. The raw cond(X) is huge because the tiny event regressors and +% the unit intercept differ in scale; the scaled condition number should be +% small and consistent with the low VIFs. +TR = 2; nscan = 100; +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', {[10 40 70 100]' [25 55 85 115]'}, 'condition_names', {'A' 'B'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = glm_map(d); g.is_timeseries = true; g = build_design(g); +g = run_diagnostics(g, 'noverbose'); + +tc.verifyGreaterThan(cond(g.X), 100); % raw cond is large (scaling) +tc.verifyLessThan(g.diagnostics.condition_number, 10); % scaled cond is small +tc.verifyLessThan(max(g.diagnostics.Variance_inflation_factors), 4); +% g.X is exposed at the top level and links to the design matrix +tc.verifyEqual(g.X, g.design.xX.X); +end + + +function test_diagnostics_hpfilter_timeseries(tc) +% Cumulative-power / high-pass-filter assessment is computed only for +% timeseries designs, and separately for regressors and contrasts. +TR = 2; nscan = 100; +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', {[10 40 70 100]' [25 55 85 115]'}, 'condition_names', {'A' 'B'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = glm_map(d); g.is_timeseries = true; g = build_design(g); +g = add_contrasts(g, [1 -1 0], {'A_vs_B'}); +g = run_diagnostics(g, 'noverbose'); + +hp = g.diagnostics.hpfilter; +tc.verifyTrue(isstruct(hp)); +tc.verifyTrue(isfield(hp, 'recommended_hpfilter_sec_regressors')); +tc.verifyTrue(isfield(hp, 'recommended_hpfilter_sec_contrasts')); +tc.verifyEqual(size(hp.cumulative_power_regressors, 2), 2); % one column per event regressor +tc.verifyNotEmpty(hp.recommended_hpfilter_sec_contrasts); % contrast was entered + +% Non-timeseries (group) design: no HP-filter assessment +g2 = glm_map('X', [zscore((1:30)') ones(30, 1)], 'level', 2); +g2 = run_diagnostics(g2, 'noverbose'); +tc.verifyEmpty(g2.diagnostics.hpfilter); +end + + +function test_glm_map_import_onsets_and_spm_flags(tc) +% glm_map.import_onsets bootstraps a design, builds it, and flags events +T = table([5; 35; 65; 20; 50], [0; 0; 0; 0; 0], {'A'; 'A'; 'A'; 'B'; 'B'}, ... + 'VariableNames', {'onset', 'duration', 'name'}); +w = warning('off', 'all'); c = onCleanup(@() warning(w)); %#ok +g = import_onsets(glm_map, T, 'TR', 2, 'nscan', 60, 'units', 'secs'); +tc.verifyEqual(size(g.X, 1), 60); +tc.verifyEqual(g.wh_interest, [true true false]); % A,B of interest; intercept not + +% Missing design + no TR/nscan -> graceful error +tc.verifyError(@() import_onsets(glm_map, T), 'glm_map:NoDesignForImport'); + +% SPM import: event semantics drive the interest flags even when SPM.xX.iH is empty +SPM = local_make_synthetic_spm(); +g2 = import_SPM(glm_map, SPM, 'noverbose'); +% names: 'Sn(1) Cue*bf(1)', 'Sn(1) Pain*bf(1)', 'Sn(1) constant' +tc.verifyEqual(g2.wh_interest, [true true false]); +tc.verifyEqual(g2.wh_intercept, [false false true]); +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 diff --git a/docs/workflows/regression_with_glm_map.m b/docs/workflows/regression_with_glm_map.m new file mode 100644 index 00000000..f2fbf496 --- /dev/null +++ b/docs/workflows/regression_with_glm_map.m @@ -0,0 +1,281 @@ +%% Mass-univariate regression with the glm_map object +% +% This walkthrough shows two equivalent ways to run a voxelwise multiple +% regression in CanlabCore and work with the result as a |glm_map| object: +% +% (1) the quick path - call fmri_data.regress directly; it now *returns* +% a glm_map, which you then access, display, and query. +% (2) the estimator path - build a glm_map, attach data, and call fit(), +% add_contrasts(), run_diagnostics(), threshold(), etc. +% (a scikit-learn-style API). +% +% It mirrors the worked examples in the fmri_data.regress help (naming, +% nuisance covariates, gray-matter masking, outlier handling, robust fits, +% residuals, brain-predicts-behavior, writing maps, diagnostics) and shows +% how the object groups related outputs into nested structs +% (.input_parameters, .input_image_metadata, .diagnostics). +% +% Requirements: CanlabCore + Neuroimaging_Pattern_Masks + SPM12 on the path. +% Every block below is copy-pasteable. +% +% See also: fmri_data.regress, glm_map, statistic_image, fmri_glm_design_matrix + + +%% 1. Load sample data and define a real predictor from the metadata table +% load_image_set('emotionreg') returns an fmri_data object with 30 first-level +% contrast images (one per participant; Wager et al. 2008, Neuron) plus a +% metadata_table. We predict each participant's brain map from their +% behavioral reappraisal success score (a 2nd-level / group regression). + +obj = load_image_set('emotionreg', 'noverbose'); % 30 images +n = size(obj.dat, 2); % number of images + +% Predictor matrix from metadata_table. A single column here; the intercept is +% added automatically if absent. Mean-centering the regressor(s) makes the +% intercept map interpretable as the activation for the average participant. +obj.X = obj.metadata_table.Reappraisal_Success; +obj.X = obj.X - mean(obj.X); % mean-center +obj.X(:, end + 1) = 1; % regressor + intercept + +X = obj.X; % keep for later sections + + +%% 2. Quick path: fmri_data.regress returns a glm_map +% Run the regression at a liberal threshold. The result "g" is a glm_map +% object. (Historically regress returned a plain struct; it now re-casts that +% struct as a glm_map. All the old struct-style field names still work - see +% section 6.) + +g = regress(obj, .05, 'unc', 'noverbose', 'nodisplay'); + +g % typing the name lists ALL properties + +% The thresholded t map (one statistic_image per regressor, incl. intercept): +g.t + + +%% 3. Name the regressors and the analysis +% Naming labels the output maps (.image_labels), which propagates to montages +% and tables; analysis_name is carried on the object and printed by summary. + +g = regress(obj, .05, 'unc', 'names', {'Reapp_Success' 'Intercept'}, ... + 'analysis_name', 'Emotion Regulation', 'noverbose', 'nodisplay'); + +g.t.image_labels % {'Reapp_Success','Intercept'} + + +%% 4. Display the result maps +% There is one image per regressor in X. Select, visualize, and re-threshold +% them with the statistic_image display methods. + +montage(g.t); % all regressor t maps + +t_reapp = get_wh_image(g.t, 1); % just the Reapp_Success map +create_figure('surface'); surface(t_reapp); + +t_reapp = threshold(t_reapp, .005, 'unc'); % re-threshold a single map +create_figure('montage P<.005'); axis off; montage(t_reapp); + +orthviews(g.t); % orthviews of the thresholded maps + +% glm_map also has object-level display wrappers that pick a map by name: +montage(g, 't'); +% table(g, 't'); % atlas-labeled results table + + +%% 5. Query the object and plot design diagnostics +% Result maps are statistic_image / fmri_data objects inside the glm_map. + +g.betas % statistic_image: one beta image per regressor +g.df % fmri_data: per-voxel error degrees of freedom +g.sigma % fmri_data: per-voxel residual standard deviation +g.dfe % scalar error df (median of g.df) + +% Design diagnostics are collected in the nested .diagnostics struct, using +% the same field names as fmri_data.regress out.diagnostics. +g = run_diagnostics(g); % Run diagnostics and return output in object +g.diagnostics.Variance_inflation_factors % VIF per regressor +g.diagnostics.Contrast_variance_inflation_factors % contrast VIFs (cVIF), if contrasts present +g.diagnostics.Leverages % per-observation leverage +g.diagnostics.condition_number % conditioning of X +g.diagnostics.Cooks_distance % per-observation influence + +% Options the fit used live in .input_parameters; input-image provenance in +% .input_image_metadata. +g.input_parameters +g.input_image_metadata + +% Plot the design diagnostics +figure; +subplot(1, 2, 1); plot(g.diagnostics.Variance_inflation_factors, 'o-'); title('VIFs'); xlabel('Regressor'); +subplot(1, 2, 2); plot(g.diagnostics.Leverages, 'o-'); title('Leverage'); xlabel('Observation'); + +% use the plot_design method to do it: +plot_design(g); + + +%% 6. The object mirrors the regress out-struct (with back-compatible aliases) +% glm_map property names match the field names of "out" in fmri_data.regress. +% Where the historical field name differs from the canonical property, an alias +% reads/writes the same data, so legacy struct-style access is unchanged: +% +% out.b -> g.betas +% out.contrast_images -> g.contrast_estimates +% out.con_t -> g.contrast_t +% out.resid -> g.residuals +% out.variable_names -> g.regressor_names +% out.C -> g.contrasts + +isequal(g.b.dat, g.betas.dat) % true: .b is an alias for .betas +isequal(g.variable_names, g.regressor_names) % true + +% You can also re-cast any regress-style struct into a glm_map yourself: +% g2 = glm_map(out_struct); + + +%% 7. Add a 2nd-level nuisance covariate (mean CSF) +% Covariates for effects of no interest stabilize results. Mean gray, white, +% and CSF signals are highly correlated; CSF likely carries image-wide +% confounds rather than task signal, so it is a reasonable nuisance covariate. + +gwcsf = extract_gray_white_csf(obj); % [n x 3]: gray, white, CSF means + +obj.X = [obj.metadata_table.Reappraisal_Success, gwcsf(:, 3)]; +obj.X = obj.X - mean(obj.X); % mean-center predictors +obj.X(:, end + 1) = 1; % intercept + +g = regress(obj, .05, 'unc', 'names', {'Reapp_Success' 'CSF_mean' 'Intercept'}, ... + 'analysis_name', 'Emotion Regulation', 'noverbose', 'nodisplay'); +montage(g.t); + + +%% 8. Apply a gray-matter mask before thresholding (FDR) +% Correcting within hypothesized gray-matter voxels reduces the comparison +% count and can make results more significant. + +gm = fmri_data(which('gray_matter_mask.nii')); +g.t = apply_mask(g.t, gm); +g.t = threshold(g.t, 0.05, 'fdr'); +montage(g.t); + + +%% 9. Find and exclude outliers, then re-run +% 'notimeseries' omits time-series-specific calculations (this is 2nd-level +% data). We flag images that do not correlate with the others (mahal_corr). + +[~, wh_outliers] = outliers(obj, 'notimeseries'); +find(wh_outliers) % indices flagged as outliers + +obj_clean = get_wh_image(obj, ~wh_outliers); % drop flagged images +obj_clean.X = obj.X(~wh_outliers, :); % and the matching design rows + +g_clean = regress(obj_clean, .05, 'unc', ... + 'names', {'Reapp_Success' 'CSF_mean' 'Intercept'}, 'noverbose', 'nodisplay'); +g_clean.t = threshold(apply_mask(g_clean.t, gm), 0.01, 'unc'); +montage(g_clean.t); + + +%% 10. More regress options +% Reset to the single-predictor design for these. +obj.X = X; + +% Robust regression (iteratively down-weights outliers; slower than OLS): +g_robust = regress(obj, .05, 'fdr', 'robust', 'names', {'Reapp_Success' 'Intercept'}, ... + 'noverbose', 'nodisplay'); + +% Save voxelwise residuals as an fmri_data object (useful for denoising / QC): +g_resid = regress(obj, .001, 'unc', 'residual', 'noverbose', 'nodisplay'); +g_resid.residuals % fmri_data [voxels x images] + +% Brain-predicts-behavior ('brainony'): a univariate map of where brain +% activity predicts obj.Y. Needs obj.Y and is slow (loops over voxels): +% obj.Y = obj.metadata_table.Reappraisal_Success; +% g_bxy = regress(obj, .05, 'unc', 'brainony'); + +% Re-threshold and re-display an existing map without refitting: +g.t = threshold(g.t, .001, 'unc'); +orthviews(g.t); + +% Write a beta image to disk (uncomment to run; writes to the current folder): +% g.betas.fullpath = fullfile(pwd, 'beta_reapp_success.nii'); +% write(g.betas); + + +%% 11. Estimator path: build, screen, fit +% The same analysis as a reusable estimator. Construct the design first (no +% data), screen its collinearity, add contrasts, then fit. fit() calls +% fmri_data.regress under the hood and stores the result maps on the object. + +g = glm_map('X', X, 'level', 2, 'regressor_names', {'Reapp_Success' 'Intercept'}, ... + 'analysis_name', 'Emotion Regulation (estimator)'); + +g = add_contrasts(g, [1 0], {'reapp_effect'}); % one row per contrast + +g = run_diagnostics(g); % VIF / cVIF / leverage / conditioning (no fit) + +g = fit(g, obj); % runs the regression + +g.is_fitted % true +summary(g); % narrative report: model, diagnostics, results +montage(g, 'contrast_t'); % thresholded contrast t map +% table(g, 'contrast'); % table for the contrast map +% +% summary(g) prints the analysis name, the model (level + input variables), +% design diagnostics (max VIF, contrast VIF, condition number, leverage, max +% Cook's distance), and -- once fitted -- how many maps, the statistical +% threshold, and the number of significant voxels per regressor and contrast. + + +%% 12. Re-threshold without refitting +% threshold() re-masks the stored t / contrast_t maps; the underlying +% statistic values are preserved, so you can sweep thresholds cheaply. + +g = threshold(g, .005, 'unc', 'k', 10); % both t and contrast_t +g = threshold(g, .05, 'fdr', 'which_map', 'contrast'); % contrast map only + + +%% 13. First-level (event) designs and AR models +% For within-run BOLD timeseries, wrap an fmri_glm_design_matrix in .design; +% glm_map builds X by HRF-convolving the onsets. Marking is_timeseries lets +% fit() use autoregressive (AR) error models. build_design runs here; the fit +% and import_SPM lines need real timeseries / an SPM.mat, so they are shown +% for reference. + +TR = 2; nscan = 100; % nscan in TRs +onsets = {[10 40 70 100]', [25 55 85 115]'}; % seconds, two conditions +d = fmri_glm_design_matrix(TR, 'nscan', nscan, 'units', 'secs', ... + 'onsets', onsets, 'condition_names', {'Event A' 'Event B'}); +g_evt = glm_map(d); % level 1, event mode +g_evt.is_timeseries = true; +g_evt = build_design(g_evt); % onsets -> X via convolution + +% Add a simple contrast comparing the two events: A vs B. The design has +% three columns (Event A, Event B, intercept), so the contrast is [1 -1 0] +% (add_contrasts takes one row per contrast; the intercept gets 0). +g_evt = add_contrasts(g_evt, [1 -1 0], {'A_vs_B'}); + +% Plot the design BEFORE diagnostics: the design matrix heat map and one panel +% per event type (the actual basis-convolved regressors). No VIF figure yet. +plot_design(g_evt); + +% g_evt = fit(g_evt, bold_timeseries_fmri_data, 'AR', 1); % AR(1) error model +% g_evt = import_SPM(glm_map, '/path/to/SPM.mat'); % import a full 1st-level model + +% Now run diagnostics and plot again. run_diagnostics computes VIFs/cVIFs and +% (because is_timeseries is true) the cumulative power by frequency and a +% recommended high-pass filter cutoff. plot_design now ALSO draws a VIF figure +% (VIFs of regressors of interest, with/without nuisance, plus contrast VIFs, +% with severity reference lines at 1/2/4/8). +g_evt = run_diagnostics(g_evt); % prints the diagnostics report +plot_design(g_evt); % design + VIF/cVIF figure +g_evt.diagnostics.hpfilter % recommended HP filter (timeseries) + +%% 14. Summary +% - fmri_data.regress returns a glm_map; access/display it via .betas/.t/ +% montage/surface/orthviews, query the nested .diagnostics / +% .input_parameters / .input_image_metadata structs, or use the historical +% out-struct aliases (.b/.con_t/.contrast_images/...). +% - Or build a glm_map as an estimator: glm_map(...) -> add_contrasts -> +% diagnostics -> fit -> summary/threshold/table/montage. +% - typing the object name lists all properties; summary(g) prints a narrative +% report; methods(glm_map) lists all operations.