Group Analysis

exploreFNIRS is Layer 2 of processFNIRS2: it takes a cell array of processed fNIRS structs (the output of processFNIRS2, typically cut into event-locked segments — see Block-Averaging-and-Epoching) and turns them into grand averages, group plots, and statistics. This page covers the exploreFNIRS.core.Experiment class end to end: building it, aggregating, behavioral/ROI/aux covariates, LME statistics, GLM workflows, and Growth Curve Analysis.

The Experiment class

% allData = {seg1, seg2, seg3, ...};  % cell array of processed/segmented structs
ex = exploreFNIRS.core.Experiment(allData);
ex = exploreFNIRS.core.Experiment(allData, ...
    'Hierarchy', {'SubjectID', 'Session', 'Condition', 'Trial', 'Block'});  % default hierarchy

Each element's .info struct (SubjectID, Group, Condition, ...) becomes a column in ex.dataTable, built internally by exploreFNIRS.dataset.buildSegmentInfoTable. A quick way to try this without your own data:

[ex, allData] = pf2.import.sampleData.group();   % ready-made, grouped, aggregated Experiment

Settings

All preprocessing/aggregation options live on ex.settings, set before calling aggregate():

Setting	Default	Description
`baseline`	`[-5, 0]`	Baseline window `[start, end]` in seconds
`taskStart`	`0`	Task onset time for bin alignment
`taskEnd`	`Inf`	Task end time (`Inf` = use full segment)
`resampleRate`	`0.5`	Seconds per bin for temporal plots (`0` = no resampling)
`barBinSize`	`0`	Seconds per bin for bar/LME data (`0` = one bar over the full task window)
`useBaseline`	`true`	Apply baseline correction
`avgMode`	`'hierarchy'`	Within-subject averaging: `'hierarchy'`, `'flat'`, `'none'`
`rawMethod` / `oxyMethod`	`''`	Reprocess segments with a named pipeline before aggregating (`''` = use data as-is)
`statWindow`	`[]`	`[start, end]` restricting LME/bar stats to a time subset (`[]` = full range)
`viewPad`	`[5, 5]`	Plot-only padding around baseline-start/task-end (does not affect bar values)
`timeModel`	`''`	Default `TimeModel` injected into LME calls (`''` = fitLME default) — see Growth Curve Analysis
`polyOrder`	`2`	Default polynomial order for GCA

ex.settings.baseline = [-5, 0];
ex.settings.resampleRate = 1;      % 1 Hz for temporal plots
ex.settings.barBinSize = 10;       % 10 s bins for bar charts / LME
ex.settings.useBaseline = true;

select, groupby, aggregate

ex.select('Group', 'Control');                          % narrow by metadata (AND logic)
ex.select('Condition', {'Easy', 'Hard'});                % match any of a list
ex.reset();                                               % clear selection

ex.groupby({'Group', 'Condition'});                       % one group per unique combination
ex.aggregate();                                            % run the full pipeline (equivalent to ex.aggregate('hierarchy'))

aggregate() runs in two stages: a cached preprocessing pass (baseline extraction via pf2.data.split, temporal/bar resampling via pf2.data.resample) and a grand-averaging pass (grandAvgFNIRS) that always re-runs. The preprocessing cache is keyed on baseline/resampleRate/barBinSize/taskStart/useBaseline; changing only avgMode and re-calling aggregate('flat') reuses the cache and just re-averages.

After aggregation, each group struct (ex.groups) carries gbyGrand (temporal resolution), gbyGrandBarFlat (bar/LME resolution, flat hierarchy), and gbyFNIRS_pp (preprocessed segments). ex.summary() prints a quick overview of the current groups.

Averaging modes

ex.settings.avgMode (or ex.aggregate(mode)) controls how repeated observations per subject are combined before the group mean:

Mode	Behavior	Use when
`'hierarchy'` (default)	Averages bottom-up through the hierarchy (Trial → Condition → Session → Subject) before the group mean; every subject contributes equally	Standard group analysis / required for valid between-group inference
`'flat'`	Pools all of a subject's observations (across conditions/trials) into one value	You want one "overall" value per subject, ignoring condition structure
`'none'`	No within-subject averaging — every segment is an independent observation	Exploratory visualization only; not valid for inference without accounting for repeated measures

