# Representing multi-stage neuroimaging models

The statistical analysis of neuroimaging data typically occurs across distinct stages of analysis, with parameter estimates from lower levels of analysis propagating to higher levels for subsequent analysis.

For example, in fMRI it is common to first fit a design matrix to run-level time series followed by a fixed-effects model to combine estimates at the subject-level. Finally, a dataset-level (or "group level") random-effects one-sample t-test can be performed to estimate population level effects. At each level of the analysis, we need to know which image inputs correspond to which design matrix, and more how to keep track of and combine outputs from the previous level at the current level of analysis.

*BIDS Stats Models* proposes a general machine-readable document
to describe multi-stage neuroimaging analyses in a precise, yet flexible manner.
We accomplish this by defining a *graph* composed of **Nodes** representing each level of the analysis
and **Edges** which define the flow of data from one Node to another.
Within each {py:class}`~bsmschema.models.Node`

we specify
a {py:class}`~bsmschema.models.Model`

to estimate,
and at least one {py:class}`~bsmschema.models.Contrast`

to define the computed outputs of each `Node`

.
Within each node we also specify how to group the incoming inputs
into analysis units using the {py:attr}`~bsmschema.models.Node.GroupBy`

directive.

## A simple example

In a Simon task,
participants were scanned for 2 runs and asked to indicate whether a diamond
that was presented to the left or right of a central fixation cross was green or red.
There were two conditions: color-spatial *congruent* and *incongruent* trials.

A basic analysis of this task is to determine
which regions showed *greater activity for incongruent versus congruent trials*,
across all participants.

We can perform this analysis by first estimating a **run-level** timeseries model
for "Incongruent" and "Congruent" trials--separately for each individual run.
We then compute a contrast comparing *Incongruent > Congruent* (IvC) trials.

Next, we pass the resulting statistical maps for the contrast to a **subject-level** estimator,
which computes the average *IvC* effect for each subject separately.

Finally, we pass the resulting estimates to a **dataset-level** estimator,
which conduts a one-sample t-test across all of the subject estimates for the *IvC* contrast.

Let's visualize this model for 3 participants:

We can formally represent this analysis as **BIDS Stats Model**:

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON

```
BIDS Stats Models *must* have a {py:attr}`~bsmschema.models.BIDSStatsModel.Name`
and {py:attr}`~bsmschema.models.BIDSStatsModel.BIDSModelVersion` defined,
and optionally can restrict input images with {py:attr}`~bsmschema.models.BIDSStatsModel.Input`.
```{note}
For this example, we have limited the model to three subjects using the `Input` key.
```

*BIDS Stats Models* defines this multi-stage analysis as a Graph,
with each level of analysis defined as a separate `Node`

object.
Let's step through each `Node`

separately.

### Run-level Model

First, we define a `Node`

for the run level analysis.

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :lines: 6-19

```
Note that the {py:attr}`~bsmschema.models.Node.Level` key is necessary
for determining which input images are available to the estimator.
At the `Run` level, there are many sources of possible variables,
most notably `_events.tsv` files which define the timing of task-related events.
Next we define a {py:attr}`~bsmschema.models.Model` for this node.
```{literalinclude} examples/model-walkthrough_smdl.json
:language: JSON
:start-at: '"Model": {"X": [1,'
:lines: 1
```

The {py:attr}`~bsmschema.models.Model.X`

parameter defines the variables in the design matrix.
Here, we are modeling the `incongruent`

and `congruent`

trial types,
in addition to an intercept (identified by the special key: `1`

).

Next, we specify an *Incongruent-Congruent (IvC)* contrast
using the {py:attr}`~bsmschema.models.Contrast`

key:

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :lines: 11-18

```
If you have used other fMRI modeling tools this should be familiar.
We have specified a t-test contrast with the weights `[1, -1]` for the conditions:
`["incongruent", "congruent"]` and given this contrast the name `IvC`.
```{attention}
`Contrasts` **define the outputs** that will be available to the next `Node`.
Since we only modeled a single contrast (`IvC`),
the next `Node` will not have access to estimates
for main effects for the `congruent` or `incongruent` conditions,
unless we explicitly define a `Contrast` for each.
```

#### How to *group* analysis inputs?

An underappreciated factor in multi-stage models is
the grouping of image inputs into analysis units.
For example, here we want to estimate a timeseries model for each `Run`

separately,
rather that concatenating all runs for each subject into one large model.

We must explicitly define this grouping structure
using the {py:attr}`~bsmschema.models.Node.GroupBy`

key for every node.
To fit a separate time series model for each individual run image, we specify:

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :start-at: ' "GroupBy": ["run", "subject"]' :lines: 1

```
Here, `GroupBy` states that for every unique combination of `run` and `subject`,
we will fit a separate model. This results in a single input image per model.
If you are familiar with tabular data such as R `DataFrames`, or `pandas`,
the `GroupBy` operation should be familiar.
For instance, given three subjects with two runs each, we can define 6 rows in a table (3x2):
```{code-cell} python3
---
tags: ["remove_input"]
---
from IPython.display import display
import pandas as pd
pd.set_option('display.max_colwidth', None)
subjects = ["01", "02", "03"]
runs = [1, 2]
# For later use
contrasts = ["IvC"]
stats = ["effect", "variance"]
def display_groups(df, groups):
for group in df.groupby(groups, as_index=False):
display(group[1].style.hide(axis="index"))
inputs = pd.DataFrame.from_records(
[{
"subject": subject,
"run": run,
"image": f"sub-{subject}_task-simon_run-{run}_bold.nii.gz",
} for subject in subjects for run in runs]
)
display(inputs.style.hide(axis="index"))
```

