diff --git a/.github/workflows/tests-walkthroughs.yml b/.github/workflows/tests-walkthroughs.yml index e1617ef2..f7664896 100644 --- a/.github/workflows/tests-walkthroughs.yml +++ b/.github/workflows/tests-walkthroughs.yml @@ -49,6 +49,7 @@ jobs: products: | Statistics_and_Machine_Learning_Toolbox Signal_Processing_Toolbox + Image_Processing_Toolbox cache: true - name: Run walkthrough tests diff --git a/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m new file mode 100644 index 00000000..cb750be2 --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_classify_environment_error.m @@ -0,0 +1,168 @@ +function category = canlab_classify_environment_error(ME, had_env_skip) +%CANLAB_CLASSIFY_ENVIRONMENT_ERROR Bucket a caught error for CI test handling. +% +% category = canlab_classify_environment_error(ME) +% category = canlab_classify_environment_error(ME, had_env_skip) +% +% Classifies an MException caught while running toolbox example/walkthrough +% code into one of a small set of buckets, so that test harnesses can decide +% whether to SKIP (mark Incomplete) or FAIL. The goal is to keep headless CI +% from going red over missing-display / missing-data / interactive-prompt +% conditions that are not real regressions, while still surfacing genuine +% code errors. +% +% This centralizes the heuristics previously hand-coded inside +% canlab_test_help_examples.m (skip_on_environment_error) so the walkthrough +% harness and the per-method help-example tests share one definition. +% +% :Inputs: +% +% **ME:** +% an MException (or any struct with .identifier, .message, and +% optionally .stack fields). +% +% **had_env_skip:** +% logical (default false). When true, "undefined variable/function" +% and "unrecognized field" errors are treated as CASCADE fallout from +% an earlier skipped section rather than genuine failures. Pass the +% running "have we already skipped something upstream" flag here. +% +% :Output: +% +% **category:** +% one of: +% 'graphics' - needs a display / OpenGL / Java / figure window, or +% the error originates inside a graphics-only code path +% (orthviews, surface rendering, etc.). +% 'input' - code tried to prompt for interactive user input, +% unavailable in batch/CI. +% 'capability' - an UndefinedFunction error for a known optional +% MATLAB toolbox function (e.g. niftiinfo) that is not +% provisioned on this runner. +% 'data' - an optional data file (signature, atlas, feature set) +% is not on the CI path. +% 'cascade' - an undefined-variable/field error following an earlier +% environment skip (downstream fallout, not a new bug). +% 'genuine' - none of the above; treat as a real failure. +% +% :See also: canlab_run_walkthrough_snapshot, canlab_test_help_examples + +% .. +% Author: CANlab. Part of the CanlabCore Unit_tests headless-CI harness. +% .. + +if nargin < 2 || isempty(had_env_skip) + had_env_skip = false; +end + +msg = lower(ME.message); +id = ''; +if isprop(ME, 'identifier') || isfield(ME, 'identifier') + id = ME.identifier; +end + +% --------------------------------------------------------------------- +% Graphics: explicit graphics error ids, telltale message tokens, or a +% stack frame inside a known graphics-only function. spm_orthviews and the +% surface renderers fail or behave erratically on a headless runner. +% --------------------------------------------------------------------- +gfx_ids = {'MATLAB:graphics:opengl:Unavailable', ... + 'MATLAB:graphics:initialize', ... + 'MATLAB:class:InvalidHandle'}; +is_gfx = any(strcmp(id, gfx_ids)) || ... + strncmp(id, 'MATLAB:hg:', 10) || ... % handle-graphics errors + contains(msg, 'opengl') || contains(msg, 'display') || ... + contains(msg, 'java') || contains(msg, 'jvm') || ... + contains(msg, 'figure window') || contains(msg, 'no display') || ... + contains(msg, 'graphics') || ... + contains(msg, 'invalid or deleted object'); % set/get on a handle + % whose graphics section + % was skipped upstream + +gfx_fns = {'spm_orthviews', 'orthviews', 'spm_check_registration', ... + 'spm_figure', 'surface', 'render_on_surface', 'addbrain', ... + 'cluster_surf', 'isosurface', 'riverplot', 'tor_3d', ... + 'render_blobs', 'cluster_orthviews'}; +if ~is_gfx && isfield_or_prop(ME, 'stack') && ~isempty(ME.stack) + stack_names = {ME.stack.name}; + if any(ismember(stack_names, gfx_fns)) + is_gfx = true; + end +end +if is_gfx + category = 'graphics'; + return +end + +% --------------------------------------------------------------------- +% Interactive input: input() / keyboard prompts are unavailable in -batch. +% MATLAB reports MissingRequiredCapability with a "support for user input" +% message in that case. +% --------------------------------------------------------------------- +is_input = strcmp(id, 'MATLAB:services:MissingRequiredCapability') || ... + contains(msg, 'support for user input') || ... + contains(msg, 'input is not available'); +if is_input + category = 'input'; + return +end + +% --------------------------------------------------------------------- +% Missing optional MATLAB toolbox. An UndefinedFunction error for a known +% MathWorks toolbox function (e.g. niftiinfo/niftiread from the Image +% Processing Toolbox) means that toolbox is not provisioned on this runner, +% not that the code is broken. Keep this list to functions that are *only* +% ever provided by a toolbox, so we never mask a genuine missing-CanlabCore +% function bug. +% --------------------------------------------------------------------- +optional_toolbox_fns = {'niftiinfo', 'niftiread', 'niftiwrite', ... + 'cfg_getfile'}; +if strcmp(id, 'MATLAB:UndefinedFunction') + undef = regexp(ME.message, "Undefined function '([^']+)'", 'tokens', 'once'); + if ~isempty(undef) && any(strcmpi(undef{1}, optional_toolbox_fns)) + category = 'capability'; + return + end +end + +% --------------------------------------------------------------------- +% Missing optional data files (NPS+ signatures, Neurosynth feature set, +% Bianciardi atlas sources, etc.). load_image_set / annotate_* abort with a +% disp() + bare error('Exiting'), so the message is literally "Exiting". +% --------------------------------------------------------------------- +bare_exiting = isempty(id) && ... + ismember(strtrim(strrep(ME.message, '.', '')), {'Exiting', 'exiting'}); +is_data = contains(msg, 'cannot find images') || ... + contains(msg, 'find and add the file') || ... + contains(msg, 'not found in matlab path') || ... + contains(msg, 'no such file or directory') || ... + contains(msg, 'cannot find') || ... + bare_exiting; +if is_data + category = 'data'; + return +end + +% --------------------------------------------------------------------- +% Cascade: an earlier section was skipped, so a variable/field it would +% have defined is now missing. Treat as fallout, not a new regression. +% --------------------------------------------------------------------- +if had_env_skip && ( ... + strcmp(id, 'MATLAB:UndefinedFunction') || ... + strcmp(id, 'MATLAB:undefinedVarOrClass') || ... + strcmp(id, 'MATLAB:nonExistentField') || ... + contains(msg, 'undefined') || ... + contains(msg, 'unrecognized field')) + category = 'cascade'; + return +end + +category = 'genuine'; + +end + + +function tf = isfield_or_prop(obj, name) +% MException is an object (use isprop); a plain struct uses isfield. +tf = (isobject(obj) && isprop(obj, name)) || (isstruct(obj) && isfield(obj, name)); +end diff --git a/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m b/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m new file mode 100644 index 00000000..e8dcce50 --- /dev/null +++ b/CanlabCore/Unit_tests/helpers/canlab_run_walkthrough_snapshot.m @@ -0,0 +1,207 @@ +function canlab_run_walkthrough_snapshot(tc, wt__snapshot_path) +%CANLAB_RUN_WALKTHROUGH_SNAPSHOT Run a frozen walkthrough snapshot under CI hardening. +% +% canlab_run_walkthrough_snapshot(tc, snapshot_path) +% +% Executes a copied CANlab_help_examples walkthrough script (a "snapshot" +% bundled under Unit_tests/walkthroughs/private/) one %%-cell at a time, in +% a single shared workspace, with headless-graphics and fault-tolerant error +% handling. Reports outcomes to the supplied matlab.unittest TestCase `tc`. +% +% Why snapshots instead of running the real tutorials. The walkthroughs in +% the CANlab_help_examples repo are teaching scripts, full of graphics +% (orthviews, surface), optional data sets, and the occasional interactive +% prompt. They should stay clean tutorials, not be retrofitted with test +% scaffolding, and CanlabCore CI should not depend on an external repo's +% un-hardened scripts. So we keep verbatim copies here and apply the +% hardening at run time, in one reusable place (this function) rather than +% editing each copy. Refresh a snapshot by overwriting the file under +% private/ from example_help_files/ — no re-hardening needed. +% +% Hardening applied: +% * Headless graphics: DefaultFigureVisible is forced 'off' for the run +% (callers also set this; doing it here makes the harness safe to call +% directly). Figures are closed between sections. +% * Section isolation: the script is split at its %% cell boundaries and +% each cell is run in its own try/catch, so a graphics-only section that +% fails on a headless runner does not abort the compute sections. +% * Error classification (see canlab_classify_environment_error): +% - graphics / input / data / cascade errors -> section SKIPPED, +% recorded, execution continues. +% - genuine errors -> the harness stops and the test FAILS with an +% informative message (which section, error id, message, and the +% offending top stack frame), because downstream cells usually +% cascade once a real error occurs. +% * Outcome mapping to matlab.unittest: +% - any genuine failure -> tc.verifyFail (test FAILS) +% - every section skipped (env) -> tc.assumeFail (test INCOMPLETE) +% - at least one section ran ok -> PASS, with a logged skip summary. +% +% Variable hygiene: every internal variable in this function is prefixed +% `wt__` because snapshot code is eval'd in this workspace and would +% otherwise clobber ordinary loop variables (i, n, t, r, ...). +% +% :Inputs: +% **tc:** a matlab.unittest.TestCase (from functiontests). +% **snapshot_path:** absolute path to the snapshot .m file to run. +% +% :See also: canlab_classify_environment_error, canlab_run_all_tests + +% .. +% Author: CANlab. Part of the CanlabCore Unit_tests headless-CI harness. +% .. + +tc.assertTrue(exist(wt__snapshot_path, 'file') == 2, ... + sprintf('Walkthrough snapshot not found: %s', wt__snapshot_path)); + +[~, wt__name] = fileparts(wt__snapshot_path); +wt__src = fileread(wt__snapshot_path); +wt__sections = local_split_sections(wt__src); +wt__nsec = numel(wt__sections); + +% Headless: callers set this too, but be self-sufficient. +wt__prev_vis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); +wt__cleanup = onCleanup(@() set(0, 'DefaultFigureVisible', wt__prev_vis)); + +wt__skipped = {}; % cellstr of human-readable skip notes +wt__ran_real = 0; % count of sections that executed without error +wt__had_env_skip = false; % drives cascade classification downstream +wt__failure = []; % first genuine failure, if any (struct) + +for wt__s = 1:wt__nsec + + wt__code = wt__sections{wt__s}; + + % Skip cells that are pure comment/whitespace, or that define a + % function (snapshots are scripts; a trailing function would not eval). + if local_is_noop(wt__code) + continue + end + + try + evalc(wt__code); % runs in THIS workspace; vars persist + wt__ran_real = wt__ran_real + 1; + catch wt__ME + wt__cat = canlab_classify_environment_error(wt__ME, wt__had_env_skip); + if strcmp(wt__cat, 'genuine') + wt__failure = struct( ... + 'section', wt__s, ... + 'id', wt__ME.identifier, ... + 'message', wt__ME.message, ... + 'where', local_top_frame(wt__ME), ... + 'snippet', local_first_stmt(wt__code)); + break % stop; later cells will cascade + else + wt__had_env_skip = true; + wt__skipped{end+1, 1} = sprintf(' section %d/%d [%s]: %s', ... + wt__s, wt__nsec, wt__cat, local_oneline(wt__ME.message)); %#ok + end + end + + close all force +end + +% --------------------------------------------------------------------- +% Map collected outcomes onto the TestCase. +% --------------------------------------------------------------------- +if ~isempty(wt__failure) + wt__report = sprintf([ ... + '%s: genuine error at section %d of %d.\n' ... + ' identifier: %s\n' ... + ' message: %s\n' ... + ' section starts: %s\n' ... + ' at: %s\n' ... + ' (%d earlier section(s) skipped for environment reasons)'], ... + wt__name, wt__failure.section, wt__nsec, ... + wt__failure.id, local_oneline(wt__failure.message), ... + wt__failure.snippet, wt__failure.where, numel(wt__skipped)); + tc.verifyFail(wt__report); + +elseif wt__ran_real == 0 + tc.assumeFail(sprintf( ... + '%s: all %d section(s) skipped for environment reasons:\n%s', ... + wt__name, wt__nsec, strjoin(wt__skipped, newline))); + +else + if ~isempty(wt__skipped) + tc.log(1, sprintf('%s: %d of %d section(s) skipped (environment):\n%s', ... + wt__name, numel(wt__skipped), wt__nsec, strjoin(wt__skipped, newline))); + end + tc.verifyTrue(true, sprintf('%s ran %d section(s) without genuine error', ... + wt__name, wt__ran_real)); +end + +end + + +% ===================================================================== +% Local helpers +% ===================================================================== + +function sections = local_split_sections(src) +% Split source at %% cell markers. Each section keeps its leading marker +% line (a comment, harmless to eval). Everything before the first marker is +% section 1. +lines = regexp(src, '\r\n|\r|\n', 'split'); +is_head = ~cellfun('isempty', regexp(lines, '^\s*%%', 'once')); +starts = find(is_head); +if isempty(starts) || starts(1) ~= 1 + starts = [1, starts]; +end +starts = unique(starts); +ends = [starts(2:end) - 1, numel(lines)]; +sections = cell(numel(starts), 1); +for k = 1:numel(starts) + sections{k} = strjoin(lines(starts(k):ends(k)), newline); +end +end + + +function tf = local_is_noop(code) +% True if the cell has no executable content: blank, comment-only, or a +% function definition (which cannot be eval'd as a statement). +stripped = regexprep(code, '%[^\n]*', ''); % drop line comments +stripped = strtrim(stripped); +tf = isempty(stripped) || ~isempty(regexp(stripped, '^function\b', 'once')); +end + + +function s = local_oneline(msg) +% Collapse a multi-line error message to a single trimmed line. +s = strtrim(regexprep(msg, '\s+', ' ')); +end + + +function s = local_top_frame(ME) +% "funcname (line N)" for the deepest stack frame that is not this harness +% (errors raised at a section's top level otherwise just point back at the +% evalc call site here, which is unhelpful). Falls back to a marker. +if isempty(ME.stack) + s = ''; + return +end +this_fn = mfilename; +for k = 1:numel(ME.stack) + if ~strcmp(ME.stack(k).name, this_fn) + s = sprintf('%s (line %d)', ME.stack(k).name, ME.stack(k).line); + return + end +end +s = ''; +end + + +function s = local_first_stmt(code) +% First non-comment, non-blank line of a section, trimmed, for the report. +lines = regexp(code, '\r\n|\r|\n', 'split'); +s = ''; +for k = 1:numel(lines) + ln = strtrim(lines{k}); + if ~isempty(ln) && ~strncmp(ln, '%', 1) + s = ln; + if numel(s) > 100, s = [s(1:97) '...']; end + return + end +end +end diff --git a/CanlabCore/Unit_tests/walkthroughs/README.md b/CanlabCore/Unit_tests/walkthroughs/README.md new file mode 100644 index 00000000..4bfbbd52 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/README.md @@ -0,0 +1,85 @@ +# Walkthrough integration tests (Tier B) + +These tests smoke-test the end-to-end CANlab tutorials ("walkthroughs") from +the [CANlab_help_examples](https://github.com/canlab/CANlab_help_examples) +repository, to catch the case where a tutorial silently bit-rots against the +current toolbox. + +They are **slow** and depend on more sibling repos and optional data than the +fast per-push unit suite, so they are a separate tier, run on the nightly +`tests-walkthroughs` GitHub Actions workflow, **not** on every push. They are +skipped by `canlab_run_all_tests` unless you ask for them: + +```matlab +canlab_run_all_tests('Walkthroughs', 'only') % just the walkthroughs +canlab_run_all_tests('Walkthroughs', 'include') % unit suite + walkthroughs +canlab_run_all_tests % default: walkthroughs skipped +``` + +## Design: snapshots + a hardening harness + +The tutorials in CANlab_help_examples are teaching scripts. They should stay +clean tutorials — full of `orthviews`, `surface`, optional signature/atlas +data, and the occasional interactive prompt — and **not** be retrofitted with +test scaffolding. CanlabCore CI should also not depend on an external repo's +un-hardened scripts at run time. So instead of executing the tutorials in +place, this folder works from frozen copies: + +``` +walkthroughs/ + canlab_test_walkthrough_*.m <- thin matlab.unittest wrappers (discovered) + private/ + canlab_help_*.m <- VERBATIM snapshots of the tutorials + README.md +``` + +- `private/` is excluded from `genpath`, so the snapshots never go on the + MATLAB path and never collide with the real tutorials when both repos are + checked out on CI. The harness reads them by absolute path. +- Each `canlab_test_walkthrough_*.m` wrapper just points at its snapshot and + calls the shared harness. The hardening lives in one place + (`Unit_tests/helpers/canlab_run_walkthrough_snapshot.m`), not duplicated into + each copy — so refreshing a snapshot is a clean file overwrite. + +### What the harness does + +`canlab_run_walkthrough_snapshot(tc, snapshot_path)`: + +1. Forces `DefaultFigureVisible 'off'` (headless). +2. Splits the snapshot at its `%%` cell boundaries and runs each cell in its + own `try/catch`, in one shared workspace, so a graphics-only section that + fails on a headless runner does not abort the compute sections. +3. Classifies any caught error via + `Unit_tests/helpers/canlab_classify_environment_error.m`: + - **graphics / input / data / cascade** → section **skipped** (recorded), + execution continues. + - **genuine** → the harness stops and the test **fails** with an informative + report (which section, the offending line, error id and message). +4. Maps outcomes to `matlab.unittest`: + - any genuine failure → **Failed** + - every section skipped for environment reasons → **Incomplete** + - at least one section ran → **Passed** (with a logged skip summary) + +So the nightly goes red only on a *genuine* regression, and tells you exactly +where; missing-display / missing-data / prompt conditions become skips. + +## Refreshing a snapshot + +When a tutorial changes upstream, re-copy it verbatim: + +```bash +cp ../../../../CANlab_help_examples/example_help_files/canlab_help_5_regression_walkthrough.m \ + private/canlab_help_5_regression_walkthrough.m +``` + +No re-hardening needed — the harness handles it. Keep the copies verbatim so +they can be diffed against upstream. + +## Why a separate nightly tier (not folded into the per-push suite) + +Measured wall-time for the full set is ~5 min on a fast workstation and an +estimated ~15–20 min on the 2-core GitHub Actions Linux runner. The per-push +`tests` suite is the fast gate that every PR blocks on; adding 15–20 min plus +graphics/data-dependent flakiness to it would slow iteration and make the +required check unreliable. End-to-end tutorial runs are exactly the kind of +slow, environment-sensitive integration test that belongs in a nightly tier. diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m index 5398e0e4..01db58ce 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_1_installing_tools.m @@ -1,34 +1,45 @@ function tests = canlab_test_walkthrough_1_installing_tools -%CANLAB_TEST_WALKTHROUGH_1_INSTALLING_TOOLS End-to-end smoke test of canlab_help_1_installing_tools. +%CANLAB_TEST_WALKTHROUGH_1_INSTALLING_TOOLS Hardened smoke test: Installing tools walkthrough. % -% Wrapper around CANlab_help_examples/example_help_files/canlab_help_1_installing_tools.m. -% Acts as a Tier B integration test: confirms the walkthrough runs end-to-end -% on the current toolbox without erroring. Does not check pixel content, -% statistics, or output values; reaching the end of the script counts as a pass. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_1_installing_tools" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. % -% Skipped (Incomplete) if CANlab_help_examples is not on the MATLAB path. +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> -% Run each walkthrough in its own temp working directory so output files -% (NIfTI writes, saved figures) don't collide with prior runs. +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_1_installing_tools.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_1_installing_tools'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m index 0ea0f7cd..ba5d6a47 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2_load_a_sample_dataset.m @@ -1,29 +1,45 @@ function tests = canlab_test_walkthrough_2_load_a_sample_dataset -%CANLAB_TEST_WALKTHROUGH_2_LOAD_A_SAMPLE_DATASET Smoke test of canlab_help_2_load_a_sample_dataset. +%CANLAB_TEST_WALKTHROUGH_2_LOAD_A_SAMPLE_DATASET Hardened smoke test: Load a sample dataset walkthrough. % -% Tier B integration test: runs the walkthrough end-to-end and counts a -% return without error as a pass. See canlab_test_walkthrough_1_installing_tools -% for shared rationale. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2_load_a_sample_dataset" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2_load_a_sample_dataset.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2_load_a_sample_dataset'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m index a836e12a..09d78caa 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2b_basic_image_visualization.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_2b_basic_image_visualization -%CANLAB_TEST_WALKTHROUGH_2B_BASIC_IMAGE_VISUALIZATION Smoke test of canlab_help_2b_basic_image_visualization. +%CANLAB_TEST_WALKTHROUGH_2B_BASIC_IMAGE_VISUALIZATION Hardened smoke test: Basic image visualization walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2b_basic_image_visualization" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2b_basic_image_visualization.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2b_basic_image_visualization'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m index f1d3a208..17258aae 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_2c_loading_datasets.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_2c_loading_datasets -%CANLAB_TEST_WALKTHROUGH_2C_LOADING_DATASETS Smoke test of canlab_help_2c_loading_datasets. +%CANLAB_TEST_WALKTHROUGH_2C_LOADING_DATASETS Hardened smoke test: Loading datasets walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_2c_loading_datasets" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_2c_loading_datasets.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_2c_loading_datasets'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m index dc36fced..3e0f3b72 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3_voxelwise_t_test.m @@ -1,30 +1,45 @@ function tests = canlab_test_walkthrough_3_voxelwise_t_test -%CANLAB_TEST_WALKTHROUGH_3_VOXELWISE_T_TEST Smoke test of canlab_help_3_voxelwise_t_test_walkthrough. +%CANLAB_TEST_WALKTHROUGH_3_VOXELWISE_T_TEST Hardened smoke test: Voxelwise t-test walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. The WorkingFolderFixture in setup -% provides a clean cwd, which matters here because the walkthrough writes -% NIfTI output files. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_3_voxelwise_t_test_walkthrough" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_3_voxelwise_t_test_walkthrough.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_3_voxelwise_t_test_walkthrough'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m index bbf8135c..8aa9a9be 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_3b_atlases_and_ROI.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_3b_atlases_and_ROI -%CANLAB_TEST_WALKTHROUGH_3B_ATLASES_AND_ROI Smoke test of canlab_help_3b_atlases_and_ROI_analysis. +%CANLAB_TEST_WALKTHROUGH_3B_ATLASES_AND_ROI Hardened smoke test: Atlases and ROI analysis walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_3b_atlases_and_ROI_analysis" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_3b_atlases_and_ROI_analysis.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_3b_atlases_and_ROI_analysis'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m index 7958e2eb..b3c54ebb 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4_masking_and_writing.m @@ -1,30 +1,45 @@ function tests = canlab_test_walkthrough_4_masking_and_writing -%CANLAB_TEST_WALKTHROUGH_4_MASKING_AND_WRITING Smoke test of canlab_help_4_masking_and_writing_nifti_files. +%CANLAB_TEST_WALKTHROUGH_4_MASKING_AND_WRITING Hardened smoke test: Masking and writing NIfTI files walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. The WorkingFolderFixture in setup -% provides a clean cwd, which matters here because the walkthrough writes -% NIfTI output files. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_4_masking_and_writing_nifti_files" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_4_masking_and_writing_nifti_files.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_4_masking_and_writing_nifti_files'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m index cc52a045..83f8d922 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_4b_3D_visualization.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_4b_3D_visualization -%CANLAB_TEST_WALKTHROUGH_4B_3D_VISUALIZATION Smoke test of canlab_help_4b_3D_visualization. +%CANLAB_TEST_WALKTHROUGH_4B_3D_VISUALIZATION Hardened smoke test: 3D visualization walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_4b_3D_visualization" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_4b_3D_visualization.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_4b_3D_visualization'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m index b0d2ca7b..a4dcfc25 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_5_regression.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_5_regression -%CANLAB_TEST_WALKTHROUGH_5_REGRESSION Smoke test of canlab_help_5_regression_walkthrough. +%CANLAB_TEST_WALKTHROUGH_5_REGRESSION Hardened smoke test: Regression walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_5_regression_walkthrough" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_5_regression_walkthrough.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_5_regression_walkthrough'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m index f8d86804..64e9cdb3 100644 --- a/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m +++ b/CanlabCore/Unit_tests/walkthroughs/canlab_test_walkthrough_7_multivariate_prediction.m @@ -1,28 +1,45 @@ function tests = canlab_test_walkthrough_7_multivariate_prediction -%CANLAB_TEST_WALKTHROUGH_7_MULTIVARIATE_PREDICTION Smoke test of canlab_help_7_multivariate_prediction_basics. +%CANLAB_TEST_WALKTHROUGH_7_MULTIVARIATE_PREDICTION Hardened smoke test: Multivariate prediction basics walkthrough. % -% Tier B integration test. See canlab_test_walkthrough_1_installing_tools -% for shared rationale and conventions. +% Tier B integration test. Runs a frozen snapshot of the CANlab_help_examples +% walkthrough "canlab_help_7_multivariate_prediction_basics" through the headless, fault-tolerant +% harness canlab_run_walkthrough_snapshot. The snapshot lives under +% Unit_tests/walkthroughs/private/ (a verbatim copy; refresh it from +% example_help_files/ when the tutorial changes). Reaching the end of the +% script without a genuine error counts as a pass; sections that need a +% display, interactive input, or optional data files are skipped, not failed. +% +% These are SLOW and are skipped by the per-push suite. Run them with +% canlab_run_all_tests('Walkthroughs', 'only') % nightly tier +% See canlab_run_walkthrough_snapshot for the hardening details. tests = functiontests(localfunctions); end -function setup(tc) %#ok<*DEFNU> +function setupOnce(tc) %#ok<*DEFNU> +here = fileparts(mfilename('fullpath')); +tc.TestData.snapshot = fullfile(here, 'private', 'canlab_help_7_multivariate_prediction_basics.m'); +end + + +function setup(tc) +% Clean cwd per run (snapshots write NIfTI/figure files) and headless graphics. tc.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); +tc.TestData.PrevFigVis = get(0, 'DefaultFigureVisible'); +set(0, 'DefaultFigureVisible', 'off'); close all force end function teardown(tc) close all force +if isfield(tc.TestData, 'PrevFigVis') + set(0, 'DefaultFigureVisible', tc.TestData.PrevFigVis); +end end function test_runs_to_completion(tc) -script_name = 'canlab_help_7_multivariate_prediction_basics'; -tc.assumeNotEmpty(which(script_name), ... - 'CANlab_help_examples (example_help_files/) must be on the MATLAB path'); -evalc(script_name); -tc.verifyTrue(true); +canlab_run_walkthrough_snapshot(tc, tc.TestData.snapshot); end diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m new file mode 100644 index 00000000..2a55efe9 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_1_installing_tools.m @@ -0,0 +1,210 @@ +%% Installing CANlab Core Tools + +%% About the CANlab Core Tools repository and CANlab toolboxes +% +% The CANlab Core Tools repository is on https://github.com/canlab/CanlabCore +% It contains core tools for MRI/fMRI/PET analysis from the Cognitive and +% Affective Neuorscience Lab (Tor Wager, PI) and our collaborators. Many of +% these functions are needed to run other toolboxes, e.g., the CAN lab?s +% multilevel mediation and Martin Lindquist?s hemodynamic response estimation +% toolboxes. +% +% The tools include object-oriented tools for doing neuroimaging analysis with +% simple commands and scripts that provide high-level functionality for +% neuroimaging analysis. For example, there is an "fmri_data" object type +% that contains neuroimaging datasets (both PET and fMRI data are ok, despite +% the name). If you have created and object called my_fmri_data_obj, then +% plot(my_fmri_data_obj) will generate a series of plots specific to neuroimaging +% data, including an interactive brain viewer (courtesy of SPM software). +% predict(my_fmri_data_obj) will perform cross-validated multivariate prediction +% of outcomes based on brain data. ica(my_fmri_data_obj) will perform independent +% components analysis on the data, and so forth. +% +% There are a wide range of functions and techniques that we've used in the +% CANlab. Most of these are supported by specific functionality in our +% object-oriented tools. You can download other toolboxes from: +% https://github.com/canlab +% +% And you can find an overview of core object-oriented tool functionality +% and demos/walkthroughs on https://canlab.github.io +% +% <> +% +% + +%% CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> + +%% Dependencies +% +% * Matlab statistics toolbox +% * Matlab signal processing toolbox +% * Statistical Parametric Mapping (SPM) software https://www.fil.ion.ucl.ac.uk/spm/ +% (this is used for image reading/writing and the orthviews function.) +% +% For full functionality, the other toolboxes below are recommended: +% + +%% Installing SPM +% +% *!* Important: You will also need to install spm12 separately, and add it +% to your matlab path. +% The latest version at the time of writing is SPM12, which can be downloaded here: +% +% +% + +%% Quick start instructions +% +% # Download SPM12 from https://www.fil.ion.ucl.ac.uk/spm/ +% # Find the SPM12 folder and drag it into your command window. +% Matlab will go to that directory (execute the cd command). +% # In Matlab, type |addpath pwd|. This adds the main SPM12 folder to your path. +% # Type |addpath canonical|. This adds the |canonical| subfolder to your +% path. +% # Go to https://github.com/canlab/ and click on CanlabCore +% # Sign into Github (top right corner) +% # Click "Clone or Download" +% # Click "Open in Desktop" (if you have Github Desktop) to clone the +% repository +% # Find the folder and drag it into your Matlab command window. +% # In Matlab, type |g = genpath(pwd)|. This lists all the subfolders in a string |g| +% # Type |addpath(g)|. This adds them to your path. +% # Type |savepath|. This saves your path for future use. +% +% Repeat the Clone/genpath/addpath steps for the Neuroimaging_Pattern_Masks +% repository, and save the path again. + + +%% +% Quick start: Batch install script +% +% The main goal of the code below is to run: +% |canlab_toolbox_setup| +% This is a script that attempts to: +% # Check for CANlab repositories on your computer +% # Download any that are missing +% # Add them to your Matlab path with subfolders. +% +% _Note_: Unfortunately as of Jan 2020 this does not work well on +% everyone's computers, due to OS-related variability. The workaround is to +% install them manually as per the instructions above. +% +% |canlab_toolbox_setup| is in the CanlabCore repository. So before you run it +% the CANlab_Core_Tools must be added to your path with +% subfolders. Otherwise, you will get errors. +% +% The script canlab_toolbox_setup can help download and install the +% toolboxes you need. +% +% First, go to a folder where you want to install +% toolboxes/repositories on your hard drive. Mine are in "Github" and I've +% already installed CanlabCore, so let's find it and go there: + +% Locate the CanlabCore Github files on your local computer: +% ------------------------------------------------------------------------ + +mypath = what(fullfile('CanlabCore', 'CanlabCore')); + +% Check if we've got it +if isempty(mypath) + disp('Download CanlabCore from Github, and go to that folder in Matlab') + disp('by dragging and dropping it from Finder or Explorer into the Matlab Command Window') + return +else + % Add CanlabCore to Matlab path with subfolders + g = genpath(mypath(1).path); + addpath(g); + +end + +mypath = mypath(1).path; + +% Set the folder in which the repositories will be installed +% ------------------------------------------------------------------------ + +% This is the location (folder) in which the repositories will be installed: +base_dir_for_repositories = fileparts(fileparts(mypath)); + +cd(base_dir_for_repositories) +fprintf('\nInstalling repositories in %s\n', base_dir_for_repositories); + +% Find the setup file, canlab_toolbox_setup.m +% ------------------------------------------------------------------------ + +% This is the file +setupfile = fullfile(mypath, 'canlab_toolbox_setup.m'); + +% Make sure we've got it +if ~exist(setupfile, 'file') + disp('Something went wrong. I can''t find canlab_toolbox_setup.m'); + disp('This file should be included in the CanlabCore Github repository.') +end + +% Run the setup script: +% ------------------------------------------------------------------------ +% +% Note: This can be system-dependent and may not work on all computers. +% If it doesn't, you can default to cloning the CANLab repositories and +% adding them with subfolders to your Matlab path. The most important ones +% starting out are CanlabCore and Neuroimaging_Pattern_Masks +% e.g., +% |git clone https://github.com/canlab/CanlabCore.git| +% |git clone https://github.com/canlab/Neuroimaging_Pattern_Masks.git| + +canlab_toolbox_setup + +% This will attempt to locate toolboxes, add them to your path, and give +% you the option to download them from Github if it can't find them. + +% It looks for and installs toolboxes in the current directory, +% which (thanks to the code above) is the path name in base_dir_for_repositories + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m new file mode 100644 index 00000000..8c8b169f --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2_load_a_sample_dataset.m @@ -0,0 +1,199 @@ +%% Load and explore a sample dataset +% This example shows how to load a pre-set example fMRI dataset into +% an fmri_data object. + +%% General instructions +% +% Before you start, the CANlab_Core_Tools must be added to your path with +% subfolders. Otherwise, you will get errors. +% +% The script canlab_toolbox_setup can help download and install the +% toolboxes you need. +% +% Sample datasets are in the "Sample_datasets" folder in CANlab_Core_Tools. +% Many tutorials apply pre-trained patterns and masks. +% These are stored in this Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% This example will use emotion regulation data in the folder: +% "Wager_et_al_2008_Neuron_EmotionReg" +% The dataset is a series of contrast images from N = 30 participants. +% Each image is a contrast image for [reappraise neg vs. look neg] +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% Here are a couple of helpful functions we will use for display: +% (you can ignore these.) +dashes = '----------------------------------------------'; +printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes); + +%% The fmri_data object +% +% The fmri_data object class is one of the most important, basic types of +% objects in the CANlab object-oriented toolbox. It stores image data in a +% matrix form, which is more space efficient and friendly for analysis with +% various software packages/algorithms. +% +% For philosophy, see: +% +% +% For more info on the object-oriented approach, see: +% +% +% When you call the _class constructor_ fmri_data(), this is what it does: +% +% <> + +%% Section 1: The quick and easy way to load a pre-specified dataset +% +% The function load_image_set has a number of pre-defined image sets that +% you can load with one simple command. This is the easiest way to load +% sample data. The images must be on your Matlab path for this to work. + +[data_obj, subject_names, image_names] = load_image_set('emotionreg'); + +% data_obj is an fmri_data object containing all 30 images. +% subject_names is a list of short names for each image. +% image_names is a list of the full image names with their path. + +%% Section 2: Manual load. Use filenames to find file names +% +% First, we need to list the file names in a string matrix. +% Then, we can load them into an fmri_data object +% We will use the filenames function to get the names, and the fmri_data +% object constructor function to load the data. + +% First, check whether images are on your path: +% We will search for one image and save the path name. + +printhdr('Check that we can find data images:') +myfile = which('con_00810001.img'); +mydir = fileparts(myfile); +if isempty(mydir), disp('Uh-oh! I can''t find the data.'), else disp('Data found.'), end + +% Now we can list all the file names. + +printhdr('Find files and get their names:') +image_names = filenames(fullfile(mydir, 'con_008100*img'), 'absolute'); +disp('Done.'); + +% Now load them into an fmri_data object. + +printhdr('Loading the image data into the object:') +data_obj = fmri_data(image_names); + +% This is the gateway to doing many other things, which are explained in +% other help files. But just to get us started, let's run through a few +% basic things we can do. We'll mainly just look at some standard plots of +% the data. + +%% Section 3: Plot the data we just loaded +% +% Operations that we can perform on fmri_data objects are called methods. +% You can see a list of methods by typing methods(data_obj). +% Here, we'll call the plot method to visualize the data. + +plot(data_obj) +snapnow + +%% Section 4: Get help on an object +% +% Objects have associated help files. +% Type: + +help fmri_data + +% ...to get help for the fmri_data() class constructor and object. +% +% Some objects also have doc files that contain additional information. +% Try: + +doc fmri_data + +%% Section 5: Explore methods +% +% Objects have methods associated with them, or things you can do with +% them. To see a list of methods for the fmri_data object, type: + +methods(fmri_data) + +%% +% You can also use the name of an _instance_ of an object, i.e., a variable +% in the workspace. |methods(data_obj)| produces the same list. +% In addition, you can get help on any of the object's methods by typing: +% |object class name_dot_method name|. e.g., + +help fmri_data.mean + +% Note: There may be trouble accessing these help files in Matlab R2019b. +% Works in 2018 and earlier. But maybe this is not a general problem...give +% it a try. + +%% Section 6: Explore on your own +% +% 1. Try to run a couple of other methods on your fmri_data object. Some +% require additional inputs, but many do not. Look at the help for more +% info. What do you see? Can you return any meaningful output, and if so, +% what did the method do? +% + +% That's it for this section!! + +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m new file mode 100644 index 00000000..b257361c --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2b_basic_image_visualization.m @@ -0,0 +1,149 @@ +%% Basic Image Visualization + +%% +% Using Basic Plot Methods (plot, orthviews, and montage) to Examine a Dataset +% +% Objects can have methods with intuitive names, some of which overlap +% with names of functions in Matlab or other toolboxes. fmri_data.plot() +% is one of these. When you call plot() and pass in an fmri_data object, +% you invoke the fmri_data object method, and a special plot for fmri_data +% objects is produced. +% +% You can list methods for an object class (e.g., fmri_data) by typing: +% +% |methods(fmri_data)| +% +% You can get help for a method by typing +% |help .| +% e.g., |help fmri_data.plot| +% +% The plot() method takes an fMRI data object as a parameter and produces +% an SPM Orthview presentation along with 6 plots of the data. +% +% The plot method expects an fMRI data object to be passed in. We can +% create an fMRI data object using the emotion regulation dataset +% via the following code: + +[image_obj, networknames, imagenames] = load_image_set('emotionreg', 'noverbose'); + +%% +% Once created, we can pass this data object to the plot function to get +% the entire set of outputs, including Matlab console output regarding +% outliers and corresponding data visualizations, using the simple command: + +plot(image_obj); + +%% Explaining the output +% +% The help file for fmri_data.plot did object has more information about the plots: +% +% e.g., |help fmri_data.plot| +% +% e.g., |help image_obj.plot| + +help fmri_data.plot + +%% Use Orthviews to visualize the mean image +% +% Some other commonly used methods to display images are the orthviews() +% and montage() methods. fmri_data objects have these methods, and other +% object classes do to, like the statistic_image, atlas, and region +% classes. +% +% First, let's create a mean image for the dataset: + +m = mean(image_obj); +orthviews(m); +drawnow, snapnow + +%% Threshold and display +% let's threshold it using the threshold() method. +% Here we'll exclude values between -1 and 1 and view the extreme values: + +m = threshold(m, [-1 1], 'raw-outside'); +orthviews(m); +drawnow, snapnow + +%% Create and display a region object +% We can create a region class object, another type of object, from the +% thresholded image. This has additional info and options about each +% contiguous 'blob' in the suprathreshold map: + +r = region(m); +orthviews(r); +drawnow, snapnow + +%% Orthviews options +% Orthviews methods have a range of options. +% They use the cluster_orthviews function, which uses spm_orthviews +% See |help cluster_orthviews| for options. +% +% Let's try one that visualizes each contiguous blob in a different, solid +% color: + +orthviews(r, 'unique', 'solid'); +drawnow, snapnow + +%% Use montage to visualize the thresholded mean image +% +% Sometimes, we want to view map that shows a canonical range of slices. +% This is really useful for producing standard output for papers +% Arguably, one should *always* view and publish montage maps showing all +% slices, so as to show the "whole picture" and not omit any results. +% +% You can customize this a lot, as it uses the fmridisplay() object class, +% which allows you to add custom montages (in axial, saggital, and coronal +% orientations, add blobs of various types, and remove them and re-plot. +% See |help fmridisplay| and |help fmridisplay.addblobs| for more details. +% +% For now, we'll just stick to a basic plot. We'll first create an empty +% figure,then plot the montage on it. + +create_figure('montage'); axis off; +montage(m); +drawnow, snapnow + +% We've already thresholded it, so it'll use the previous threshold. +% however, we can re-threshold the image and redisplay it as well. + +%% Use montage to visualize each blob in a thresholded map +% +% A really useful thing to do is to take a region object, often from a +% thresholded map, and visualize each region. +% the montage() methods also have a number of options. +% Let's try one, |'regioncenters'|, that plots each blob (region) on a +% separate slice. +% +% Furthermore, we can use the |'colormap'| option to view the regions with +% colors mapped to their associated values, e.g., hot colors for positive +% values and cool colors for negative values. +% +% Finally, we might want to assign names to each region based on an atlas. +% We'll do that before plotting, so that the names appear on the plots. +% These names are saved in the r(i).shorttitle field, for each region i. +% The region.table() method automatically labels them as well. +% You can customize the atlas used; the default is the 'canlab2018_2mm' +% atlas (see |help load_atlas| for more info.) + +r = autolabel_regions_using_atlas(r); + +montage(r, 'regioncenters', 'colormap'); +drawnow, snapnow + +% Note that some large regions may span multiple areas. +% This can happen if the various regions are connected by contiguous +% suprathreshold voxels. + +%% Explore on your own +% +% 1. Try to re-threshold the image using some values you choose and re-plot. +% Look at the help for more options on thresholding. and pick one. What do you see? +% Try to plot only voxels with positive values in the mean image. +% +% 2. Try to bring up only *one* region from the region object (r) in orthviews. +% Can you visualize it in all three views? +% +% 3. Try using a couple of other display options in the montage() and orthviews() +% methods. What do you see? + +% That's it for this section!! \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m new file mode 100644 index 00000000..11aed9bd --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_2c_loading_datasets.m @@ -0,0 +1,276 @@ +%% Datasets used in CANlab tutorials +% +% Note: this report was generated from |canlab_help_2c_loading_datasets.m| +% in the repository +% + +%% The Neuroimaging_pattern_masks Github repository and website +% +% For these walkthroughs, we'll load several image datasets. +% Some are stored in Github, if the files are small enough. +% The Wager 2008 emotion regulation dataset, for example, is +% in the CANlab Core repository. Others are on figshare or Neurovault. +% +% Many tutorials apply pre-trained patterns and masks. +% These are stored in this Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% To run this tutorial and many others in this series, you'll need to +% download the Neuroimaging_Pattern_Masks Github repository and add it to +% your Matlab path with subfolders. +% The Github site has three main types of datasets, shown here: +% +% <> +% +%% Registries for easy loading +% There are several functions that contain sets of images that you can load +% by name. For example: +% +% |load_image_set()| includes a registry of datasets that you can access by +% name. Try |help load_image|set| for a list of images you can load by +% keyword. +% +% |load_atlas()| also has a named set of atlases/brain parcellations you can load. +% +% |canlab_load_ROI| has a registry of many named regions derived from previous +% studies. This is particularly useful for loading a subcortical ROI and +% visualizing it or applying it to data. + +%% Emotion regulation dataset +% "Wager_et_al_2008_Neuron_EmotionReg" +% The dataset is a series of contrast images from N = 30 participants. +% Each image is a contrast image for [reappraise neg vs. look neg] +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% To load, try the following: + +[data_obj, subject_names, image_names] = load_image_set('emotionreg', 'noverbose'); + +montage(mean(data_obj), 'trans', 'compact2'); +drawnow, snapnow + +% data_obj is an fmri_data object containing all 30 images. +% subject_names is a list of short names for each image. +% image_names is a list of the full image names with their path. + +%% Buckner lab resting-state connectivity maps +% +% A popular series of masks is a set of "networks" developed from resting-state fMRI +% connectivity in 1,000 people. For our purposes, the "network" maps +% consist of a set of voxels that load most highly on each of 7 ICA +% components from the 1,000 Functional Connectomes project. +% These were published by Randy Buckner's lab in 2011 in three papers: +% Choi et al., Buckner et al., and Yeo et al. We'll use the cortical +% networks from Yeo et al. here. +% + +[bucknermaps, mapnames] = load_image_set('bucknerlab', 'noverbose'); % loads 7 masks from Yeo et al. +disp(char(mapnames)) + +% Create a montage plot +o2 = montage(get_wh_image(bucknermaps, 1), 'trans', 'compact2', 'color', rand(1, 3)); +for i = 2:7, o2 = addblobs(o2, region(get_wh_image(bucknermaps, i)), 'trans', 'color', rand(1, 3)); end +drawnow, snapnow + +%% Pain, Cognition, Emotion balanced N = 270 dataset +% +% This dataset is an fmri_data object class object created using CANlab tools (canlab.github.io). +% It contains 270 single-participant fMRI contrast maps across 18 studies (with 15 participants per study). +% Studies are grouped into three domains: Pain, Cognitive Control, and Negative Emotion, with 9 studies each. +% Each domain includes 3 subdomains, with 3 studies each. +% +% The dataset is from Kragel et al. 2018, Nature Neuroscience. +% It is too large for Github, and it's stored on Neurovault.org +% as collection #504. You could get it using CANlab tools like this: +% +% |[files_on_disk, url_on_neurovault, mycollection, myimages] = retrieve_neurovault_collection(504);| +% +% It is also available on Figshare, with DOI 10.6084/m9.figshare.24033402 +% https://figshare.com/ndownloader/files/42143352 +% +% If you download from Neurovault, you'd +% have to add metadata for the study category labels to use it. +% Therefore, we suggest you use the CANLab load_image_set() function, as below. +% It saves a file on your hard drive the first time you run it: +% +% |kragel_2018_nat_neurosci_270_subjects_test_images.mat| +% +% This file also has metadata that is not necessarily on Neurovault. +% If you save this file somewhere on your Matlab path, you'll be able to +% reload and reuse the dataset easily. + +[test_images, names] = load_image_set('kragel18_alldata', 'noverbose'); + +% This field contains a table object with metadata for each image: +metadata = test_images.metadata_table; + +disp('Metadata available in test_images.metadata_table:') +metadata(1:5, :) + +% Make a plot of the spatial correlation of the average image for each study + +imgs = cellstr(test_images.image_names); +m = mean(test_images); + +for i = 1:length(names) + + % Create a mean image for each study and store in "m" object. + + wh = metadata.Studynumber == i; + studymean = mean(get_wh_image(test_images, find(wh))); + m.dat(:, i) = studymean.dat; + +end + +disp('Map of spatial correlations across the mean images for each study'); + +plot_correlation_matrix(m.dat, 'names', names, 'partitions', [ones(1, 6) 2*ones(1, 6) 3*ones(1, 6)], 'partitionlabels', {'Pain', 'Cognition', 'Emotion'}); +drawnow, snapnow + +%% Pain across six intensity levels per person (BMRK3) +% +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). Each image is the average over several +% (4-8) trials of heat delivered at a single stimulus intensity, ranging +% from 44.3 - 49.3 degrees C in one-degree increments. Each image is also +% paired with an average reported pain value for that set of trials, rated +% immmediately after heat experience. +% +% This dataset is interesting for mixed-effects and predictive analyses, as +% it has both within-person and between-person sources of variance. +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj, an fmri_data class object. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + +subject_id = image_obj.additional_info.subject_id; +ratings = reshape(image_obj.Y, 6, 33)'; +temperatures = image_obj.additional_info.temperatures; + +% Plot the ratings + +create_figure('ratings'); +title('Pain ratings by stimulus intensity in BMRK3 dataset') +hold on; plot(ratings', '-', 'Color', [.7 .7 .7], 'LineWidth', .5); +lineplot_columns(ratings, 'color', [.7 .3 .3], 'markerfacecolor', [1 .5 0]); +xlabel('Temperature'); +ylabel('Rating'); +set(gca, 'XTickLabel', [44.3:49.3], 'FontSize', 18); +drawnow, snapnow + +%% CANlab 2018 Combined Atlas object +% An atlas-class object is a specialized subclass of fmri_data that stores +% information about a series of parcels, or pre-defined regions, and in +% some cases the probalistic maps that underlie the final parcellation. +% +% A list of pre-defined atlases is contained in the function |load_atlas|. +% The default atlas for some CANlab functions is the "CANlab combined 2018" +% This was created by Tor Wager from other published atlases. It includes: +% +% * Glasser 2016 Nature 180-region cortical parcellation (in MNI space, not indivdidualized) +% * Pauli 2016 PNAS basal ganglia subregions +% * Amygdala/hippocampal and basal forebrain regions from SPM Anatomy Toolbox +% * Thalamus regions from the Morel thalamus atlas +% * Subthalamus/Basal forebrain regions from Pauli "reinforcement learning" HCP atlas +% * Diedrichsen cerebellar atlas regions +% * Multiple named brainstem nuclei localized based on individual papers +% * Additional Shen atlas parcels to cover areas (esp. brainstem) not otherwise named +% +% References for the corresponding papers are stored in the atlas object, and +% printed in tables generated with the region.table() method. + +atlas_obj = load_atlas('canlab2018'); + +% Create a figure showing cortical parcels +create_figure('iso'); isosurface(atlas_obj, 'alpha', .5); +drawnow, snapnow + +% Create a figure of the thalamus regions only +thal = select_atlas_subset(atlas_obj, {'Thal'}); +create_figure('iso'); isosurface(thal); +hh = addbrain('hires'); +set(hh, 'FaceAlpha', 0.05); +set(gca, 'XLim', [-40 40], 'YLim', [-40 20]); +drawnow, snapnow + + +%% Kragel 2015 emotion category-predictive patterns +% +% This loads a series of 7 multivariate patterns from Kragel et al. 2015. +% These patterns were developed to predict emotion categories. +% For our purposes here, we'll just treat them as data images. + +[test_dat, names] = load_image_set('kragelemotion', 'noverbose'); % loads 7 masks from Kragel et al. +disp(char(names)) + +% Plot spatial correlations +plot_correlation_matrix(test_dat.dat, 'names', names); +drawnow, snapnow + +%% Subcortical and brainstem regions +% |canlab_load_ROI| has a registry of many named regions derived from previous +% studies. This is particularly useful for loading a subcortical ROI and +% visualizing it or applying it to data. + +help canlab_load_ROI + +vmpfc = canlab_load_ROI('vmpfc'); +nac = canlab_load_ROI('nac'); +pag = canlab_load_ROI('pag'); + +create_figure('brain'); +isosurface(vmpfc, 'color', {[1 0 1]}); +isosurface(nac, 'color', {[1 .5 0]}); +isosurface(pag, 'color', {[1 0 0]}); +bhan = addbrain('hires right'); +view(-119, 5); +drawnow, snapnow + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m new file mode 100644 index 00000000..c3658b4a --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3_voxelwise_t_test_walkthrough.m @@ -0,0 +1,224 @@ +%% Perform a voxel-wise t-test +% +% In this example we perform a second level analysis on first level +% statistical parametric maps. Specifically we use a t-test to obtain an +% ordinary least squares estimate for the group level parameter. +% +% The example uses the emotion regulation data provided with +% CANlab_Core_Tools. This dataset consists of a set of contrast images for +% 30 subjects from a first level analysis. The contrast is [reappraise neg vs. look neg], +% reappraisal of negative images vs. looking at matched negative images. +% +% These data were published in: +% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., +% Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating +% successful emotion regulation. Neuron, 59, 1037-50. +% +% By the end of this example we will have regenerated the results of figure +% 2A of this paper. + +%% Executive summary of the whole analysis +% +% As a summary, here is a complete set of commands to load the data and run +% the entire analysis: +% +% img_obj = load_image_set('emotionreg'); % Load a dataset +% t = ttest(img_obj); % Do a group t-test +% t = threshold(t, .05, 'fdr', 'k', 10); % Threshold with FDR q < .05 and extent threshold of 10 contiguous voxels +% +% % Show regions and print a table with labeled regions: +% montage(t); drawnow, snapnow; % Show results on a slice display +% r = region(t); % Turn t-map into a region object with one element per contig region +% table(r); % Print a table of results using new region names +% +% Here a graphical description of what it does: +% +% <> +% +% Now, let's walk through it step by step, with a few minor additions. + +%% Load sample data + +% We can load this data using load_image_set(), which produces an fmri_data +% object. Data loading exceeds the scope of this tutorial, but a more +% indepth demonstration may be provided by load_a_sample_dataset.m + +[image_obj, networknames, imagenames] = load_image_set('emotionreg'); + +%% Do a t-test + +% Voxel-wise tests estimate parameters for each voxel independently. The +% fmri_data/ttest function performs this just like the native matlab +% ttest() function, and similarly returns a p-value and t-stat. Unlike the +% matlab ttest() function, fmri_data/ttest performs this for every voxel in +% the fmri_data object. All voxels are evaluated independently, but all are +% evaluated nonetheless. A statistic_image is returned. + +t = ttest(image_obj); + +%% Visualize the results +% There are many options. See methods(statistic_image) and methods(region) + +orthviews(t) +drawnow, snapnow; + +% As can be seen, the ttest() result is unthresholded. The threshold +% function can be used to apply a desired alpha level using any of a number +% of methods. Here FDR is used to control for alpha=0.05. Note that no +% information is erased when performing thresholding on a statistic_image. + +t = threshold(t, .05, 'fdr'); +orthviews(t) +drawnow, snapnow; + +% Many neuroimaging packages (e.g., SPM and FSL) do one-tailed tests +% (with one-tailed p-values) and only show you positive effects +% (i.e., relative activations, not relative deactivations). +% All the CANlab tools do two-sided tests, report two-tailed p-values. +% By default, hot colors (orange/yellow) will show activations, and cool +% colors (blues) show deactivations. + +% We can also apply an "extent threshold" of > 10 contiguous voxels: + +t = threshold(t, .05, 'fdr', 'k', 10); +orthviews(t) +drawnow, snapnow; + +% montage is another visualization method. This function may require a +% relatively large amount of memory, depending on the resolution of the +% anatomical underlay image you use. We recommend around 8GB of free memory. +% +create_figure('montage'); axis off +montage(t) +drawnow, snapnow; + +%% Print a table of results + +% First, we'll have to convert to another object type, a "region object". +% This object groups voxels together into "blobs" (often of contiguous +% voxels). It does many things that other object types do, and inter-operates +% with them closely. See methods(region) for more info. + +% Create a region object called "r", with contiguous blobs from the +% thresholded t-map: + +r = region(t); + +% Print the table: + +table(r); + +%% More on getting help + +% Help info is available on https://canlab.github.io/ +% and in the CANlab help examples repository: https://github.com/canlab/CANlab_help_examples +% This contains a series of walkthroughs, including this one. + +% You can get help for the main object types by typing this in Matlab: +% doc fmri_data (or doc any other object types) + +% Also, more detailed help is available for each function using the "help" command. +% You can get help and options for any object method, like "table". But +% because the method names are simple and often overlap with other Matlab functions +% and toolboxes (this is OK for objects!), you will often want to specify +% the object type as well, as follows: + +% help region.table + +%% Write the t-map to disk + +% Now we need to save our results. You can save the objects in your +% workspace or you can write your resulting thresholded map to an analyze +% file. The latter may be useful for generating surface projections using +% Caret or FreeSurfer for instance. +% +% Thresholding did not actually eliminate nonsignificant voxels from our +% statistic_image object (t). If we simply write out that object, we will +% get t-statistics for all voxels. + +t.fullpath = fullfile(pwd, 'example_t_image.nii'); +write(t, 'overwrite') + +% If we use the 'thresh' option, we'll write thresholded values: +% (write() does not overwrite an existing file unless asked, and we just +% wrote this path above, so pass 'overwrite' here too.) +write(t, 'thresh', 'overwrite') + +t_reloaded = statistic_image(t.fullpath, 'type', 'generic'); +orthviews(t_reloaded) + + +%% Some useful optional extensions + +figure; surface(t); +figure; surface(t, 'foursurfaces'); +figure; surface(t, 'coronal_slabs_5'); + +table_of_atlas_regions_covered(t); + + +%% Explore on your own +% +% 1. Click around the thresholded statistic image. What lobes of the brain +% are most of the results in? What can we say about areas that are not +% active, if anything? +% +% 2. Why might some of the results appear to be outside the brain? What +% does this mean for the validity of the analysis? Should we consider these +% areas in our interpretation, or should they be "masked out" (excluded)? +% +% That's it for this section!! + +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> +%% +% % Bogdan Petre on 9/14/2016, updated by Tor Wager July 2018, Jan 2020 +% \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m new file mode 100644 index 00000000..4574e9c7 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_3b_atlases_and_ROI_analysis.m @@ -0,0 +1,319 @@ +%% Using the atlas object for Region of Interest (ROI) analyses + +%% About the atlas object +% An atlas-class object is a specialized subclass of fmri_data that stores +% information about a series of parcels, or pre-defined regions, and in +% some cases the probalistic maps that underlie the final parcellation. +% +% Some common uses of atlas objects include: +% * Labeling regions by best-matching atlas regions, as in region.table() +% * Analysis within specific ROIs, or on ROI averages +% * Defining regions for connectivity and graph theoretic analyses +% +% A list of pre-defined atlases is contained in the function |load_atlas|. +% The default atlas for some CANlab functions is the "CANlab combined 2018" +% This was created by Tor Wager from other published atlases. It includes: +% +% * Glasser 2016 Nature 180-region cortical parcellation (in MNI space, not indivdidualized) +% * Pauli 2016 PNAS basal ganglia subregions +% * Amygdala/hippocampal and basal forebrain regions from SPM Anatomy Toolbox +% * Thalamus regions from the Morel thalamus atlas +% * Subthalamus/Basal forebrain regions from Pauli "reinforcement learning" HCP atlas +% * Diedrichsen cerebellar atlas regions +% * Multiple named brainstem nuclei localized based on individual papers +% * Additional Shen atlas parcels to cover areas (esp. brainstem) not otherwise named +% +% References for the corresponding papers are stored in the atlas object, and +% printed in tables generated with the region.table() method. +% +% There are many methods for atlas, including all the fmri_data methods. +% You can see those with |methods(atlas)|. For example: +% +% Using the |select_atlas_subset()| method: +% You can select any subset of atlas regions by name or number, and return +% them in a new atlas object. You can also 'flatten' regions, combining +% them into a single new region. +% +% % Using the |extract_data()| method: +% You can extract average data from every image in an atlas, for a series of target images. +% +% We will explore these here. + +%% load an atlas +% 'atlas' objects are a class of objects specially designed for brain +% atlases. Here is more information on this class (also try |doc atlas|) + +help atlas + +% The function load_atlas in the CANlab toolbox loads a number of named +% atlases included with the toolbox. Here is a list of named atlases: + +help load_atlas + +% Now load the "CANlab combined 2018" atlas: +atlas_obj = load_atlas('canlab2018_2mm'); + + +%% visualize the atlas regions + +orthviews(atlas_obj); + +o2 = montage(atlas_obj); + +drawnow, snapnow +%% Selecting regions of interest +% There are several ways of selecting regions of interest. You can do it: +% * By name or part of name (in the .labels, .labels_2, or other fields) +% * By number +% * Graphically, or by entering a coordinate and a search radius + +%% +% *Select regions graphically or by coordinate* + +% First bring up the orthviews display: +orthviews(atlas_obj); + +% Now click on a coordinate in the center of vmPFC +% We'll use the SPM code below to do this here: + +spm_orthviews('Reposition', [0 38 -11]) + +% Running this will create a subset atlas containing only regions with centers +% within 20 mm of that coordinate. + +[obj_within_20mm,index_vals] = select_regions_near_crosshairs(atlas_obj, 'thresh', 20); + +orthviews(obj_within_20mm); +drawnow, snapnow + +% To reproduce this later in a script, you can enter the coordinates as +% input instead: + +[obj_within_20mm,index_vals] = select_regions_near_crosshairs(atlas_obj, 'coords', [0 38 -11], 'thresh', 20); + +% This function returns index values, so you can use these directly to +% specify the regions in select_atlas_subset as well. + +find(index_vals) + +test_obj = select_atlas_subset(atlas_obj, find(index_vals)); + +% % This should yield the same map of vmPFC regions: orthviews(test_obj) + +%% +% *Select regions by name* +% Let's select all the regions in the thalamus. All regions are labeled in +% the atlas object, so we can select them by name. + +% Select all regions with "Thal" in the label: +thal = select_atlas_subset(atlas_obj, {'Thal'}) + +% Print the labels: +thal.labels + +% Select a few thalamus/epithalamus regions of interest: +thal = select_atlas_subset(atlas_obj, {'Thal_Intra', 'Thal_VL', 'Thal_MD', 'Thal_LGN', 'Thal_MGN', 'Thal_Hb'}); +thal.labels + +% Select all the regions with "Thal" in the label, and collapse them into a single region: +whole_thal = select_atlas_subset(atlas_obj, {'Thal'}, 'flatten'); + +%% Extract and analyze data from atlas regions +% + +%% +% First, we'll load a dataset to extract data from for an ROI analysis +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj +% This is an fmri_data object from the CANlab Core Tools repository for neuroimaging data analysis. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed +% +% Alternative sample datasets: +% -------------------------------------------- +% This dataset will take time to download from figshare if you haven't yet +% downloaded it. An alternative is to use the dataset included with the +% CANlab Core toolbox (though the results are less interesting for this +% example). +% image_obj = load_image_set('emotionreg'); + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + + +%% +% *Now, we'll extract data from each atlas region* +% - "r" is a region-class object. see >> help region +% - r(i).dat contains averages over voxels within region i. for n images, it contains an n x 1 vector with average data. +% - r(i).data contains an images x voxels matrix of all data within region i. + +r = extract_roi_averages(image_obj, thal); + +%% +% You could also extract data from the whole thalamus into a single +% variable using: +% |r = extract_roi_averages(data_obj, whole_thal);| + +%% +% An alternative is the extract_data method: +% |region_avgs = extract_data(thal, image_obj);| +% The values may differ slightly based on which image is being resampled to +% which. Resampling matches the voxels in two sets of images that may have different +% voxel sizes and voxel boundaries by interpolating values from one to +% match the space of the other. + +%% Do a group ROI analysis and make a "violin plot" (or barplot) for each region +% +% We'll use the |barplot_columns| function, which is a versatile function +% for plotting columns of matrices. It has many options for displaying or +% hiding individual data points, "violins" (distribution estimates), +% standard errors, and more. +% +% The |barplot_columns| function also does a basic one-sample t-test +% analysis on each column, which corresponds to the average of each region +% here. This ROI analysis can be used for inference about a region. +% +% With the pain dataset, for a real analysis we might worry that multiple +% images come from the same subject. We can't do stats across all the +% iamges without accounting for the correlations among images belonging to +% the same subject. But first, let's plot the averages and see waht we get: + +% Concatenate the region averages into a matrix: +roi_averages = cat(2, r.dat); + +% Get the labels for each region from the atlas object, \ +% and replace some characters to make them look a bit nicer: +thal_labels = format_strings_for_legend(thal.labels); + +% Now make the figure: +create_figure('Thalamus regions', 1, 2); + +barplot_columns(roi_averages, 'nofig', 'colors', scn_standard_colors(length(r)), 'names', thal_labels); +xlabel('Region') +ylabel('ROI activity') + +subplot(1, 2, 2) + +barplot_columns(roi_averages, 'nofig', 'colors', scn_standard_colors(length(r)), 'names', thal_labels, 'noind', 'noviolin'); +xlabel('Region') +ylabel('ROI activity') + +%% +% Look at the results. Do they make sense for a pain dataset? +% The sensory/nociceptive regions of the thalamus (VL and intralaminar) respond to painful +% stimuli, but not the regions in visual and auditory pathways (LGN, MGN). The +% habenula also responds. + +%% +% To do a proper ROI analysis, let's just select a subset of the images +% responding to one stimulus intensity, with one image per person. +% the stimulus temperatures are stored in +% |image_obj.additional_info.temperatures| +% So we can use basic Matlab code to select one temperature. + +wh = image_obj.additional_info.temperatures == 48; % A logical vector for 48 degrees +indx_list = find(wh); % a list of which images + +image_obj_48 = get_wh_image(image_obj, indx_list); % Get the images we want + +% You can try to extract data from the thalamus and repeat the t-test. + +%% Selecting multiple regions using custom label fields +% In addition to the .labels field, atlas objects have fields for additional +% labels (.labels_2, etc.). These can be used to group, and select, regions +% at multiple levels of a spatial hierarchy. For example, ?canlab2018_2mm? +% atlas has macro-network and anatomical group labels in the ?labels_2? field. +% You can use these to select all the regions associated with a particular +% resting-state network or group like the brainstem or basal ganglia. +% Here?s a code snippet to locate/pull out all the default mode network +% regions using the select_atlas_subset() method. This method can also return +% a vector that you can use to select elements of associated connectivity matrices. +% +% You can create and/or annotate your own atlas labels and add them to the +% .labels, labels_2, or other fields, and use select_atlas_subset to +% extract and manipulate them. + +atlas_obj = load_atlas('canlab2018_2mm'); + +% List all the macro-level labels in .labels_2 field +disp(unique(atlas_obj.labels_2)') + +% Pull out all regions with ?Def? (Default) as part of the label in labels_2 +atlas_obj = select_atlas_subset(atlas_obj, {'Def'}, 'labels_2'); +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); +% let's add a cortical surface for context +% We'll make the back surface (right) opaque, and the front (left) transparent +han = addbrain('hires right'); +set(han, 'FaceAlpha', 1); +han2 = addbrain('hires left'); +set(han2, 'FaceAlpha', 0.1); +axis image vis3d off +material dull +view(-88, 31); + +lightRestoreSingle +drawnow, snapnow + +%% Explore on your own +% +% 1. One of the ideas about emotion regulation is that positive appraisal involves +% activity increases in the nucleus accumbens (NaC) in the basal ganglia. +% Load the emotion regulation dataset from previous walkthroughs. +% Try to locate and extract data from this region, and do an ROI analysis. +% What do you see? +% Hint: The NAc is part of the ventral striatum. +% To extract the ventral striatum from the canonical CANlab atlas, try: +% NAc = select_atlas_subset(atlas_obj, {'V_Str'}); +% +% 2. try the t-test on 48 degree only images above. What do you see? Which +% regions show sigificant activation? +% +% 3. Try to identify which prefrontal regions show the greatest emotion +% regulation effect, and which show the greatest response during pain. +% Are they the same? + +% That's it for this section!! + + + + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m new file mode 100644 index 00000000..d05620fd --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4_masking_and_writing_nifti_files.m @@ -0,0 +1,163 @@ +%% Applying masks and writing image data to disk + +% One common operation in neuroimaging is masking. Masking is defining a set +% of voxels to include or exclude from analysis or presentation. Some uses +% include: +% +% * Apply a generous brain mask to store only in-brain voxels, saving space, +% (CANlab data objects do this by default, with a liberal mask) +% * Identify gray matter voxels to calculate global mean signal or other metrics +% * Identify voxels in ventricles or out-of-brain to try to isolate sources +% of noise or artifacts to remove +% * Identify gray matter voxels of interest to limit the space within which +% you correct for multiple comparisons, thus increasing power +% * Limiting presentation of results to an a priori set of ROIs +% * Limiting the scope of training or testing a multivariate pattern or +% predictive model, to help test hypotheses about which systems/regions +% contribute to prediction. +% +%% The Neuroimaging_pattern_masks Github repository and website +% +% For this walkthrough, we'll load some new image datasets to use as masks. +% A popular one is a set of "networks" developed from resting-state fMRI +% connectivity in 1,000 people. For our purposes, the "network" maps +% consist of a set of voxels that load most highly on each of 7 ICA +% components from the 1,000 Functional Connectomes project. +% These were published by Randy Buckner's lab in 2011 in three papers: +% Choi et al., Buckner et al., and Yeo et al. We'll use the cortical +% networks from Yeo et al. here. +% +% This set of images and other broadly useful sets of images are stored in +% the CANlab function |load_image_set| in a registry that you can access by +% name. Try |help load_image|set| for a list of images you can load by +% keyword. These datasets, if they can be stored on Github, are in this +% Github repository: +% +% +% +% In addition, you can explore these and find more information here: +% +% +% To run this tutorial and many others in this series, you'll need to +% download the Neuroimaging_Pattern_Masks Github repository and add it to +% your Matlab path with subfolders. +% The Github site has three main types of datasets, shown here: +% +% <> +% + +%% Load and select a sample mask +% You will need the Neuroimaging_Pattern_Masks repository installed to find +% these images. + +[bucknermaps, mapnames] = load_image_set('bucknerlab'); % loads 7 masks from Yeo et al. + +%% +% This loads 7 images. try |orthviews(dat)| to see them all. +% What values do the images have? +% You can also see descriptive statistics for the images like this: + +desc = descriptives(bucknermaps); + +%% +% The variable |desc| is a structure that contains various descriptive +% stats. +% + +%% +% load_image_set will generally store and return names for each image, +% returned as outputs if requested (here in |mapnames|). Let's print +% those names: + +disp(char(mapnames)) + +%% +% Now we can select one image from those 7 to use as a mask. + +mask_image = bucknermaps.get_wh_image(1); + +%% +% Now we'll load a sample dataset to apply the mask to. +% This loads a series of 7 multivariate patterns from Kragel et al. 2016. +% These patterns were developed to predict emotion categories. +% For our purposes here, we'll just treat them as data images. + +test_dat = load_image_set('kragelemotion', 'noverbose'); % loads 7 masks from Kragel et al. + +% Try |orthviews(dat)| or use |montage| to see these patterns. +% This command shows the first 4 (a limit in |montage| to avoid messy +% displays): + +montage(test_dat); + +%% Apply a mask +% +% We can apply the mask using the apply_mask() method for fmri_data +% objects. It has various options, but the simplest one applies a single +% mask to an fmri_data object. As a rule, you should pass in the object you're +% operating on first, followed by additional modifying arguments/objects. +% Here, this means data object followed by mask object. + +test_dat = apply_mask(test_dat, mask_image); + +% Now let's view the new, masked dataset: + +montage(test_dat); + +% Now, the dataset contains only voxels in the "visual network", and +% subsequent analyses or operations would test only that. + +%% A few exercises +% +% 1. Here's a question to answer: What are some potential uses of the masking +% operation we did above? Why would you want to do it? +% +% 2. Try loading the emotion regulation dataset from previous tutorials, +% and do a one-sample t-test group analysis only within the +% "frontoparietal" network. Try thresholding the results with FDR +% correction for both the full dataset and the analysis within the +% restricted mask. How are the results different? Why? + +%% Writing image data in objects to Nifti files on disk +% +% canlabcore fmri data objects can be easily written to standard image file +% formats. The main function for this is write() +% +% This walkthrough shows how to do this for 1) a mask image, and 2) a +% statistic image. + +%% +% Create a Nifti image called dummy.nii in the current folder. Write the +% mask image to this file. See help for write() for more options + +fname = 'dummy.nii'; % outname file name. the extension (.nii, .img) determines the format. + +% |write(mask_image, 'fname', fname);| + +%% +% If you've already done this, |write()| will return an error. This is a +% safeguard to keep users from inadvertently overwriting important files. +% To force overwriting the existing file, try: + +write(mask_image, 'fname', fname, 'overwrite'); + +%% +% Now let's read it back in from disk and view it again: + +orthviews(fmri_data(fname)) + + +%% Exercises +% +% 1. Select a sample image from |test_dat|, threshold it at an arbitrary threshold +% of absolute values > .0005, and write it to disk. Load it again from +% disk and display it. See the thresholding walkthrough for more help if +% needed. +% Hint: help fmri_data.threshold will give you help on how to do this +% e.g., Retain voxels with absolute value > .0001 +% obj = threshold(obj, [-.0001 0001], 'raw-outside') +% +% 2. Do the one-sample t-test on the emotion regulation dataset (see the +% t-test walkthrough if needed), and threshold the results at 0.05 +% FDR-corrected. Write the thresholded results map to disk. Reload it from +% disk and display it using the montage() method. diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m new file mode 100644 index 00000000..ad99f6bf --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_4b_3D_visualization.m @@ -0,0 +1,468 @@ +%% 3-D volume visualization + +%% About volume visualisation +% +% Creating 3-D renderings of brain data and results can help localize brain +% activity and show patterns in a way that allows meaningful comparisons +% across studies. +% +% The CANlab object-oriented tools have methods built in for rendering +% volume data on brain surfaces. They generally use Matlab's isosurface and +% isocaps functions. +% + +%% Prepare sample statistical results maps +% +% For this walkthrough, we'll use the "Pain, Cognition, Emotion balanced N = 270 dataset" +% from Kragel et al. 2018, Nature Neuroscience. +% More information about this dataset and how to download it is in an +% earlier walkthrough on "loading datasets" +% +% The code below loads it; but you could also use any other brain dataset +% for this purpose (e.g., the 'emotionreg' dataset in load_image_set. +% +% Our goal is to select subject-level images corresponding to a particular task type and +% do a t-test on these, and visualize the results in various 3-D ways. + +[test_images, names] = load_image_set('kragel18_alldata', 'noverbose'); + +% This field contains a table object with metadata for each image: +metadata = test_images.metadata_table; +metadata(1:5, :) % Show the first 5 rows + +% Here are the 3 domains: +disp(unique(metadata.Domain)) + +% Find all the images of "Neg_Emotion" and pull those into a separate +% fmri_data object: + +wh = find(strcmp(metadata.Domain, 'Neg_Emotion')); +neg_emo = get_wh_image(test_images, wh); + +% Do a t-test on these images, and store the results in a statistic_image +% object called "t". Threshold at q < 0.001 FDR. This is an unusual +% threshold, but we have lots of power here, and want to restrict the +% number of significant voxels for display purposes below. + +t = ttest(neg_emo, .001, 'fdr'); + +% Find the images for other conditions and save them for later: + +t_emo = t; + +wh = find(strcmp(metadata.Domain, 'Cog_control')); +cog_control = get_wh_image(test_images, wh); +t_cog = ttest(cog_control, .001, 'fdr'); + +wh = find(strcmp(metadata.Domain, 'Pain')); +pain = get_wh_image(test_images, wh); +t_pain = ttest(pain, .001, 'fdr'); + + +%% The surface() method +% Most data objects, including fmri_data, statistic_image, and region objects, +% have a surface() method. Entering it with no arguments creates a surface +% or series of surfaces (depending on the object) +% +% The surface() method lets you easily render blobs +% on a number of pre-set choices for brain surfaces or 3-D cutaway surfaces. + +% This generates a series of 6 surfaces of different types, including a +% subcortical cutaway +surface(t); +drawnow, snapnow + +% This generates a different plot +create_figure('lateral surfaces'); +surface_handles = surface(t, 'foursurfaces', 'noverbose'); +snapnow + +% You can see more options with: +% |help statistic_image.surface| or |help t.surface| + +% You can also use the |render_on_surface()| to method to change the colormap: + +render_on_surface(t, surface_handles, 'colormap', 'summer'); +snapnow + +%% Using the fmridisplay object to create a registry of montages/surfaces +% The fmridisplay object is an object class that stores information about a +% set of display items you have created. It stores information about a +% series of montages (in the same figure or different figures) and surfaces +% along with their handles. +% +% You can use this to render activations on a whole set of slices/surfaces +% at once. You can create your own custom views, or use pre-packaged sets +% of montages/surfaces created by the function |canlab_results_fmridisplay| +% +% The code snippet below displays our t-map on a relatively complete series +% of slices and surfaces, so you can see the "whole picture". + +my_display_obj = canlab_results_fmridisplay([], 'montagetype', 'full'); +my_display_obj = addblobs(my_display_obj, region(t)); +snapnow + +%% +% You can see some prepackaged sets in |help canlab_results_fmridisplay| +% These include: +% 'full' Axial, coronal, and saggital slices, 4 cortical surfaces +% 'compact' Midline saggital and two rows of axial slices [the default] +% 'compact2' A single row showing midline saggital and axial slices +% 'multirow' A series of 'compact2' displays in one figure for comparing different images/maps side by side +% 'regioncenters' A series of separate axes, each focused on one region +% +% In addition, you can pass in a series of optional keywords that will be passed onto the +% |addblobs| object method that controls rendering of the blobs on slices. +% (They don't all work for surfaces). These control colors, inclusion of an +% *outline* around each activation blob/region, and more. +% See |help fmridisplay.addblobs| for more information. +% Here is an example using a few of the options. + +my_display_obj = canlab_results_fmridisplay(t, 'compact2', 'color', [1 0 0], 'nooutline', 'trans'); + +%% +% After creating the display, other maps can be added. +% The addblobs() method requires a region object, so we use region(t) to convert to a region object. +% This shows us regions activated during negative emotion in red and pain +% in blue. + +my_display_obj = addblobs(my_display_obj, region(t_pain), 'color', [0 0 1], 'nooutline', 'trans'); +snapnow + +%% +% my_display_obj is an fmridisplay-class object. It has its own methods: + +methods(my_display_obj) + +% One particularly useful feature is that you can add or _remove_ activation blobs. +% This allows you to create a canonical set of views and then render one +% activation map on them, remove the blobs, and render other maps on the +% same set of display items. You can add points, spheres (e.g., from +% different studies), and more. +% +% When you operate on an fmridisplay object, pass it back out as an output +% argument so that the information you updated will be available. + +my_display_obj = removeblobs(my_display_obj); +snapnow + +%% Using addbrain to load or build surfaces +% You can also build any surface object you want and render colors onto it. +% This is a very versatile method. +% +% Matlab creates handles to surfaces (and +% other graphics objects). You can use the handles to manipulate the +% objects in a many ways. They have, for example, color, line, face color, +% edge color, and many more properties. If |h| is a figure handle, |get(h)| +% shows a list of its properties. |set(h, 'MyPropertyName', myvalue)| sets +% the property |'MyPropertyName'| to |myvalue|. +% +% Here, we'll use |addbrain| to build up a set of surfaces with handles, +% and then render colored blobs on them using |surface()|. Instead of the +% object method, you can also use its underlying function: |cluster_surf| +% +% Try |help addbrain| for a list of surfaces, and |help cluster_surf| for +% other color/display/scaling options. + +% Let's build a surface by starting with a group of thalamic nuclei, and +% adding the parabrachial complex and the red nucleus. +% Then, we'll render our t-statistics in colors onto those surfaces. + +create_figure('cutaways'); axis off +surface_handles = addbrain('thalamus_group'); + +surface_handles = [surface_handles addbrain('pbn')]; +surface_handles = [surface_handles addbrain('rn')]; +surface_handles = [surface_handles addbrain('pag')]; + +drawnow, snapnow + +%% +% Now render the statistic image stored in _t_ onto those surfaces. +% All non-activated areas turn gray. + +render_on_surface(t, surface_handles, 'colormap', 'summer'); +set(surface_handles, 'FaceAlpha', .85); +set(surface_handles(5:6), 'FaceAlpha', .15); % Make brainstem and thalamus shell more transparent +drawnow, snapnow + +% Emphasize positive values, in a hot colormap +render_on_surface(t, surface_handles, 'clim', [0.01 3]); +drawnow, snapnow + +% Show positive and negative values, in a hot colormap +render_on_surface(t, surface_handles, 'clim', [-3 3]); +drawnow, snapnow + +%% +% Now render the statistic image stored in _t_ onto those surfaces. +% All non-activated areas turn gray. + +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +%% +% Transparent surfaces are a little hard to see. +% We can use Matlab's powerful handle graphics system to inspect and change +% all kinds of properties. See the Matlab documentation for more details. +% Let's just make the surfaces solid: + +set(surface_handles, 'FaceAlpha', 1); +drawnow, snapnow + +%% +% These surfaces can be rotated too (or zoom in/out, pan, set lighting, +% etc.) Let's shift the angle. + +view(222, 15); % rotate +camzoom(0.8); % zoom out + +% A helpful CANlab function to re-set the lighting after rotating is: + +lightRestoreSingle; + +drawnow, snapnow + +%% +% addbrain also has keywords for composites of multiple surfaces + +create_figure('cutaways'); axis off +surface_handles = addbrain('limbic'); +t.surface('surface_handles', surface_handles, 'noverbose'); + +%% Removing blobs and re-adding new colors +% There is a special command to re-set the surface colors to gray. +% Then we can add a t-test for different contrasts without re-drawing +% surfaces. + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Cognitive Control'); +t = ttest(cog_control, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Pain'); +t = ttest(pain, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + +surface_handles = addbrain('eraseblobs',surface_handles); + +title('Negative Emotion'); +t = ttest(neg_emo, .001, 'unc'); +t.surface('surface_handles', surface_handles, 'noverbose'); +drawnow, snapnow + + +%% Prepackaged cutaway surfaces +% Addbrain has many 3-D surfaces, including cutaways that show various 3-D +% sections. The code below visualizes the various options for cutaways. +% +% You can also pass any of these keywords into surface() to generate the +% surface and render colored blobs. +% Let's render these with and without an unthresholded pain map + +t = ttest(pain, 1, 'unc'); +t_thr = threshold(t, .05, 'fdr'); + +keywords = {'left_cutaway' 'right_cutaway' 'left_insula_slab' 'right_insula_slab' 'accumbens_slab' 'coronal_slabs' 'coronal_slabs_4' 'coronal_slabs_5'}; + +for i = 1:length(keywords) + + create_figure('cutaways'); axis off + + surface_handles = surface(t_thr, keywords{i}); + + % This essentially runs the code below: + + % surface_handles = addbrain(keywords{i}, 'noverbose'); + % render_on_surface(t, surface_handles, 'clim', [-7 7]); + + % Alternative: This command creates the same surfaces: + % surface_handles = canlab_canonical_brain_surface_cutaways(keywords{i}); + % render_on_surface(t, surface_handles, 'clim', [-7 7]); + + drawnow, snapnow + +end + +%% +% Plot thresholded and unthresholded maps side by side + +create_figure('cutaways',1, 2); axis off + +surface_handles = surface(t, 'right_cutaway'); +title('Pain, Unthresholded') + +subplot(1, 2, 2); +surface_handles = surface(t_thr, 'right_cutaway'); +title('Pain, thresholded') + +drawnow, snapnow + + +%% +% *Controlling options in render_on_surface()* +% +% When the t map has positive and negative values, render_on_surface creates a special +% bicolor split colormap that has warm colors for positive vals, and cool colors +% for negative vals. You can set the color limits (here, in t-values) + +create_figure('cutaways'); axis off +surface_handles = surface(t_thr, 'coronal_slabs'); + +render_on_surface(t_thr, surface_handles, 'clim', [-4 4]); +drawnow, snapnow + +% You can set the colormap to any Matlab colormap: +render_on_surface(t_thr, surface_handles, 'colormap', 'summer'); +drawnow, snapnow + +% If your image is positive-valued only (or negative-valued only), a +% bicolor split map will not be created: + +t = threshold(t_thr, [2 Inf], 'raw-between'); +render_on_surface(t, surface_handles, 'colormap', 'winter', 'clim', [2 6]); +drawnow, snapnow + +%% Creating isosurfaces +% The fmri_data.isosurface() method allows you to create and save +% isosurfaces for any image, parcel, or blob. +% +% You can also save the meshgrid output and surface structure that lets you +% render this surface easily in the future. +% +% Here's a simple example visualizing all of the parcels in an atlas as 3-D +% blobs: + +t = ttest(pain, .01, 'unc'); + +atlas_obj = load_atlas('thalamus'); + +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); + +% Now we'll set some lighting and figure properties +axis image vis3d off +material dull +view(210, 20); +lightRestoreSingle + +drawnow, snapnow + +% You can render blobs on these surfaces too. + +render_on_surface(t, surface_handles, 'colormap', 'hot'); +drawnow, snapnow + + +%% +% Let's load another atlas, and pull out all the "default mode" regions +% We'll use isosurface to visualize these + +atlas_obj = load_atlas('canlab2018_2mm'); +atlas_obj = select_atlas_subset(atlas_obj, {'Def'}, 'labels_2'); + +create_figure('isosurface'); +surface_handles = isosurface(atlas_obj); + +% let's add a cortical surface for context +% We'll make the back surface (right) opaque, and the front (left) transparent +han = addbrain('hires right'); +set(han, 'FaceAlpha', 1); + +han2 = addbrain('hires left'); +set(han2, 'FaceAlpha', 0.1); + +axis image vis3d off +material dull +view(-88, 31); +lightRestoreSingle + +drawnow, snapnow + +%% Using isosurface to create a custom surface +% You can build a cutaway surface or a series of them with the isosurface method, from +% any suitable image. The image used here is a mean anatomical image that +% renders nicely. +% As always, you can use Matlab's rendering to change the look of the surfaces + +anat = fmri_data(which('keuken_2014_enhanced_for_underlay.img'), 'noverbose'); + +create_figure('cutaways'); axis off + +surface_handles = isosurface(anat, 'thresh', 140, 'nosmooth', 'xlim', [-Inf 0], 'YLim', [-30 Inf], 'ZLim', [-Inf 40]); + +render_on_surface(t, surface_handles, 'colormap', 'hot'); + +view(-127, 33); +set(surface_handles, 'FaceAlpha', .85); +snapnow + +%% Explore on your own +% +% 1. Try to create your own custom brain surface to visualize one of the 3 +% statistic image maps we've been working on. Can you make a plot that +% really shows the important results in a clear way? +% +% 2. Try exploring movie_tools to create a movie where you rotate, zoom, +% and change the transparency of your surfaces in some way. +% You can write movie files to disk for use in presentations, too. +% +% 3. Try rendering another atlas, or subset of an atlas, with isosurface() +% Maybe you can pull out and render all the "ventral attention" network +% regions. +% +%% Explore More: CANlab Toolboxes +% Tutorials, overview, and help: +% +% Toolboxes and image repositories on github: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
CANlab Core Toolshttps://github.com/canlab/CanlabCore
CANlab Neuroimaging_Pattern_Masks repositoryhttps://github.com/canlab/Neuroimaging_Pattern_Masks
CANlab_help_exampleshttps://github.com/canlab/CANlab_help_examples
M3 Multilevel mediation toolboxhttps://github.com/canlab/MediationToolbox
M3 CANlab robust regression toolboxhttps://github.com/canlab/RobustToolbox
M3 MKDA coordinate-based meta-analysis toolboxhttps://github.com/canlab/Canlab_MKDA_MetaAnalysis
+% +% +% Here are some other useful CANlab-associated resources: +% +% +% +% +% +% +% +% +% +% +% +% +% +% +% +%
Paradigms_Public - CANlab experimental paradigmshttps://github.com/canlab/Paradigms_Public
FMRI_simulations - brain movies, effect size/powerhttps://github.com/canlab/FMRI_simulations
CANlab_data_public - Published datasetshttps://github.com/canlab/CANlab_data_public
M3 Neurosynth: Tal Yarkonihttps://github.com/neurosynth/neurosynth
M3 DCC - Martin Lindquist's dynamic correlation tbxhttps://github.com/canlab/Lindquist_Dynamic_Correlation
M3 CanlabScripts - in-lab Matlab/python/bashhttps://github.com/canlab/CanlabScripts
+% +% +% *Object-oriented, interactive approach* +% The core basis for interacting with CANlab tools is through object-oriented framework. +% A simple set of neuroimaging data-specific objects (or _classes_) allows you to perform +% *interactive analysis* using simple commands (called _methods_) that +% operate on objects. +% +% Map of core object classes: +% +% <> \ No newline at end of file diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m new file mode 100644 index 00000000..88c13409 --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_5_regression_walkthrough.m @@ -0,0 +1,561 @@ +%% Regression Walkthrough + +%% +% This walkthrough works through multiple aspects of a basic 2nd-level (group) +% analysis. We are interested in activity related to emotion regulation +% for a contrast comparing reappraisal of negative images vs. looking at +% negative images without reappraisal, i.e., a [Reappraise Negative - Look Negative] +% contrast. This walkthrough will regress brain +% contrast values in each voxel (y) onto individual differences in +% "Reappraisal Success", defined as the drop in negative emotion ratings +% for the [Reappraise Negative - Look Negative] contrast. +% +% We are interested in both the brain correlates of individual differences +% and the group-average [Reappraise Negative - Look Negative] effect in the +% brain, which is the intercept in our regression model. + +%% +% This walkthrough includes multiple considerations you'll need when doing +% a practical data analysis. This goes far beyond just running a model! +% It's also important to understand the data and results. +% Understanding the data allows you to make informed a priori choices about +% the analysis and avoid "P-Hacking". Understanding the results helps +% ensure they are valid, and helps you interpret them. +% +% A) Understanding the data +% +% # Loading the dataset into an fmri_data object +% # Examining the brain data and checking +% +% outliers +% global signal +% which voxels are analyzed (coverage) +% homogeneity across subjects (images) +% +% # Examining the behavioral (individual differences) data and leverages +% # Use the info gained to make data preparation and analysis decisions +% +% Masking +% Transformation of the y variable +% Rescaling/normalizing image data +% Robust vs. OLS estimation +% +% B) Fitting the model +% +% # Running the regression and writing results images to disk +% +% C) Understanding the results +% +% # Explore variations in statistical threshold +% # Compare with other studies by comparing with a meta-analysis mask +% # Explore the impact of variations in data scaling/transformation +% # Localizing results, labeling regions, and making tables using a brain atlas +% # Extract and plot data from regions of interest +% # Explore the impact of plotting data in biased vs. unbiased regions of interest +% +% D) Multivariate prediction (overlaps with later walkthroughs) +% # Extra: Multivariate prediction within a set of unbiased ROIs + +%% Part 1: Load sample data +% Load a sample dataset, including a standard brain mask +% And a meta-analysis mask that provides a priori regions of interest. +% --------------------------------------------------------------- + +% Load sample data using load_image_set(), which produces an fmri_data +% object. Data loading exceeds the scope of this tutorial, but a more +% indepth demosntration may be provided by canlab_help_2_load_a_sample_dataset.m + +% These are [Reappraise - Look Neg] contrast images, one image per person + +[image_obj, networknames, imagenames] = load_image_set('emotionreg'); + +% Summarize and print a table with characteristics of the data: +desc = descriptives(image_obj); + +% Load behavioral data +% This is "Reappraisal success", one score per person, in our example +% If you do not have the file on your path, you will get an error. + +beh = importdata('Wager_2008_emotionreg_behavioral_data.txt') +success = beh.data(:, 2); % Reappraisal success + +% Load a mask that we would like to apply to analysis/results + +mask = which('gray_matter_mask.img') + +maskdat = fmri_data(mask, 'noverbose'); + +% Map of emotion-regulation related regions from neurosynth.org +metaimg = which('emotion regulation_pAgF_z_FDR_0.01_8_14_2015.nii') + +metamask = fmri_data(metaimg, 'noverbose'); + +%% Part 2: Examine the setup and coverage + +% Visualize the mask +% --------------------------------------------------------------- +% BLOCK 2: Check mask + +% This is an underlay brain with three separate montages so we can compare images: + +o2 = canlab_results_fmridisplay([], 'noverbose', 'multirow', 3); +drawnow, snapnow; + +% This is a basic gray-matter mask we will use for analysis: +% It can help reduce multiple comparisons relative to whole-image analysis +% but we should still look at what's happening in ventricles and +% out-of-brain space to check for artifacts. + +o2 = addblobs(o2, region(maskdat), 'wh_montages', 1:2); + +% Add a title to montage 2, the axial slice montage: +[o2, title_han] = title_montage(o2, 2, 'Mask image: Gray matter with border'); +set(title_han, 'FontSize', 18) + +% Visualize summary of brain coverage +% --------------------------------------------------------------- +% Check that we have valid data in most voxels for most subjects +% Always look at which voxels you are analyzing. The answer may surprise +% you! Some software programs mask out voxels in individual analyses, +% which may then be excluded in the group analysis. + +% Show summary of coverage - how many images have non-zero, non-NaN values in each voxel + +o2 = montage(region(desc.coverage_obj), o2, 'trans', 'transvalue', .3, 'wh_montages', 3:4); + +[o2, title_han] = title_montage(o2, 4, 'Coverage: Voxels included in images'); +set(title_han, 'FontSize', 18) + +% Visualize meta-analysis mask +% --------------------------------------------------------------- +o2 = montage(region(metamask), o2, 'colormap', 'wh_montages', 5:6); + +[o2, title_han] = title_montage(o2, 6, 'Meta-analysis mask'); +set(title_han, 'FontSize', 18) + +% There are many other alternatives to the montage - e.g.,: +% orthviews(desc.coverage_obj, 'continuous'); + +%% Part 3: Examine the brain data +% This is a really important step to understand what your images look like! +% Good data tend to produce good results. + +% The method plot() for an fmri_data object shows a number of plots that +% are helpful in diagnosing some problems. + +plot(image_obj); + +% Check histograms of individual subjects for global shifts in contrast values + +% The 'histogram' object method will allow you to examine a series of +% images in one panel each. See help fmri_data.histogram for more options, +% including breakdown by tissue type. + +hist_han = histogram(image_obj, 'byimage', 'color', 'b'); + +% This shows us that some of the images do not have the same mean as the +% others. This is fairly common, as individual subjects can often have +% global artifacts (e.g., task-correlated head motion or outliers) that +% influence the whole contrast image, even when baseline conditions are +% supposed to be "subtracted out". +% +% It suggests that we may want to do an outlier analysis and/or standardize +% the scale of the images. We'll return to this below. + +% If we want to visualize this separated by tissue class, try this: + +create_figure('histogram by tissue type') +hist_han = histogram(image_obj, 'byimage', 'by_tissue_type'); + +% Our previous peek at the data suggested that subjects 6 and 16 have +% whole-brain shifts towards acftivation and deactivation, respectively. + +% This plot shows us that there is a substantial shift in gray matter values +% across the brain for some individuals. +% Also, the means across gray and white matter are correlated. So are the standard deviations. +% We might consider a rescaling or nonparametric procedure (e.g., ranking) +% to standardize the images. + +% Some options are here: +help image_obj.rescale + + +%% Part 4: Examine the behavioral/task predictor(s) +% Examine predictor distribution and leverages +% Leverage is a measure of how much each point influences the regression +% line. The more extreme the predictor value, the higher the leverage. +% Outliers will have very high leverage. High-leverage behavioral observations +% can strongly influence, and sometimes invalidate, an analysis. + +X = scale(success, 1); X(:, end+1) = 1; % A simple design matrix, behavioral predictor + intercept +H = X * inv(X'* X) * X'; % The "Hat Matrix", which produces fits. Diagonals are leverage + +create_figure('levs', 2, 1); +plot(success, 'o', 'MarkerFaceColor', [0 .3 .7], 'LineWidth', 3); +set(gca, 'FontSize', 24); +xlabel('Subject number'); +ylabel('Reappraisal success'); + +subplot(2, 1, 2); +plot(diag(H), 'o', 'MarkerFaceColor', [0 .3 .7], 'LineWidth', 3); +set(gca, 'FontSize', 24); +xlabel('Subject number'); +ylabel('Leverage'); + +drawnow, snapnow + +% The distribution suggests that there are some high-leverage values to +% watch out for. They will influence the results disproportionately at every voxel in the brain! +% Ranking predictors and/or robust regression are good +% mitigation/prevention procedures in many cases. + +%% Part 5: Make data prep and analysis decisions +% Before we've "peeked" at the results, let's make some sensible decisions +% about what to do with the analysis. We'll compare the results to the +% "standard analysis" if we hadn't examined the data at all at end. + +% Choices: + +% 1. Masking: Apply a gray-matter mask. This will limit the area subjected to multiple-comparisons +% correction to a meaningful set. I considered {gray matter, no mask} + +image_obj = apply_mask(image_obj, maskdat); % Mask - analyze gray-matter + +% 2. Image scaling: I considered arbitrary rescaling of images, e.g., by +% the l2norm or z-scoring images. Z-scoring can cause artifactual +% de-activations when we are looking at group effects. L2-norming is +% sensible if there are bad scaling problems... but I'd rather minimize ad +% hoc procedures applied to the data. Ranking each voxel is also sensible, and +% if the predictor is ranked as well, this would implement Spearman's correlations. +% I considered {none, l2norm images, zscoreimages, rankvoxels}, and decided +% on ranking voxels. For the group average (intercept), this is not -- I'll use robust regression to help mitigate issues. + +image_obj = rescale(image_obj, 'rankvoxels'); % rescale images within this mask (relative activation) + +% 3. Predictor values: I considered dropping subject 16, which is a global +% outlier and also has high leverage, or windsorizing values to 3 sd. +% I also considered ranking predictor values, a very sensible choice. +% So I considered {do nothing, drop outlier, windsorize outliers, rank values} +% But robust regression (below) is partially redundant with this. I ultimately decided +% on ranking. + +% This is done below: +% image_obj.X = scale(rankdata(image_obj.X), 1); + +% 4. Regression method: I considered {standard OLS and robust IRLS +% regression}. I chose robust regression because it helps mitigate the +% effects of outliers. However, this is partly redundant with ranking, so +% it would have been sensible to choose ranking OR robust regression. + +%% Part 6: Run the regression analysis and get results +% The regress method takes predictors that are attached in the object's X +% attribute (X stands for design matrix) and regresses each voxel's +% activity (y) on the set of regressors. +% This is a group analysis, in this case correlating brain activity with +% reappraisal success at each voxel. + +% .X must have the same number of observations, n, in an n x k matrix. +% n images is the number of COLUMNS in image_obj.dat + +% mean-center success scores and attach them to image_obj in image_obj.X +image_obj.X = scale(success, 1); + +image_obj.X = scale(rankdata(image_obj.X), 1); + +% runs the regression at each voxel and returns statistic info and creates +% a visual image. regress = multiple regression. +% % Track the time it takes (about 45 sec for robust regression on Tor's 2018 laptop) +tic + +% out = regress(image_obj); % this is what we'd do for standard OLS regression +out = regress(image_obj, 'robust'); + +toc + +% out has statistic_image objects that have information about the betas +% (slopes) b, t-values and p-values (t), degrees of freedom (df), and sigma +% (error variance). The critical one is out.t. +% out = +% +% b: [1x1 statistic_image] +% t: [1x1 statistic_image] +% df: [1x1 fmri_data] +% sigma: [1x1 fmri_data] + +% This is a multiple regression, and there are two output t images, one for +% each regressor. We've only entered one regressor, why two images? The program always +% adds an intercept by default. The intercept is always the last column of the design matrix + +% Image 1 <--- brain contrast values correlated with "success" +% Positive effect: 0 voxels, min p-value: 0.00001192 +% Negative effect: 0 voxels, min p-value: 0.00170529 +% Image 2 FDR q < 0.050 threshold is 0.003193 +% +% Image 2 <--- intercept, because we have mean-centered, this is the +% average group effect (when success = average success). "Reapp - Look" +% contrast in the whole group. +% Positive effect: 3133 voxels, min p-value: 0.00000000 +% Negative effect: 51 voxels, min p-value: 0.00024068 + +% Now let's try thresholding the image at q < .05 fdr-corrected. +t = threshold(out.t, .05, 'fdr'); + +% Select the predictor image we care about: (the 2nd/last image is the intercept) +t = select_one_image(t, 1); + +% ...and display on a series of slices and surfaces +% There are many options for displaying results. +% montage and orthviews methods for statistic_image and region objects provide a good starting point. + +o2 = montage(t, 'trans', 'full'); +o2 = title_montage(o2, 2, 'Behavioral predictor: Reappraisal success (q < .05 FDR)'); +snapnow + +% Or: +% orthviews(t) + + +%% Write the thesholded t-image to disk + +% % t.fullpath = fullfile(pwd, 'Reapp_Success_005_FDR05.nii'); +% % write(t) +% % +% % disp(['Written to disk: ' t.fullpath]) + +%% Part 7. Explore: Display regions at a lower threshold, alongside meta-analysis regions +% Display the results on slices another way: +% multi_threshold lets us see the blobs with significant voxels at the +% highest (most stringent) threshold, and voxels that are touching +% (contiguous) down to the lowest threshold, in different colors. +% This shows results at p < .001 one-tailed (p < .002 two-tailed), with 10 contiguous voxels. +% Contiguous blobs are more liberal thresholds, down to .02 two-tailed = +% .01 one tailed, are shown in different colors. + +% Set up a display object with two separate brain montages: +o2 = canlab_results_fmridisplay([], 'multirow', 2); + +% Use multi-threshold to threshold and display the t-map: +% Pass the fmridisplay object with montages (o2) in as an argument to use +% that for display. Pass out the object with updated info. +% 'wh_montages', 1:2 tells it to plot only on the first two montages (one sagittal and one axial slice series) + +o2 = multi_threshold(t, o2, 'thresh', [.002 .01 .02], 'sizethresh', [10 1 1], 'wh_montages', 1:2); +o2 = title_montage(o2, 2, 'Behavioral predictor: Reappraisal success (p < .001 one-tailed and showing extent'); + +% Get regions from map from neurosynth.org +% Defined above: metaimg = which('emotion regulation_pAgF_z_FDR_0.01_8_14_2015.nii') + +r = region(metaimg); + +% Add these to the bottom montages +o2 = addblobs(o2, r, 'maxcolor', [1 0 0], 'mincolor', [.7 .2 .5], 'wh_montages', 3:4); + +o2 = title_montage(o2, 4, 'Neurosynth mask: Emotion regulation'); + +% Here, the regions predictive of success are largely different from those +% that respond to reappraisal demand in past literature. + +%% Part 8: Compare to a standard analysis with no scaling or variable transformation + +% Start over and re-load the images: +[image_obj_untouched] = load_image_set('emotionreg', 'noverbose'); + +% attach success scores to image_obj in image_obj.X +% Just mean-center them so we can interpret the group result +image_obj_untouched.X = scale(success, 1); + +% Do a standard OLS regression and view the results: +out = regress(image_obj_untouched, 'nodisplay'); + +% Treshold at q < .05 FDR and display +t_untouched = threshold(out.t, .05, 'fdr'); + +o2 = montage(t_untouched); +o2 = title_montage(o2, 2, 'No scaling: Behavioral predictor: Reappraisal success (q < .05 FDR)'); +o2 = title_montage(o2, 4, 'No scaling: Intercept: Group average activation'); + +% Since we have no results for the predictor, re-threshold at p < .005 +% uncorrected and display so we can see what's there: + +t_untouched = threshold(t_untouched, .005, 'unc'); + +o2 = montage(t_untouched); +o2 = title_montage(o2, 2, 'No scaling: Behavioral predictor: Reappraisal success (p < .005 uncorrected)'); +o2 = title_montage(o2, 4, 'No scaling: Intercept: Group average activation'); +snapnow + +% select just the image for reappraisal success: +t_untouched = select_one_image(t_untouched, 1); + +%% Make plots comparing the analysis with our choices and a naive analysis directly + +o2 = canlab_results_fmridisplay([], 'multirow', 4); + +t = threshold(t, .05, 'fdr'); +o2 = montage(region(t), o2, 'wh_montages', 1:2, 'colormap'); +o2 = title_montage(o2, 2, 'Reappraisal success with informed choices (q < .05 corrected)'); + +t_untouched = threshold(t_untouched, .05, 'fdr'); +o2 = montage(region(t_untouched), o2, 'wh_montages', 3:4, 'colormap'); +o2 = title_montage(o2, 4, 'Reappraisal success with no scaling (q < .05 corrected)'); + +% Here we'll use "multi-threshold" to display blobs at lower thresholds +% contiguous with those at higher thresholds. We'll use uncorrected p < +% .005 so we can see results in both maps: + +o2 = multi_threshold(t, o2, 'thresh', [.002 .01 .02], 'k', [10 1 1], 'wh_montages', 5:6); +o2 = title_montage(o2, 6, 'Reappraisal success with informed choices (p < .001 one-tailed unc, showing extent to .01 one-tailed)'); + +o2 = multi_threshold(t_untouched, o2, 'thresh', [.002 .01 .02], 'k', [10 1 1], 'wh_montages', 7:8); +o2 = title_montage(o2, 8, 'Reappraisal success with no scaling (p < .001 unc, showing extent to .01 one-tailed)'); + +% One-tailed results are the default in SPM and FSL - so this is what you'd +% get with 0.001 in SPM. + + + +%% Part 9: Localize results using a brain atlas +% The table() command is a simple way to label regions +% In addition, we can visualize each region in table to check the accuracy +% of the labels and the extent of the region. + +% Select the Success regressor map +r = region(t); + +% Autolabel regions and print a table +r = table(r); + +% Make a montage showing each significant region +montage(r, 'colormap', 'regioncenters'); + +%% Part 10: Extract and plot data from regions of interest +% Here, we'll explore data extracted from two kinds of regions: + +% 1 - biased regions based on the significant results. This is good for +% visualizing the distribution of the data in a region and checking +% assumptions, but examining the correlation in these regions is +% "circular", as the correlation was used to select voxels in the first +% place. The observed correlations in these regions will be upwardly +% biased as a result. In general, you cannot estimate the effect size +% (correlation or other metric) in a region or voxel after selecting it +% using criteria that are non-independent of the effect of interest. + +% 2 - unbiased regions from the meta-analysis. These are fair. + + + +%% Extract and plot from (biased) regions of interest based on results +% Here, we'll compare the implications of picking BIASED ROIs (which is +% often done in the literature) vs. picking UNBIASED ROIs. +% +% First, let's visualize the correlation scatterplots in the areas we've +% discovered as related to Success + +% Extract data from all regions +% r(i).dat has the averages for each subject across voxels for region i +r = extract_data(r, image_obj); + +% Select only regions with 3+ voxels +% wh = cat(1, r.numVox) < 3; +% r(wh) = []; + +% Set up the scatterplots +nrows = floor(sqrt(length(r))); +ncols = ceil(length(r) ./ nrows); +create_figure('scatterplot_region', nrows, ncols); + +% Make a loop and plot each region +for i = 1:length(r) + + subplot(nrows, ncols, i); + + % Use this line for non-robust correlations: + %plot_correlation_samefig(r(i).dat, datno16.X); + + % Use this line for robust correlations: + plot_correlation_samefig(r(i).dat, image_obj.X, [], 'k', 0, 1); + + set(gca, 'FontSize', 12); + xlabel('Reappraise - Look Neg brain response'); + ylabel('Reappraisal success'); + + % input('Press a key to continue'); + +end + + +%% Extract and plot data from unbiased regions of interest +% Let's visualize the correlation scatterplots in some meta-analysis +% derived ROIs + +% Select the Neurosynth meta-analysis map +r = region(metaimg); + +% Extract data from all regions +r = extract_data(r, image_obj); + +% Select only regions with 20+ voxels +wh = cat(1, r.numVox) < 20; +r(wh) = []; + +r = table(r); + +% Make a montage showing each significant region +montage(r, 'colormap', 'regioncenters'); + +% Set up the scatterplots +nrows = floor(sqrt(length(r))); +ncols = ceil(length(r) ./ nrows); +create_figure('scatterplot_region', nrows, ncols); + +% Make a loop and plot each region +for i = 1:length(r) + + subplot(nrows, ncols, i); + + % Use this line for non-robust correlations: + %plot_correlation_samefig(r(i).dat, datno16.X); + + % Use this line for robust correlations: + plot_correlation_samefig(r(i).dat, image_obj.X, [], 'k', 0, 1); + + set(gca, 'FontSize', 12); + xlabel('Brain'); + ylabel('Success'); + title(' '); + +end + +%% Part 11: (bonus) Multivariate prediction from unbiased ROI averages +% Predict reappraisal success using brain images +% with 5-fold cross-validation + +% Since we are predicting outcome data from brain images, we have to attach +% data in the .Y field of the fmri_data object. +% .Y is the outcome to be explained, in this case success scores +image_obj.Y = image_obj.X; + +% 5-fold cross validated prediction, stratified on outcome + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', 'cv_lassopcrmatlab', 'nfolds', 5); + +% Plot cross-validated predicted outcomes vs. actual + +% Critical fields in stats output structure: +% stats.Y = actual outcomes +% stats.yfit = cross-val predicted outcomes +% pred_outcome_r: Correlation between .yfit and .Y +% weight_obj: Weight map used for prediction (across all subjects). + +% Continuous outcomes: + +create_figure('scatterplot'); +plot(stats.yfit, stats.Y, 'o') +axis tight +refline +xlabel('Predicted reappraisal success'); +ylabel('Observed reappraisal success'); + +% Though many areas show some significant effects, these are not strong +% enough to obtain a meaningful out-of-sample prediction of Success + diff --git a/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m new file mode 100644 index 00000000..1f27baae --- /dev/null +++ b/CanlabCore/Unit_tests/walkthroughs/private/canlab_help_7_multivariate_prediction_basics.m @@ -0,0 +1,391 @@ +%% Multivariate prediction of a continuous outcome + +%% +% This walkthrough uses the |predict()| method for fmri_data objects to +% predict a continuous outcome using cross-validated principal component +% regression (PCR / LASSO-PCR). +% +% The |predict()| method can run multiple algorithms with a range of options. +% The main ones used in the CANlab are SVM for two-choice classification and PCR for +% regression. + +%% About the pain dataset +% ------------------------------------------------------------------------- +% +% The dataset contains data from 33 participants, with brain responses to six levels +% of heat (non-painful and painful). Each image is the average over several +% (4-8) trials of heat delivered at a single stimulus intensity, ranging +% from 44.3 - 49.3 degrees C in one-degree increments. Each image is also +% paired with an average reported pain value for that set of trials, rated +% immmediately after heat experience. +% +% This dataset is interesting for mixed-effects and predictive analyses, as +% it has both within-person and between-person sources of variance. +% +% Aspects of this data appear in these papers: +% Wager, T.D., Atlas, L.T., Lindquist, M.A., Roy, M., Choong-Wan, W., Kross, E. (2013). +% An fMRI-Based Neurologic Signature of Physical Pain. The New England Journal of Medicine. 368:1388-1397 +% (Study 2) +% +% Woo, C. -W., Roy, M., Buhle, J. T. & Wager, T. D. (2015). Distinct brain systems +% mediate the effects of nociceptive input and self-regulation on pain. PLOS Biology. 13(1): +% e1002036. doi:10.1371/journal.pbio.1002036 +% +% Lindquist, Martin A., Anjali Krishnan, Marina López-Solà, Marieke Jepma, Choong-Wan Woo, +% Leonie Koban, Mathieu Roy, et al. 2015. ?Group-Regularized Individual Prediction: +% Theory and Application to Pain.? NeuroImage. +% http://www.sciencedirect.com/science/article/pii/S1053811915009982. +% +%% Load dataset with images and print descriptives +% This dataset is shared on figshare.com, under this link: +% https://figshare.com/s/ca23e5974a310c44ca93 +% +% Here is a direct link to the dataset file with the fmri_data object: +% https://ndownloader.figshare.com/files/12708989 +% +% The key variable is image_obj +% This is an fmri_data object from the CANlab Core Tools repository for neuroimaging data analysis. +% See https://canlab.github.io/ +% +% image_obj.dat contains brain data for each image (average across trials) +% image_obj.Y contains pain ratings (one average rating per image) +% +% image_obj.additional_info.subject_id contains integers coding for which +% Load the data file, downloading from figshare if needed + +fmri_data_file = which('bmrk3_6levels_pain_dataset.mat'); + +if isempty(fmri_data_file) + + % attempt to download + disp('Did not find data locally...downloading data file from figshare.com') + + fmri_data_file = websave('bmrk3_6levels_pain_dataset.mat', 'https://ndownloader.figshare.com/files/12708989'); + +end + +load(fmri_data_file); + +descriptives(image_obj); + +%% Collect some variables we need + +% subject_id is useful for cross-validation +% +% this used as a custom holdout set in fmri_data.predict() below will +% implement leave-one-subject-out cross-validation + +subject_id = image_obj.additional_info.subject_id; + +% ratings: reconstruct a subjects x temperatures matrix +% so we can plot it +% +% The command below does this only because we have exactly 6 conditions +% nested within 33 subjects and no missing data. + +ratings = reshape(image_obj.Y, 6, 33)'; + +% temperatures are 44 - 49 (actually 44.3 - 49.3) in order for each person. + +temperatures = image_obj.additional_info.temperatures; + +%% Plot the ratings +% ------------------------------------------------------------------------- + +create_figure('ratings'); +hold on; plot(ratings', '-', 'Color', [.7 .7 .7], 'LineWidth', .5); +lineplot_columns(ratings, 'color', [.7 .3 .3], 'markerfacecolor', [1 .5 0]); +xlabel('Temperature'); +ylabel('Rating'); +set(gca, 'XTickLabel', [44.3:49.3], 'FontSize', 18); + +%% Prediction options and choices +% There are many choices for algorithm and process/parameter tuning +% Here, we describe a constrained set of some of the most common choices +% +% Base model: Linear, whole brain prediction +% +% Outcome distribution: +% Continuous: Regression (LASSO-PCR) Categorical/2 choice: SVM, logistic regression +% +% Model Options: +% - Feature selection (part of the brain) +% - Feature integration (new features) +% +% Testing options: +% - Input data (what level) +% - Cross-validation +% - What testing metric (classification accuracy? RMSE?) and at what level +% (single-trial, condition averages, one-map-per-subject) +% +% Inference options: +% - Component maps, voxel-wise maps +% - Thresholding (Bootstrapping, Permutation) +% - Global null hypothesis testing / sanity checks (Permutation) + +%% Run the Base model +% +% *relevant functions:* +% predict method (fmri_data) +% predict_test_suite method (fmri_data) +% + +%% +% *Define the holdout set for cross-validation* +% We want to define custom holdout set. If we use subject_id, which is a vector of +% integers with a unique integer per subject, then we are doing +% leave-one-subject-out cross-validation. +% +% Let's build five-fold cross-validation set that leaves out ALL the images +% from a subject together. That way, we are always predicting out-of-sample +% (new individuals). If not, dependence across images from the same +% subjects may invalidate the estimate of predictive accuracy. + +holdout_set = zeros(size(subject_id)); % holdout set membership for each image +n_subjects = length(unique(subject_id)); +C = cvpartition(n_subjects, 'KFold', 5); +for i = 1:5 + teidx = test(C, i); % which subjects to leave out + imgidx = ismember(subject_id, find(teidx)); % all images for these subjects + holdout_set(imgidx) = i; +end + +%% +% *Run the prediction* + +algoname = 'cv_lassopcr'; % cross-validated penalized regression. Predict pain ratings + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +%% Plot cross-validated predicted vs. actual outcomes +% +% Critical fields in stats output structure: +% stats.Y = actual outcomes +% stats.yfit = cross-val predicted outcomes +% pred_outcome_r: Correlation between .yfit and .Y +% weight_obj: Weight map used for prediction (across all subjects). + +% Continuous outcomes: + +create_figure('scatterplots', 1, 2); +plot(stats.yfit, stats.Y, 'o') +axis tight +refline +xlabel('Predicted pain'); +ylabel('Observed pain'); + +% Consider that each subject has their own series of multiple observations. +% Separate subjects into multiple cells, which will be plotted as separate +% individual lines: + +n = max(subject_id); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +subplot(1, 2, 2); +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +% Yfit is a cell array, one cell per subject, with fitted values for that subject + +%% Summarize within-person classification accuracy +% +% Classification can turn a continuous outcome, such as the magnitude of +% activity in an integrated model/pattern, into a binary choice (Type A or +% B, Yes/No, Pain/No pain). Classification accuracy is an easy-to-interpret +% metric, which makes it appealing to report. +% +% But though it seems easy to understand, it is perhaps tricker to +% interpret than we might think. +% +% Classification accuracy is not one metric, but a family of them. The accuracy +% you obtain depends on what type of classification you're doing, and what the +% units of analysis are. For example, is your model based on a single person +% (individualized for that person), or based on other participants' data? +% Are you classifying images corresponding to a single trial, or an average across a group of trials? +% Do you have knowledge about which conditions/images belong to the same person, or not? +% Here, we will: +% (a) use a model based on other people's data (cross-validated across participants) +% (b) classify images summarizing a group of trials, and +% (c) assume we have 2 images belonging to the same person, and we classify +% which is which +% +% On this last point, there are two basic senses of classification, which assume different +% levels of background knowledge: +% - In _single-interval classification_, we do not know which images belong to which participants, +% and we are classifying a single image as Type A or B based on whether +% the response in the model is above or below a threshold. +% - In _forced-choice_ classification, we have two images that we know come +% from the same person, and we know one is Type A and one is Type B...or, +% equivalently, that one is "more A" than the other (e.g., more painful). +% The classifier tries to pick out which is the best answer. +% +% Forced-choice classification is an inherently easier problem, with higher +% resulting accuracy, because (a) each participant serves as their own control +% (many sources of noise are canceled out), and (b) you assume more +% background knowledge, conferring an advantage. +% +% In this case, we want to summarize each subject with two images: One high +% pain, and one low pain. We'll use stimulus intensity as the variable to +% classify here, and ask if we can predict which intensity is highest given +% pairs of images one degree apart for each participant. + +% If the stim intensity increases by 1 degree, how often do we correctly +% predict an increase? +diffs = cellfun(@diff, Yfit, 'UniformOutput', false); % successive increases in predicted values with increases in stim intensity + +% Subjects (rows) x comparisons (cols), each comparison 1 degree apart +diffs = cat(2, diffs{:})'; + +acc_per_subject = sum(diffs > 0, 2) ./ size(diffs, 2); +acc_per_comparison = (sum(diffs > 0) ./ size(diffs, 1))'; +cohens_d_per_comparison = (mean(diffs) ./ std(diffs))'; + +fprintf('Mean acc is : %3.2f%% across all successive 1 degree increments\n', 100*mean(acc_per_subject)); + +% Create and display a table: + +comparison = {'45-44' '46-45' '47-46' '48-47' '49-48'}'; +results_table = table(comparison, acc_per_comparison, cohens_d_per_comparison); +disp(results_table) + +% 45-44 degrees, 46-45, 47-46, etc. + +%% Visualize the classifier weight map +% +% The weight map is stored in |stats.weight_obj|, as a statistic_image +% class object. If we have requested bootstrapping, the image will be +% thresholded (0.05 uncorrected by default). Otherwise, it will be unthresholded, +% and we'll see weights everywhere within the analysis mask. + +w = stats.weight_obj; % This is an image object that we can view, etc. + +orthviews(w) + +create_figure('montage'); +o2 = montage(w); +o2 = title_montage(o2, 5, 'Predictive model weights'); + + +%% Bootstrap weights: Get most reliable weights and p-values for voxels +% +% Here is an exercise: +% Re-run the predictive model, adding a 'bootstrap' flag. +% For a final analysis, at least 5K bootstrap samples is a good idea. +% For now, try it with just 1,000 bootstrap samples (and be prepared to +% wait a bit). What does the resulting weight map look like? +% +% (This is not run in the example; see |help fmri_data.predict| for how to +% add the bootstrap option.) + +%% Try normalizing/scaling data +% +% One of the biggest sources of noise can be whole-image (or +% large-spatial-scale shifts) in image values across participants. These can +% be apparent even in contrast or "subtraction" images where various +% sources of noise and artifacts are supposed to have been "subtracted +% out". This may result from different chance correlations between task +% regressors and physiological noise or artifacts/outliers across +% individuals. +% +% Rescaling the data can remove whole-brain signal. This should be used +% with caution, as you're also removing signal from the images, and the +% resulting weight maps should be interpreted as effects *relative* to the mean +% signal across the image. i.e., a task may be associated with whole-brain +% increases, and *relative* decreases in a local area that appear only +% after the global increase has been accounted for. +% +% Another reason to rescale images is to ask whether the *pattern* of +% activity across the brain (or within a local area) is predictive of task +% state, after removing the overall intensity of activation. This relates +% to pattern information. +% +% Here, we'll rescale the images and ask if this improves our +% cross-validated predictions or not. If it does, the global signal we've +% removed is likely more noise than signal. If not, the global signal may +% contain information about pain. + +% z-score every image, removing the mean and dividing by the std. +% Remove some scaling noise, but also some signal... + image_obj = rescale(image_obj, 'zscoreimages'); + +% re-run prediction and plot + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +create_figure('scatterplots'); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +fprintf('Overall prediction-outcome correlation is %3.2f\n', stats.pred_outcome_r) + +figure; o2 = montage(stats.weight_obj); +o2 = title_montage(o2, 5, 'Predictive model weights'); + +%% Try selecting an a priori 'network of interest' + +% re-load +load('bmrk3_6levels_pain_dataset.mat', 'image_obj') + +% load an a priori pain-related mask and apply it + +mask = fmri_data(which('pain_2s_z_val_FDR_05.img.gz')); % in Neuroimaging_Pattern_Masks repository + +image_obj = apply_mask(image_obj, mask); + +% show mean data with mask +m = mean(image_obj); +figure; o2 = montage(m); +o2 = title_montage(o2, 5, 'Mean image data, masked'); + + +%% re-run prediction and plot + +[cverr, stats, optout] = predict(image_obj, 'algorithm_name', algoname, 'nfolds', holdout_set); + +create_figure('scatterplots'); + +for i = 1:n + YY{i} = stats.Y(subject_id == i); + Yfit{i} = stats.yfit(subject_id == i); +end + +line_plot_multisubject(Yfit, YY, 'nofigure'); +xlabel('Predicted pain'); +ylabel('Observed pain'); +axis tight + +fprintf('Overall prediction-outcome correlation is %3.2f\n', stats.pred_outcome_r) + +figure; o2 = montage(stats.weight_obj); +o2 = title_montage(o2, 5, 'Predictive model weights'); + +%% Other ideas +% There are many ways to build from this basic analysis. + +% Data +% ------------ +% Examine/clean outliers +% Transformations of data: Scaling/centering to remove global mean, +% ventricle mean +% Identify unaccounted sources of variation/subgroups (experimenter sex, +% etc.) + +% Features +% ------------ +% A priori regions / meta-analysis +% Feature selection in cross-val loop: Univariate, recursive feature elim +