ex.aggregate('hierarchy');   % default
ex.aggregate('flat');        % cache hit if preprocessing settings unchanged — only re-averages
ex.aggregate('none');

The risk with 'none' is pseudoreplication: a subject with 20 trials gets 4x the weight of one with 5 in the group mean, and error bands mix within- and between-subject variability. LME models (plotLME/statsFitLME) sidestep this at the modeling level via the (1|SubjectID) random intercept regardless of avgMode, but bar/scatter plots and simple group means do not. See example_averaging_modes.m for a side-by-side comparison.

Behavioral variables

Metadata-only variables (reaction time, accuracy, age, questionnaire scores, ...) can be plotted and modeled directly from .info — no aggregate() required:

ex.groupby({'Condition'});
fig = ex.plotInfoBar('reactionTime', 'YLabel', 'Reaction Time (ms)');
fig = ex.plotInfoScatter('reactionTime', 'accuracy', 'FitLine', true, 'ErrorBand', true);

% Single LME model with the behavioral variable as the response
% (unlike statsFitLME, which iterates over fNIRS channels, this fits one model)
rtResults = ex.statsInfoLME('reactionTime');
fprintf('%s\n', rtResults.formula);      % e.g. 'reactionTime ~ Condition + (1|SubjectID)'
disp(rtResults.anova_pval);

plotInfoBar/plotInfoScatter respect avgMode just like fNIRS plots (default 'hierarchy': one point per subject).

ROI-level analysis

Define regions of interest over channels, then run the same plotting/stats API at the ROI level instead of per-channel:

nCh = size(allData{1}.HbO, 2);
roiDef   = {1:6; 7:12; 13:18};                         % optode index lists (column cell array)
roiNames = {'Left_PFC', 'Center_PFC', 'Right_PFC'};
for i = 1:numel(allData)
    allData{i} = pf2.probe.roi.defineROI(allData{i}, roiDef, roiNames);
    allData{i} = pf2_build_nanmean_ROI(allData{i});     % adds ROI-averaged biomarker fields
end

ex = exploreFNIRS.core.Experiment(allData);
ex.groupby({'Condition'});
ex.aggregate();

fig = ex.plotTemporal('Biomarkers', {'HbO'}, 'ROIs', 'all');
fig = ex.plotBar('Biomarker', 'HbO', 'ROIs', 'all', 'ShowIndividual', true);
[fig, results] = ex.plotLME('Biomarkers', {'HbO'}, 'DataType', 'ROI');
roiResults = ex.statsROILME('Biomarkers', {'HbO'});      % stats only, convenience wrapper
[fig, ~] = ex.plotTopoROILME('Biomarkers', {'HbO'});      % ROI-level F-stats on the 3D brain

T = ex.toLongTable({'HbO'}, [], [], 'IncludeROI', true);  % export with ROI column

ROIs on plotTemporal/plotBar/plotLME is mutually exclusive with Channels and automatically switches DataType to 'ROI'.

Auxiliary signals at the group level

Aux fields (heart rate, accelerometer, ...) flow through the same preprocessing + grand-average pipeline as fNIRS biomarkers once aggregate() has run (see Auxiliary-Signals for signal-level details):

ex.auxFields();                                    % list available Aux fields
fig = ex.plotAux('accelerometer', 'Layout', 'grid');
fig = ex.plotAux('heartRate', 'ErrorType', 'SEM');
fig = ex.plotAuxBar('heartRate', 'TimeWindow', [5, 20], 'ShowIndividual', true);

hrResults = ex.statsAuxLME('heartRate');            % channel-wise LME on an Aux field

Aux signals can also enter fNIRS/behavioral LME models as covariates:

% Promote an aux_* signal to a fixed-effect predictor in statsAutoLME
results = ex.statsAutoLME('Biomarkers', {'HbO'}, 'AuxCovariates', {'heartRate'});

% Force a per-trial behavioral covariate (e.g. RT) into every candidate model
results = ex.statsAutoLME('Covariates', {'reactionTime'});