If we `GroupBy`

*subject*, there would be three groups of images--one for each subject:

## ```{code-cell} python3

## tags: ["remove_input"]

display_groups(inputs, "subject")

```
If we `GroupBy` *run*, all images with the same *run* ID would be grouped together,
resulting in two groups, one for each distinct group ID:
```{code-cell} python3
---
tags: ["remove_input"]
---
display_groups(inputs, "run")
```

However, since we want to model each `BOLD`

image separately,
we must `GroupBy`

**both subject and run**,
resulting in six groups with a single image each.

## ```{code-cell} python3

## tags: ["remove_input"]

display_groups(inputs, ["run", "subject"])

```
### Subject level Node
At this point, we have defined a `Model` that will be fit separate to each grouping--
in this case a separate time-series model for each `run`.
Next, we want to define a `subject` level node
to pool together estimates from this `Node` for each `subject`
using a fixed-effects model.
```{literalinclude} examples/model-walkthrough_smdl.json
:language: JSON
:lines: 20-26
```

```
By default, `Nodes` are linked in sequential order,
with all the `Contrast` outputs from a `Node` available to the subsequent `Node`.
```

#### From Run Outputs to Subject Inputs

We need to use `GroupBy`

to define how to group the outputs from the `Run`

node
as inputs to the `Subject`

level:

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :lines: 23

```
Here we are specifying that all images belonging to a single `subject`
and from a single `contrast` should be grouped into a unit of analysis.
Note that with 3 subjects and 2 runs, we will have 6 groups of output images from the `Run` node.
Given two types of images (`variance` and `effect`),
this results in 12 images that would be grouped as follows:
```{code-cell} python3
---
tags: ["remove_input"]
---
outputs = pd.DataFrame.from_records(
[{
"subject": subject,
"run": run,
"contrast": contrast,
"image": f"sub-{subject}_task-simon_run-{run}_"
f"contrast-{contrast}_stat-{stat}_statmap.nii",
}
for subject in subjects for run in runs
for contrast in contrasts for stat in stats]
)
display_groups(outputs, ["subject", "contrast"])
```

```
Although there is only one `contrast`, we include `contrast` as a grouping variable to be explicit.
```

#### Subject Model

We can now specify the `Subject`

level `Model`

.

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :start-at: '"Model": {"X": [1], "Type": "meta"}' :lines: 1-2

```
Since our intent is to estimate the _mean effect_ for each subject,
we only need an intercept in our model.
We specify the `"Type"` to be `Meta`, which is a special type to identify fixed-effects models.
```{note}
`1` is a special variable used to represent the intercept.
```

Remember that we must specify `Contrasts`

in order to produce outputs for the next `Node`

.
`DummyContrasts`

is a convenience function which will create contrasts
with the weights `[1]`

for all modeled inputs.
Since we are not comparing anything at the subject-level
and simply want to pass forward the estimates generated by the fixed-effects model,
this is useful and saves us from specifying a more verbose (but identical) `Contrast`

.

### Dataset level Node

We are ready to perform a one-sample t-test to estimate population-level effects
for the *IvC* `Contrast`

.
We refer to this level as the `Dataset`

level.

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :start-at: '"Level": "Dataset",' :lines: 1-3

```
Here we only need to `GroupBy: ['contrast']` as we want to compute a separate estimate
for each contrast, but want to include all subjects in the same analysis.
Since we only have one `contrast`, all the incoming subject-level images will be grouped together:
```{code-cell} python3
---
tags: ["remove_input"]
---
outputs = pd.DataFrame.from_records(
[{
"subject": subject,
"contrast": contrast,
"image": f"sub-{subject}_task-simon_"
f"contrast-{contrast}_stat-{stat}_statmap.nii",
}
for subject in subjects
for contrast in contrasts for stat in stats]
)
display_groups(outputs, ["contrast"])
```

As before, we can specify an intercept-only model, but of type `glm`

since we want to perform a random-effects analysis.
We can again use `DummyContrasts`

to specify a simple one-sample t-test contrast
on the incoming `IvC`

subject-level contrasts.

```{literalinclude} examples/model-walkthrough_smdl.json :language: JSON :lines: 31-32

```
The outputs of this node collapse across subjects, leaving a single effect/variance pair:
```{code-cell} python3
---
tags: ["remove_input"]
---
outputs = pd.DataFrame.from_records(
[{
"contrast": contrast,
"image": f"task-simon_contrast-{contrast}_stat-{stat}_statmap.nii",
}
for contrast in contrasts for stat in stats]
)
display_groups(outputs, ["contrast"])
```

## Ready to run ðŸš€

At this point, we have a fully specified model three-stage fMRI model.

Our model will compute a run level incongruent-congruent contrast,
pass forward the estimates to a fixed-effects model to pool subject estimates,
and compute a dataset-level random-effects model and one-sample t-test
to estimate population effects for the `IvC`

contrast.

We can now pair this *BIDS Stats Models* specification with a pre-processed derivative
from the original raw dataset and hand these to a tool that supports *BIDS Stats Models*
for **fully automated execution**.

## Next up

- Read the next section to dive deeper into advanced usage of
*BIDS Stats Models*to enable more complex models. - Check out tools like FitLins
which implement
*BIDS Stats Models*execution to learn how to run a*BIDS Stats Model* - Take a look at more example models in the
**.**