Continuous covariates enter uncentered — consider mean-centering (e.g. zscore) beforehand so intercepts stay interpretable. A forced covariate that is constant within every subject is confounded with the (1|SubjectID) random intercept and triggers a warning.

Statistics: LME modeling

Two entry points fit per-channel linear mixed-effects models using the current groupby() variables as fixed effects and (1|SubjectID) as the default random effect:

Method	Model	Use when
`ex.statsFitLME(...)`	Fixed formula: additive main effects by default, or `'AllInteractions', true` for the full interaction model	You know the model you want
`ex.statsAutoLME(...)`	Forward stepwise selection per channel (AIC/BIC)	You want per-channel automatic model discovery; `'Criterion'`, `'DeltaThreshold'`, `'Covariates'` control the search

ex.plotLME(...) wraps statsFitLME-style fitting with a bar-chart visualization of F-statistics per channel/term in one call; use statsFitLME/statsAutoLME directly for statistics only (no plot).

ex.groupby({'Condition'});
ex.aggregate();

% Visual: bar chart of F-stats per channel + fitted models
[fig, results] = ex.plotLME('Biomarkers', {'HbO'}, 'Channels', [1, 2, 3]);
fprintf('%s\n', results.formula);
disp(results.anova_pval);

% Pure stats path (no plot)
lmeResults = ex.statsFitLME('Biomarkers', {'HbO'}, 'Channels', 1:3);
% or: lmeResults = ex.statsAutoLME('Biomarkers', {'HbO'}, 'Channels', 1:3);

results (from either plotLME or statsFitLME/statsAutoLME) contains .models, .anova, .anova_pval / .anova_Fstat ([channels x terms] tables), .contrasts, .AIC, .formula, and .mergedTable.

Contrasts, FDR, and publication tables

contrastResults = ex.statsRunContrasts(lmeResults);          % post-hoc contrasts, FDR-corrected across channels
% ex.statsRunContrasts(lmeResults, 'FDRThreshold', 0.05, 'FDRMethod', 'bh');  % 'bh' (default) or 'twostep'

anovaT     = ex.statsSummarize(lmeResults, 'Type', 'anova', 'Format', 'apa');       % or 'console'/'latex'/'table'
contrastT  = ex.statsSummarize(contrastResults, 'Type', 'contrasts', 'Format', 'console');

statsSummarize handles Type = 'anova', 'contrasts', 'coefficients', 'fit', 'correlations', each in Format = 'table' (default), 'console', 'latex', or 'apa'.

For small-N designs, statsPermTest (sign-flip permutation, paired 2-condition) and statsEffectSize (Hedges' g / Cohen's d with bootstrap CIs) are non-parametric/effect-size alternatives to LME:

permResults = ex.statsPermTest('Biomarkers', {'HbO'}, 'NumPerm', 1000, 'Statistic', 'tstat');
esResults   = ex.statsEffectSize('Biomarkers', {'HbO'}, 'Method', 'hedges_g', 'NumBoot', 2000);

Bridging group statistics to brain projection

plotLME/statsFitLME/statsAutoLME compute the statistics; pf2.probe.project.* only renders a stat vector you already have. Two ways to combine them:

A — manual bridge (when you need the raw p-values/F-stats to threshold, FDR-correct, or export before rendering):

[ex, allData] = pf2.import.sampleData.group();
seg  = allData{1};
nOpt = size(seg.HbO, 2);          % full optode count, including short-separation

[~, results] = ex.plotLME('Biomarkers', {'HbO'}, 'Visible', 'off');

% results.channels are absolute optode indices for the modeled (non-SS) channels
pCondition = results.anova_pval.Condition;      % [nModeled x 1]
pFull = nan(1, nOpt);
pFull(results.channels) = pCondition;           % unmodeled optodes stay NaN -> render transparent

pf2.probe.project.pvalues(pFull, seg, 'pThreshold', 0.05, 'includeSS', true, ...
    'titleString', 'Condition effect (p)');

Fcondition = results.anova_Fstat.Condition;
Ffull = nan(1, nOpt);
Ffull(results.channels) = Fcondition;
pf2.probe.project.fstats(Ffull, seg, 'includeSS', true, 'titleString', 'Condition effect (F)');

B — all-in-one (plotTopoLME computes and projects in a single call; use it when you don't need the raw vectors):

ex.plotTopoLME('Biomarkers', {'HbO'});                                   % 3D F-stat topo per term
ex.plotTopoLME('Biomarkers', {'HbO', 'HbR'});                            % multi-biomarker rows
ex.plotTopoLME('Biomarkers', {'HbO'}, 'SigType', 'q', 'SigThreshold', 0.05);  % FDR-corrected
ex.plotTopoLME('Biomarkers', {'HbO'}, 'PlotMetric', 'p');                % -log10(p) topo
ex.plotTopoROILME('Biomarkers', {'HbO'});                                % ROI-level equivalent

3D renders need the 'savePath' option (or capture imgOut/imwrite) to save headlessly — see Visualization.

GLM workflows

The epoch-averaging path (defineBlocks → extractBlocks → Experiment, see Block-Averaging-and-Epoching) cuts continuous recordings around events and averages. GLMExperiment instead keeps the continuous recording intact and fits a linear model with HRF-convolved regressors — handling overlapping events, explicit drift/nuisance modeling, and producing per-condition beta weights as the group-analysis input.

GLMExperiment (automated)

[subjects, blockDefs] = pf2.import.sampleData.experiment('blocks');   % continuous recordings + block defs
[rawMethod, oxyMethod] = pf2_base.examples.addDemoPipelines();

gx = exploreFNIRS.core.GLMExperiment(subjects, blockDefs);
% gx = exploreFNIRS.core.GLMExperiment(subjects)   % uses subjects{i}.blocks if already embedded

gx.settings.rawMethod = rawMethod;
gx.settings.oxyMethod = oxyMethod;
gx.glm.conditions  = {'Easy', 'Hard'};    % which regressors to include (empty = auto-detect)
gx.glm.auxFields   = {'heartRate'};       % also fit a GLM on an Aux signal
% gx.glm.driftOrder, gx.glm.biomarkers, gx.glm.fitMethod ('OLS' default or 'AR-IRLS')

gx.fit();   % process -> design matrix -> per-subject GLM -> pseudo-segments -> behavioral aggregation

r1 = gx.getSubjectResult(1);
fprintf('Sub01 mean R2 (HbO): %.3f\n', mean(r1.results.HbO.R2));
fig = gx.plotDesignMatrix(1);

GLMExperiment extends Experiment, so after fit() every plotting/stats/export method above works unchanged on the beta pseudo-segments (useBaseline=false, resampleRate=0, barBinSize=0 are set automatically):

gx.groupby({'Condition'});
gx.aggregate();
fig = gx.plotBar('Biomarker', 'HbO', 'Channels', 1:4, 'ShowIndividual', true);
[fig, lmeResults] = gx.plotLME('Biomarkers', {'HbO'}, 'Channels', 1:4);
T = gx.betaTable('Channels', 1:4, 'IncludeStats', true);   % direct beta export (bypasses Experiment)

Manual design-matrix path

For full control over design construction, per-subject fitting, and first-level contrasts:

events = pf2.data.blocksToEvents(blockDefs{1});
[X, regressorNames] = pf2_base.fnirs.buildDesignMatrix(d.time, d.fs, events, 'DriftOrder', 3);
% 'DriftType' ('legendre'|'dct'), 'HRF', 'IncludeDerivative', 'IncludeDispersion', 'ShortChannels'

glmResult = pf2_base.fnirs.fitGLM(d.HbO, X, regressorNames);

% First-level contrast: Hard > Easy
easyIdx = find(strcmp(regressorNames, 'Easy'));
hardIdx = find(strcmp(regressorNames, 'Hard'));
C = zeros(1, numel(regressorNames)); C(hardIdx) = 1; C(easyIdx) = -1;
contrastResult = pf2_base.fnirs.fitGLM(d.HbO, X, regressorNames, ...
    'Contrasts', C, 'ContrastNames', {'Hard > Easy'});

% Package betas into Experiment-compatible pseudo-segments
segs = pf2.data.betasToSegments(glmResult, d, 'Conditions', {'Easy', 'Hard'}, ...
    'BiomarkerResults', glmResult);
ex = exploreFNIRS.core.Experiment([allSegments, segs]);

See example_glm_advanced.m for the complete multi-subject loop.

Growth Curve Analysis

Growth Curve Analysis (GCA) models the shape of the response over time — rather than collapsing each block to a single mean — using orthogonal polynomial time terms inside the LME. It answers "do groups/conditions differ in how the response unfolds?" as opposed to a GLM ("is there a response at all?") or a plain mean ("do groups differ in amplitude?").

Set barBinSize to divide each block into several time bins (rule of thumb: 4-15 bins; fewer than 3 auto-clamps the polynomial order), then let statsFitLME/plotLME add polynomial time terms automatically:

ex.settings.barBinSize = 10;   % e.g. 60s block -> 6 time points

% Default TimeModel='polynomial', PolynomialOrder=2 produces a formula like:
%   HbO ~ Condition + ot1 + ot2 + Condition:ot1 + Condition:ot2 + (1+ot1|SubjectID)
results = ex.statsFitLME('Biomarkers', {'HbO'}, 'Channels', 1:5);
disp(results.anova_pval);   % Condition:ot1 -> slope differs; Condition:ot2 -> curvature differs

TimeModel options (also settable persistently via ex.settings.timeModel/ex.settings.polyOrder):

`TimeModel`	Behavior
`'polynomial'` (default)	Orthogonal polynomial time terms (`ot1`, `ot2`, ...) — best for smooth response shapes
`'discrete'`	Each time bin as an independent categorical level
`'continuous'`	Single centered numeric time predictor (linear trend only)
`'none'`	Drops time from the model (equivalent to `barBinSize = 0`)

r = ex.statsFitLME('Biomarkers', {'HbO'}, 'Channels', 1, 'TimeModel', 'discrete');
r = ex.statsFitLME('Biomarkers', {'HbO'}, 'Channels', 1, 'PolynomialOrder', 3);
r = ex.statsFitLME('Biomarkers', {'HbO'}, 'Channels', 1, 'StatWindow', [10, 50]);  % restrict the time range

Reference: Mirman, D. (2017). Growth Curve Analysis and Visualization Using R. Chapman and Hall/CRC. DOI: 10.1201/9781315373218

Group plotting entry points

plotTemporal, plotBar, plotScatter, and plotHeatmap are the core group visualizations, each usable directly or through the dimension-mapping ex.plot proxy (ex.plot.bar('X', 'Condition', 'Color', 'Group', ...)). See Visualization for the full option set (error bands, color schemes, topo maps, saving).

fig = ex.plotTemporal('Biomarkers', {'HbO', 'HbR'}, 'Channels', 1:5, 'ErrorType', 'SEM');
fig = ex.plotBar('Biomarker', 'HbO', 'Channels', 1:5, 'PlotBy', 'Condition', 'ShowIndividual', true);
[fig, stats] = ex.plotScatter('reactionTime', 'Biomarkers', {'HbO'}, 'Channels', 1, 'FitLine', true);
fig = ex.plotHeatmap('Biomarker', 'HbO', 'Channels', 1:8);

Export

longT = ex.toLongTable({'HbO', 'HbR'}, 1:5);          % long format for R/lme4, tidyverse
wideT = ex.toWideTable({'HbO'}, 1:5);                 % wide format for SPSS/Excel
T     = ex.infoTable();                                % full metadata table
T     = ex.demographicsTable('Variables', {'Age', 'Group'}, 'GroupBy', 'Group');  % Table 1
writetable(longT, 'results.csv');

Navigation Menu

Repository navigation

Group Analysis

Group Analysis

The Experiment class

Settings

select, groupby, aggregate

Averaging modes

Behavioral variables

ROI-level analysis

Auxiliary signals at the group level

Statistics: LME modeling

Contrasts, FDR, and publication tables

Bridging group statistics to brain projection

GLM workflows

GLMExperiment (automated)

Manual design-matrix path

Growth Curve Analysis

Group plotting entry points

Export

See also

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally