Repository: bids-standard/pybids

Issue:
modeling: CLI pre-run report
From discussion with @effigies @jmumford :
 
Contents of the report:
 - model in JSON
 - graph representation of the model
 - DM pre/post convolution 
      - raw numbers (combined df with n/as) + plot (scale for viz purposes0
 - VIF
 - pair wise correlation--- pearson--- post convolution

warnings:
- VIF > 5 > 10
- n/a becomes a zero
- n/a prior to scale transformation?

Add to CLI ability to run this report

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/bids-standard/pybids/pull/918?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard) Report
Base: **86.24**% // Head: **63.24**% // Decreases project coverage by **`-23.00%`** :warning:
> Coverage data is based on head [(`e00afcc`)](https://codecov.io/gh/bids-standard/pybids/pull/918?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard) compared to base [(`de95cb5`)](https://codecov.io/gh/bids-standard/pybids/commit/de95cb580d838c7a08bf436ad2c8f94cc0cfd058?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard).
> Patch coverage: 5.08% of modified lines in pull request are covered.

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           master     #918       +/-   ##
===========================================
- Coverage   86.24%   63.24%   -23.01%     
===========================================
  Files          32       33        +1     
  Lines        3904     3961       +57     
  Branches      947      942        -5     
===========================================
- Hits         3367     2505      -862     
- Misses        346     1316      +970     
+ Partials      191      140       -51     
```


| [Impacted Files](https://codecov.io/gh/bids-standard/pybids/pull/918?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard) | Coverage Δ | |
|---|---|---|
| [bids/modeling/report.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9yZXBvcnQucHk=) | `3.22% <3.22%> (ø)` | |
| [bids/modeling/statsmodels.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9zdGF0c21vZGVscy5weQ==) | `4.53% <5.88%> (-81.46%)` | :arrow_down: |
| [bids/cli.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9jbGkucHk=) | `6.61% <9.09%> (-50.67%)` | :arrow_down: |
| [bids/modeling/auto\_model.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9hdXRvX21vZGVsLnB5) | `0.00% <0.00%> (-87.50%)` | :arrow_down: |
| [bids/modeling/\_\_init\_\_.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9fX2luaXRfXy5weQ==) | `25.00% <0.00%> (-75.00%)` | :arrow_down: |
| [bids/modeling/transformations/base.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy90cmFuc2Zvcm1hdGlvbnMvYmFzZS5weQ==) | `21.42% <0.00%> (-65.97%)` | :arrow_down: |
| [bids/modeling/model\_spec.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9tb2RlbF9zcGVjLnB5) | `26.92% <0.00%> (-54.62%)` | :arrow_down: |
| [bids/modeling/hrf.py](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard#diff-Ymlkcy9tb2RlbGluZy9ocmYucHk=) | `15.15% <0.00%> (-52.53%)` | :arrow_down: |
| ... and [8 more](https://codecov.io/gh/bids-standard/pybids/pull/918/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard) | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/bids-standard/pybids/pull/918?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=bids-standard).

adelavega: As an example, this is an output of `run_graph`:

```
{'run': {'outputs': [<BIDSStatsModelsNodeOutput(level=run, entities={'run': 1, 'subject': '001'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 2, 'subject': '001'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 3, 'subject': '001'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 4, 'subject': '001'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 1, 'subject': '002'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 2, 'subject': '002'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 3, 'subject': '002'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 4, 'subject': '002'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 1, 'subject': '003'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 2, 'subject': '003'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 3, 'subject': '003'})>,
   <BIDSStatsModelsNodeOutput(level=run, entities={'run': 4, 'subject': '003'})>],
  'children': {'subject': {'outputs': [<BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'gain', 'subject': '001'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'loss', 'subject': '001'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'trials', 'subject': '001'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'gain', 'subject': '002'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'loss', 'subject': '002'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'trials', 'subject': '002'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'gain', 'subject': '003'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'loss', 'subject': '003'})>,
     <BIDSStatsModelsNodeOutput(level=subject, entities={'contrast': 'trials', 'subject': '003'})>],
    'children': {'positive': {'outputs': [<BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'gain', 'group': 'equalIndifference'})>,
       <BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'loss', 'group': 'equalIndifference'})>,
       <BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'trials', 'group': 'equalIndifference'})>,
       <BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'gain', 'group': 'equalRange'})>,
       <BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'loss', 'group': 'equalRange'})>,
       <BIDSStatsModelsNodeOutput(level=positive, entities={'contrast': 'trials', 'group': 'equalRange'})>],
      'children': {}},
     'negative-loss': {'outputs': [<BIDSStatsModelsNodeOutput(level=negative-loss, entities={'contrast': 'loss', 'group': 'equalIndifference'})>,
       <BIDSStatsModelsNodeOutput(level=negative-loss, entities={'contrast': 'loss', 'group': 'equalRange'})>],
      'children': {}},
     'between-groups': {'outputs': [<BIDSStatsModelsNodeOutput(level=between-groups, entities={'contrast': 'loss'})>],
      'children': {}}}}}}}
```

I'm starting to think it would make more sense to allows `BIDSStatsModelsNode` objects to be run in a different mode where they store their own ouputs in `outputs_`, and also recursively run their children. 

That way, everything could just be stored with the original `BIDSStatsModelsGraph` object, but in a pre-populated manner. 
adelavega: I think regardless of how we decide to store these outputs, next to dos are:

- [ ] Add `run_report` mode to Output objects
- [ ] Create HTML report 
- [ ] Generate w/ CLI
effigies: @adelavega Do you know why the tests are failing? I don't see anything obvious in your code...
adelavega: Looks like `nilearn` is a requirement for reports but not used anywhere else in pybids. I can add to reqs if desired
adelavega: - Fill out a template (default to pybids template, but fitlins could pass custom template + extra vars)
adelavega: @effigies This is the new API for running a complete graph:

```
graph.load_collections()
graph.run_graph()

graph.root_node.outputs_
>>[<BIDSStatsModelsNodeOutput(name=run, entities={'run': 1, 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 2, 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 3, 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 4, 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 1, 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 2, 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 3, 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 4, 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 1, 'subject': '003'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 2, 'subject': '003'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 3, 'subject': '003'})>,
 <BIDSStatsModelsNodeOutput(name=run, entities={'run': 4, 'subject': '003'})>]

graph.root_node.children[0].destination.outputs_
>>[<BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'loss', 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'trials', 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'gain', 'subject': '001'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'loss', 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'trials', 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'gain', 'subject': '002'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'loss', 'subject': '003'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'trials', 'subject': '003'})>,
 <BIDSStatsModelsNodeOutput(name=subject, entities={'contrast': 'gain', 'subject': '003'})>]
```

If running a single node, `run` still returns a list of `BIDSStatsModelsNodeOutput` objects, but `run_graph` stores them for you under `outputs_` in the `BIDSStatsModelsNodeOutput` object.
adelavega: If running with `run_report=True` the `BIDSStatsModelsNodeOutput` will contain altair plots, such as:

![image](https://user-images.githubusercontent.com/2774448/201794184-cbf6ba95-57bd-4d77-b74c-f93b3d167a7a.png)

Q: Perhaps these should be in a `reports` dictionary object to group them together?

adelavega: @jmumford I addded your VIF code to this. if you pass `run_reports=True` to `.run()` of a Node, you will have `VIF` as an attribute with the contrast VIFs:

![image](https://user-images.githubusercontent.com/2774448/201800580-2010f752-e339-4943-9691-7d71eafc6679.png)

I modified your functions somewhat, so that `get_all_contrast_vif` takes a `BIDSStatsModelNodeOutput` object as input, and `get_eff_reg_vif` takes the vector representation `contrast_def` as input, instead of the string representation, since we don't actually use string representation in `statsmodels`. Instead, I convert the compress weights found in `ContrastInfo` (e.g. only the weights for non-zero items), into a complete vector prior to calling `get_eff_reg_vif`. 

adelavega: This is now generating a basic HTML report copied from FitLins:

TODO:
- [ ] Modify HTML report (e.g. add VIF, remove FitLins stuff)
- [ ] Add altair plots
- [ ] Make it possible to extend template (modify jinja environment loading)

adelavega: This is almost ready to go

Here's an example report (NARPS): https://adelavega.github.io/example_report.html

cc: @jmumford @effigies 
effigies: Wow, that looks great!
adelavega: Feedback:

- VIF cutoffs
- VIF only for t-test or, subset of F-test
- Add note about high pass filtering in VIF section?
- Illustrate frequency of regressors for future high pass filtering
- Only make plot for non binary matrices?
- Add transformation history
- Add CLI flag to make report only for first node

FitLins related feedback:
- Don't automatically add high-ass filering (e.g cosine)

Patch:
diff --git a/MANIFEST.in b/MANIFEST.in
index 84924737..e7e08f6e 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -4,3 +4,4 @@ recursive-include bids/reports/config *.json
 recursive-include bids/tests/data *
 recursive-include bids/layout/tests/data *
 include long_description.rst
+include bids/modeling/report/report_template.jinja
\ No newline at end of file
diff --git a/bids/cli.py b/bids/cli.py
index 6417e667..2c932444 100644
--- a/bids/cli.py
+++ b/bids/cli.py
@@ -7,6 +7,7 @@ import click
 from . import __version__
 from .layout import BIDSLayoutIndexer, BIDSLayout
 from .utils import validate_multiple as _validate_multiple
+from .modeling.report.base import generate_report
 
 # alias -h to trigger help message
 CONTEXT_SETTINGS = {'help_option_names': ['-h', '--help']}
@@ -242,3 +243,31 @@ def upgrade_filenames(root, description):
                 os.rename(bidsfile.path, new_path)
             else:
                 click.echo(f"Not renaming {bidsfile.path}")
+
+
+@cli.command(context_settings=CONTEXT_SETTINGS)
+@click.argument('model', type=click.Path(file_okay=True, exists=True, dir_okay=False))
+@click.argument('root', type=click.Path(file_okay=False, exists=True))
+@click.argument('derivatives', type=click.Path(file_okay=False, exists=True), nargs=-1)
+@click.argument('output_dir', type=click.Path(file_okay=False, exists=False))
+@click.option('--scan_length', type=click.FLOAT)
+@click.option('--sub_label', '-s', multiple=True, 
+    help='The label(s) of the subjects(s) to create report for')
+@click.option('--ses_label', '-ss', multiple=True, 
+    help='The label(s) of the session(s) to create report for')
+@click.option('--task_label', '-t', multiple=True,
+    help='The label(s) of the task(s) to create report for')
+@click.option('--run_label', '-r', multiple=True,
+    help='The label(s) of the run(s) to create report for')
+def model_report(model, root, derivatives, output_dir, scan_length, 
+    sub_label, ses_label, task_label, run_label):
+    """
+    Generate a report of a BIDS StatsModel.
+    """
+    entities = {'subject': sub_label, 'session': ses_label,
+                'task': task_label, 'run': run_label}
+
+    entities = {key:val for key, val in entities.items() if val}
+    generate_report(
+            model, root, derivatives, output_dir, scan_length, **entities
+    )
\ No newline at end of file
diff --git a/bids/modeling/report/__init__.py b/bids/modeling/report/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/bids/modeling/report/base.py b/bids/modeling/report/base.py
new file mode 100644
index 00000000..2c3a31e8
--- /dev/null
+++ b/bids/modeling/report/base.py
@@ -0,0 +1,168 @@
+import os
+import datetime
+from pathlib import Path
+from bids.layout import BIDSLayout
+from bids.modeling import BIDSStatsModelsGraph
+from bids import __version__ as bids_version
+import pkg_resources as pkgr
+from .utils import deroot, snake_to_camel, displayify, to_alphanum, generate_contrast_matrix
+
+PATH_PATTERNS = [
+    'reports/[sub-{subject}/][ses-{session}/][level-{level}_][sub-{subject}_][ses-{session}_]'
+    '[run-{run}_]model-{model}.html'
+]
+
+def _build_node_dict(node, all_entities):
+    report_node = {
+        'name': node.name, 
+        'analyses': [],
+        'group_by': node.group_by,
+        'model': node.model,
+        'level': node.level,
+        'transformations': node.transformations,
+        }
+    for out in node.outputs_:
+        analysis_dict = {'entities': {}, 'contrasts': []}
+
+        for contrast_info in out.contrasts:
+            cents = contrast_info.entities.copy()
+            cents["level"] = out.node.level
+            cents["name"] = out.node.name
+
+            for key in ('datatype', 'desc', 'suffix', 'extension'):
+                cents.pop(key, None)
+            for key in all_entities:
+                cents.pop(key, None)
+
+            for k, v in cents.items():
+                if k in ("name", "contrast"):
+                    cents.update({k: to_alphanum(str(v))})
+
+            analysis_dict['entities'] = {
+                key: val
+                for key, val in cents.items()
+                if key in ('subject', 'session', 'task', 'run') and val
+            }
+
+            analysis_dict['contrasts'].append(
+                {
+                    'name': displayify(contrast_info.name),
+                }
+            )
+
+        ents = out.entities.copy()
+        # Space doesn't apply to design/contrast matrices, or resolution
+        for k in ['space', 'res']:
+            ents.pop(k, None)
+
+        analysis_dict['X'] = out.X
+
+        # If reports were run-level
+        if out.report_ is not None:
+            analysis_dict['design_matrix_plot'] = out.report_['design_matrix_plot'].to_json()
+            analysis_dict['design_matrix_corrplot'] = out.report_['design_matrix_corrplot'].to_json()
+            analysis_dict['VIF'] = out.report_['VIF']
+            analysis_dict['trans_hist'] = out.trans_hist
+            
+        analysis_dict['contrast_matrix'] = generate_contrast_matrix(out.contrasts, out.X.columns)
+        report_node['analyses'].append(analysis_dict)
+
+    return report_node
+
+
+def _build_report_dict(graph):
+    report = {
+        'dataset': {
+            'name': graph.layout.description['Name'],
+        },
+        'model': graph.model,
+        'nodes': [],
+        'version': bids_version,
+        'timestamp': datetime.datetime.now(),
+        'graph_plot': graph.write_graph(format='svg', pipe=True)
+    }
+
+    if 'DatasetDOI' in graph.layout.description:
+        report['dataset']['doi'] = graph.layout.description['DatasetDOI']
+
+    all_entities = graph.layout.get_entities(metadata=True)
+    for name, node in graph.nodes.items():
+        report['nodes'].append(_build_node_dict(node, all_entities))
+
+    # Get subjects hackily
+    report['subjects'] = sorted(
+        {analysis_dict['entities']['subject'] for analysis_dict in report['nodes'][0]['analyses']}
+    )
+
+    return report
+
+
+def _write_report(report_dict, out_dir, template_path=None):
+    try:
+        import jinja2
+    except ImportError:
+        raise ImportError(
+            "Jinja2 must be installed to generate reports. "
+            "You can install it with pip install jinja2."
+        )   
+
+    if template_path is None:
+        searchpath = pkgr.resource_filename('bids', '/')
+        template_file = 'modeling/report/report_template.jinja'
+    else:
+        searchpath, template_file = os.path.split(template_path)
+
+    env = jinja2.Environment(
+        loader=jinja2.FileSystemLoader(searchpath=searchpath))
+    tpl = env.get_template(template_file)
+
+    model = snake_to_camel(report_dict['model']['name'])
+    target_file = os.path.join(
+        out_dir, f"{model}_report.html"
+    )
+
+    report_dict = deroot(report_dict, os.path.dirname(target_file))
+
+    html = tpl.render(report_dict)
+    Path(target_file).parent.mkdir(parents=True, exist_ok=True)
+    Path(target_file).write_text(html)
+
+
+def generate_report(
+    model, dataset_path, derivatives, output_dir, scan_length=None,
+    **entities):
+    """
+    Generate a report for a model.
+
+    Parameters
+    ----------
+    model : str
+        Path to BIDS-StatsModels model.
+    dataset_path : str
+        Path to BIDS dataset.
+    derivatives : str or list
+        Path to BIDS derivatives.
+    output_dir : str
+        Path to output directory.
+    scan_length : float (optional)
+        Length of scan in seconds. If not provided, will be read from
+        dataset (requires image files to be present).
+    entities : dict (optional)
+        A dictionary of BIDS entities to filter the layout on.
+    """
+    # Initialize BIDSLayout
+    layout = BIDSLayout(dataset_path, derivatives=derivatives)
+
+    # Build BIDSStatsModelsGraph instance
+    graph = BIDSStatsModelsGraph(layout, model)
+
+    graph.load_collections(scan_length=scan_length, **entities)
+
+    # Run entire graph
+    graph.run_graph(
+        scan_length=scan_length, entities=entities, transformation_history=True, 
+        node_reports=True)
+
+    report_dict =  _build_report_dict(graph)
+
+    _write_report(report_dict, output_dir)
\ No newline at end of file
diff --git a/bids/modeling/report/report_template.jinja b/bids/modeling/report/report_template.jinja
new file mode 100644
index 00000000..95f74e58
--- /dev/null
+++ b/bids/modeling/report/report_template.jinja
@@ -0,0 +1,325 @@
+<?xml version="1.0" encoding="utf-8"?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en" data-theme="light">
+
+    <head>
+        <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
+        <title>PyBIDS StatsModel Report</title>
+        <style type="text/css">
+        .warning {
+            border: 1px solid #ffaaaa;
+            background: #ffe8e8;
+            padding: 0.8em;
+        }
+        .vega-embed.has-actions {
+            width: 90%;
+        }
+        .overflow-scroll {
+            overflow-x: hidden;
+            overflow-y: auto;
+            max-height: 700px;
+            padding: 1rem;
+        }
+        .card-spacing {
+            padding: 2rem;
+            margin: 2rem 0;
+            box-shadow: 2px 2px 6px 2px lightgrey;
+        }
+    </style>
+        <!-- Import Vega & Vega-Lite -->
+        <script src="https://cdn.jsdelivr.net/npm/vega@5"></script>
+        <script src="https://cdn.jsdelivr.net/npm/vega-lite@4"></script>
+        <!-- Import vega-embed -->
+        <script src="https://cdn.jsdelivr.net/npm/vega-embed@6"></script>
+        <link rel="stylesheet" href="https://unpkg.com/@picocss/pico@latest/css/pico.min.css">
+            <link rel="stylesheet" href="https://cdn.datatables.net/1.13.1/css/jquery.dataTables.min.css">
+                <link rel="stylesheet" href="https://cdn.datatables.net/buttons/2.3.2/css/buttons.dataTables.min.css">
+                    <script src="https://code.jquery.com/jquery-3.5.1.js"></script>
+                    <script src="https://cdn.datatables.net/1.13.1/js/jquery.dataTables.min.js"></script>
+                    <script src="https://cdn.datatables.net/buttons/2.3.2/js/dataTables.buttons.min.js"></script>
+                    <script src="https://cdn.datatables.net/buttons/2.3.2/js/buttons.html5.min.js"></script>
+    </head>
+
+    <body onload="handleSelectNode('{{ nodes[0].name }}')">
+        <script>
+            $.extend( $.fn.dataTable.defaults, {
+            searching: false,
+            ordering:  false
+            } );
+        </script>
+        <script>
+            function handleSelectNode(element) {
+                const nodeButtonList = document.getElementsByClassName('node_button');
+                // change styling of button that we have selected, but remove styling if that button is already selected
+                for(let nodeButton of nodeButtonList) {
+                    if (nodeButton.id === `${element}_node_button` && !nodeButton.classList.contains('secondary')) {
+                        nodeButton.classList.add('secondary');
+                    } else {
+                        nodeButton.className = 'node_button';
+                    }
+                }
+
+                // make node report visible for the associated button we have clicked, hide if corresponding button was already selected and then clicked
+                const nodeReportList = document.getElementsByClassName('node_report');
+                for(let nodeReport of nodeReportList) {
+                    if (nodeReport.id === `${element}_node_report` && nodeReport.style.display === 'none') {
+                        nodeReport.style.display = 'initial';
+                    } else {
+                        nodeReport.style.display = 'none';
+                    }
+                }
+            }
+        </script>
+        <main class="container-fluid", style="width:70%">
+            <div id="summary">
+                <h3>Stats Model Report</h3>
+                <div class="grid">
+                    <div>
+                            <b>Dataset</b>
+                            :
+                    </div>
+                    <div>
+                            {{ dataset.name }}{% if dataset.doi %} (doi:
+                            <a href="https://doi.org/{{ dataset.doi }}">{{ dataset.doi }}</a>
+                            ){% endif %}
+                    </div>
+                    <div></div>
+                </div>
+                <div class="grid">
+                    <div>
+                            <b>Model</b>
+                            :
+                    </div>
+                    <div>
+                        {{ model.name }}
+                    </div>
+                    <div></div>
+                </div>
+                <div class="grid">
+                    <div>
+                            <b>Participants ({{ subjects|count }})</b>
+                            :
+                    </div>
+                    <div>
+                        {{ subjects|join(', ') }}
+                    </div>
+                    <div></div>
+                </div>
+                <br>
+
+                <details>
+                    <summary class="secondary">Specification</summary>
+                    <pre>
+                        <code class="overflow-scroll">{{ model|tojson(indent=2) }}</code>
+                    </pre>
+                </details>
+
+                <details>
+                    <summary class="secondary">Graph</summary>
+                    <div id="model-graph">
+                    {{ graph_plot }}
+                </div>
+                </details>
+            </div>
+            <div id="model">
+                <h3 style="margin-top: 3rem">Node Reports</h3>
+                <div class="grid" style="margin-bottom: 3rem">
+                    {% for node in nodes %}
+                    <span onclick="handleSelectNode('{{node.name}}')" role="button" id="{{node.name}}_node_button" class="node_button" style="margin-bottom: 10px">{{ node.name }}</span>
+                    {% endfor %}
+                </div>
+
+                <div>
+                    {% for node in nodes %}
+                    {% set node_loop = loop %}
+                    <div style="display: none" id="{{node.name}}_node_report" class="node_report">
+                        <div>
+                            <div style="margin-bottom: 3rem;">
+                                <div class="grid">
+                                    <div>
+                                        <b>Level</b>
+                                    </div>
+                                    <div>{{ node.level }}</div>
+                                    <div></div>
+                                </div>
+                                <div class="grid">
+                                    <div>
+                                        <b>Group By</b>
+                                    </div>
+                                    <div>{{ node.group_by | join(", ") }}</div>
+                                    <div></div>
+                                </div>
+                                <div class="grid">
+                                    <div>
+                                        <b>Model</b>
+                                    </div>
+                                    <div>X: {{ node.model['x'] | join(", ") }}</div>
+                                    <div></div>
+                                </div>
+                                {% if node.model['type'] %}
+                                <div class="grid">
+                                    <div></div>
+                                    <div>Type: {{ node.model['type'] }}</div>
+                                    <div></div>
+                                </div>
+                                {% endif %}
+                            </div>
+                            {% if loop.first %}
+                            <div style="margin-bottom: 1.5rem; color: #bbbbbb">
+                                <i>A report was generated for each node instance. All but the first are collapsed, but each should be inspected.</i>
+                            </div>
+                            {% endif %}
+
+                            {% for analysis in node.analyses %}
+                            {% set analysis_loop = loop %}
+                            <details {% if loop.first %}open{% endif %}>
+                                <summary role="button" class="secondary">
+                                    <b>{{ analysis.entities.items()|map('join', ': ')|map('capitalize')|join(', ')}}</b>
+                                </summary>
+                                <div class="overflow-scroll">
+                                    <article class="card-spacing">
+                                        <h5>Design Matrix</h5>
+                                        {% if 'design_matrix_plot' in analysis %}
+                                        <div id="{{ node.name }}{{ loop.index }}_dm"></div>
+                                        <script type="text/javascript">
+                                            var spec = {{ analysis.design_matrix_plot }};
+                                            vegaEmbed('#{{ node.name }}{{ loop.index }}_dm', spec).then(function (result) {}).catch(console.error);
+                                        </script>
+                                        <details>
+                                            <summary role="button" class="contrast">Design Matrix Table</summary>
+                                            <div>
+                                        {% endif %}
+                                        <figure>{{ analysis.X.to_html(border=0, bold_rows=True, classes=["display", "compact"], table_id="{}{}_dmraw".format(node.name,loop.index)) | safe }}</figure>
+                                        <script>
+                                            $(document).ready( function () {
+                                                $('#{{ node.name }}{{ loop.index }}_dmraw').DataTable({
+                                                    dom: 'Bfrtip',
+                                                    pageLength : 5,
+                                                    buttons: [
+                                                        {
+                                                            extend: 'csvHtml5',
+                                                            text: 'Download TSV',
+                                                            fieldSeparator: '\t',
+                                                        }
+                                                    ]
+                                                });
+                                            });
+                                        </script>
+                                        {% if 'design_matrix_plot' in analysis %}
+                                        </div>
+                                        </details>
+                                        {% endif %}
+                                    </article>
+                                    {% if 'trans_hist' in analysis %}
+                                    <article class="card-spacing">
+                                        <h5>Transformation History</h5>
+                                        {% for trans in analysis.trans_hist %}
+                                        {% if loop.first %}
+                                        <p>Each transformation is applied to the design matrix in the order shown. The final design matrix is used in the model estimation.</p>
+                                        <details>
+                                            <summary>Original Variables</summary>
+                                            {% else %}
+                                            <details>
+                                                <summary>{{ trans.index }}. {{ trans.transformation_name}}</summary>
+                                                <ul class="elem-desc">
+                                                    <li>
+                                                        <strong>Input Columns:</strong>
+                                                        {{ trans.input_cols|join(', ') }}
+                                                    </li>
+                                                    <li>
+                                                        <strong>Arguments:</strong>
+                                                        {{ trans.transformation_kwargs.items()|map('join', ': ')|map('capitalize')|join(', ')}}
+                                                    </li>
+                                                </ul>
+                                        {% endif %}
+                                        {% for var in trans.output.variables|sort() %}
+                                        <div id="{{ node.name }}{{ analysis_loop.index }}{{ trans.index }}{{ var }}_thdiv">
+                                            {{ trans.output.variables[var].to_df(entities=False, condition=False).to_html(index=False, border=0, classes=["display", "compact"], table_id="{}{}{}{}_thvar".format(node.name, analysis_loop.index, trans.index, var.replace('.', '_')))}}
+                                            <script>
+                                                $(document).ready( function () {
+                                                    $('#{{ node.name }}{{ analysis_loop.index }}{{ trans.index }}{{ var.replace('.', '_') }}_thvar').DataTable({
+                                                        dom: 'Bfrtip',
+                                                        pageLength : 5,
+                                                        buttons: [
+                                                            {
+                                                                extend: 'csvHtml5',
+                                                                text: 'Download TSV',
+                                                                fieldSeparator: '\t',
+                                                            }
+                                                        ]
+                                                    });
+                                                });
+                                            </script>
+                                        </div>
+                                        {% endfor %}
+                                        </details>
+                                        {% endfor %}
+                                    </article>
+                                    {% endif %}
+                                    {% if 'design_matrix_corrplot' in analysis %}
+                                    <article class="card-spacing">
+                                        <h5>Correlation Matrix</h5>
+                                        <div class="grid">
+                                            <div id="{{ node.name }}{{ loop.index }}_corrdm"></div>
+                                            <div>
+                                                {% if node_loop.first %}
+                                                <p>
+                                                    <strong>Pairwise correlation between regressors.</strong>
+                                                    Very high or low correlations among variables of interest or between variables of interest and nuisance regressors can indicate deficiency in the design. High correlations among nuisance regressors will generally have little effect on the model. It is recommended to interpret these correlations in the context of the VIF for each contrast.
+                                                </p>
+                                                {% endif %}
+                                            </div>
+                                        </div>
+                                        <script type="text/javascript">
+                                            var spec = {{ analysis.design_matrix_corrplot }};
+                                            vegaEmbed('#{{ node.name }}{{ loop.index }}_corrdm', spec).then(function (result) {}).catch(console.error);
+                                        </script>
+                                    </article>
+                                    {% endif %}
+                                    <article class="card-spacing">
+                                        <h4>Contrasts</h4>
+                                        {% if node_loop.first %}
+                                        <i style="color: #bbbbbb">A contrast matrix was generated for each {{ node.name }}. Except in rare cases, these should be identical.</i>
+                                        {% endif %}
+                                        <figure>{{ analysis.contrast_matrix.style.to_html(border=0) | safe }}</figure>
+                                    </article>
+                                    {% if 'VIF' in analysis %}
+                                    <article class="card-spacing">
+                                        <h5>Variance Inflation Factor</h5>
+                                        <div class="grid">
+                                            <div>
+                                                <figure>{{ analysis.VIF.style.hide(axis='index').to_html(border=0) | safe }}</figure>
+                                            </div>
+                                            <div>
+                                                {% if node_loop.first %}
+                                                <p>
+                                                    <strong>Estimated VIF for each contrast. </strong>
+                                                    This is done by extending the effective regressor definition from
+                                                    <em>Smith et al (2007) Meaningful design and contrast estimability (NeuroImage).</em>
+                                                    Regressors involved in the contrast estimate are rotated to span the same space as the original space consisting of the effective regressor and and an orthogonal basis.  The rest of the regressors are unchanged. VIF values below 5 are generally considered acceptable, but values above 10 are problematic. VIF can be used to determine if high pass filter range is appropriate.  If the VIF for a contrast is high, it may be necessary to increase the high pass filter cutoff.  If the VIF for a contrast is low, it may be possible to decrease the high pass filter cutoff.
+                                                </p>
+                                                {% endif %}
+                                            </div>
+                                        </div>
+                                    </article>
+                                    {% endif %}
+                                </div>
+                            </details>
+                            {% endfor %}
+                        </div>
+                    </div>
+                    {% endfor %}
+                </div>
+            </div>
+
+            <div id="about">
+                <h2 style="margin-top: 3rem" class="sub-report-title">About</h2>
+                <ul>
+                    <li>PyBIDS version: {{ version }}</li>
+                    <li>Date processed: {{ timestamp }}</li>
+                </ul>
+            </div>
+        </main>
+    </body>
+
+</html>
\ No newline at end of file
diff --git a/bids/modeling/report/utils.py b/bids/modeling/report/utils.py
new file mode 100644
index 00000000..5e743052
--- /dev/null
+++ b/bids/modeling/report/utils.py
@@ -0,0 +1,199 @@
+import numpy as np
+import pandas as pd
+import re
+from scipy.linalg import null_space
+import warnings
+
+def node_report(node_output):
+    """
+    Generate a report for a single node.
+    """
+    # If only intercept model, don't run reports
+    if node_output.node.level != 'run':
+        return {}
+
+    _report = {
+        'VIF': get_all_contrast_vif(node_output),
+    }
+    try:
+        from .viz import plot_design_matrix, plot_corr_matrix
+    except ImportError:
+        warnings.warn(
+            'altair failed to import, and is required for StatsModel report plots.', 
+            ImportWarning)
+    else:
+        _report['design_matrix_plot'] = plot_design_matrix(
+            node_output.X, timecourse=True)
+
+        _report['design_matrix_corrplot'] = plot_corr_matrix(node_output.X)
+
+    return _report
+
+
+def est_vif(desmat):
+    '''
+    General variance inflation factor estimation.  Calculates VIF for all 
+    regressors in the design matrix.
+
+    Parameters
+    ----------
+        desmat (DataFrame): design matrix. Intercept not required.
+
+    Returns
+    -------
+        vif_data (DataFrame): Variance inflation factor for each regressor in the design matrix
+                generally goal is VIF<5
+    '''
+    from statsmodels.stats.outliers_influence import variance_inflation_factor
+    desmat_with_intercept = desmat.copy()
+    desmat_with_intercept['intercept'] = 1
+    vif_data = pd.DataFrame()
+    vif_data['regressor'] = desmat_with_intercept.columns.drop('intercept')
+    vif_data['VIF'] = [variance_inflation_factor(desmat_with_intercept.values, i)
+                          for i in range(len(desmat_with_intercept.columns))
+                          if desmat_with_intercept.columns[i] != 'intercept']
+    return vif_data
+
+
+def get_eff_reg_vif(desmat, contrast_def):
+    '''
+    The goal of this function is to estimate a variance inflation factor for a contrast.
+    This is done by extending the effective regressor definition from Smith et al (2007)
+    Meaningful design and contrast estimability (NeuroImage).  Regressors involved
+    in the contrast estimate are rotated to span the same space as the original space
+    consisting of the effective regressor and and an orthogonal basis.  The rest of the 
+    regressors are unchanged.
+
+    Parameters
+    ----------
+        desmat (DataFrame): design matrix.  Assumed to be a pandas dataframe with column  
+             headings which are used define the contrast of interest
+        contrast_def (array or Series): a single contrast defined as a vector
+
+    Returns
+    -------
+        vif (float): a single VIF for the contrast of interest  
+    '''
+    des_nuisance_regs = desmat[desmat.columns[contrast_def == 0]]
+    des_contrast_regs = desmat[desmat.columns[contrast_def != 0]]
+
+    con = np.atleast_2d(contrast_def[contrast_def != 0])
+    con2_t = null_space(con)
+    con_t = np.transpose(con)
+    x = des_contrast_regs.copy().values
+    q = np.linalg.pinv(np.transpose(x)@ x)
+    f1 = np.linalg.pinv(con @ q @ con_t)
+    pc = con_t @ f1 @ con @ q
+    con3_t = con2_t - pc @ con2_t
+    f3 = np.linalg.pinv(np.transpose(con3_t) @ q @ con3_t)
+    eff_reg = x @ q @ np.transpose(con) @ f1
+    eff_reg = pd.DataFrame(eff_reg, columns = [0])
+
+    other_reg = x @ q @ con3_t @ f3 
+    other_reg_names = [f'orth_proj{val}' for val in range(other_reg.shape[1])]
+    other_reg = pd.DataFrame(other_reg, columns = other_reg_names)
+
+    des_for_vif = pd.concat([eff_reg, other_reg, des_nuisance_regs], axis = 1)
+    vif_dat = est_vif(des_for_vif)
+    vif_output = vif_dat[vif_dat.regressor == 0].VIF.values[0]
+    return vif_output
+
+
+def generate_contrast_matrix(contrasts, cols):
+    """Generate a contrast matrix from a list of contrast definitions.
+
+    Parameters
+    ----------
+    contrasts : list of ContrastInfo objects
+    cols : list of design matrix columns
+
+    Returns
+    -------
+    contrast_matrix : 2D array
+
+    """
+
+    rows = []
+    ix = []
+    for con in contrasts:
+        vec = cols.map(dict(zip(con.conditions, con.weights))).fillna(0)
+        rows.append(vec.tolist())
+        ix.append(con.name)
+
+    df = pd.DataFrame(rows, columns=cols.tolist(), index=ix)
+    return df
+
+
+def get_all_contrast_vif(node_output):
+    '''
+    Calculates the VIF for multiple contrasts
+
+    Parameters
+    ----------
+        node_output (BIDSStatsModelNodeOutput): Node output to compute VIF for.
+
+    Returns
+    -------
+        vif_contrasts (DataFrame): Data frame containing the VIFs for all contrasts.
+    '''
+    vif_contrasts = {'contrast': [],
+                      'VIF': []}
+    
+    con_matrix = generate_contrast_matrix(
+        node_output.contrasts, node_output.X.columns)
+    for name, weights in con_matrix.iterrows():
+        # Transform weights to vector matching X's columns
+        vif_out = get_eff_reg_vif(node_output.X, weights)
+        vif_contrasts['contrast'].append(name)
+        vif_contrasts['VIF'].append(vif_out) 
+    vif_contrasts = pd.DataFrame(vif_contrasts)
+    return vif_contrasts     
+
+
+def deroot(val, root):
+    if isinstance(val, str):
+        if val.startswith(root):
+            idx = len(root)
+            if val[idx] == '/':
+                idx += 1
+            val = val[idx:]
+    elif isinstance(val, list):
+        val = [deroot(elem, root) for elem in val]
+    elif isinstance(val, dict):
+        val = {key: deroot(value, root) for key, value in val.items()}
+
+    return val
+
+def snake_to_camel(string):
+    words = string.split('_')
+    return words[0] + ''.join(word.title() for word in words[1:])
+
+def displayify(contrast_name):
+    for match, repl in (('_gt_', ' &gt; '), ('_lt_', ' &lt; '), ('_vs_', ' vs. ')):
+        contrast_name = contrast_name.replace(match, repl)
+    return contrast_name
+
+def to_alphanum(string):
+    """Convert string to alphanumeric
+
+    Replaces all other characters with underscores and then converts to camelCase
+
+    Examples
+    --------
+
+    >>> to_alphanum('abc123')
+    'abc123'
+    >>> to_alphanum('a four word phrase')
+    'aFourWordPhrase'
+    >>> to_alphanum('hyphen-separated')
+    'hyphenSeparated'
+    >>> to_alphanum('object.attribute')
+    'objectAttribute'
+    >>> to_alphanum('array[index]')
+    'arrayIndex'
+    >>> to_alphanum('array[0]')
+    'array0'
+    >>> to_alphanum('snake_case')
+    'snakeCase'
+    """
+    return snake_to_camel(re.sub("[^a-zA-Z0-9]", "_", string))
\ No newline at end of file
diff --git a/bids/modeling/report/viz.py b/bids/modeling/report/viz.py
new file mode 100644
index 00000000..8bdc7e77
--- /dev/null
+++ b/bids/modeling/report/viz.py
@@ -0,0 +1,104 @@
+import altair as alt
+import numpy as np
+
+def plot_corr_matrix(dm_wide):
+    dm_corr = dm_wide.corr()
+    dm_corr = dm_corr.where(
+        np.triu(np.ones(dm_corr.shape)).astype(np.bool)).reset_index()
+    dm_corr_long = dm_corr.melt(
+        var_name='index_2',
+        value_vars=dm_wide.columns,  id_vars='index', value_name='r')
+
+    corr_mat = alt.Chart(dm_corr_long).mark_rect().encode(
+        alt.X('index',  sort=None, axis=alt.Axis(title=None)),
+        alt.Y('index_2', sort=None, axis=alt.Axis(title=None)),
+        tooltip=['r'],
+        fill=alt.Color(
+            'r:Q', sort='descending', legend=alt.Legend(title='r'),
+            scale=alt.Scale(scheme='redyellowblue', domain=[-1, 1]))
+    ).properties(
+        width=375,
+        height=350
+    ).configure_scale(
+        bandPaddingInner=0.015
+    ).configure_view(
+        strokeWidth=0
+    ).interactive()
+
+    return corr_mat
+
+alt.data_transformers.enable('default', max_rows=None)
+
+
+def melt_dm(dm):
+    dm = dm.reset_index().rename(columns={'index': 'scan_number'})
+    return dm.melt('scan_number', var_name='regressor', value_name='value')
+
+
+def plot_design_matrix(dm_wide, scale=False, timecourse=True):
+    if scale:
+        dm_wide = (dm_wide - dm_wide.mean()) / (dm_wide.max() - dm_wide.min())
+    dm = melt_dm(dm_wide)
+
+    time_labels = list(range(0, dm.scan_number.max(), 50))
+
+    base_color = alt.Color(
+        'value:Q', sort='ascending',
+        legend=None, scale=alt.Scale(scheme='viridis'))
+
+    heat = alt.Chart(dm).mark_rect().encode(
+        alt.Y('scan_number:O',
+              axis=alt.Axis(title='Time (TRs)', ticks=False,
+                            values=time_labels, labels=True)),
+        alt.X('regressor:N', sort=None, axis=alt.Axis(
+            labelAngle=-45, title=None, ticks=False)),
+        fill=base_color,
+        stroke=base_color
+    ).properties(
+        width=750,
+        height=450,
+    )
+
+    selection = alt.selection_multi(fields=['regressor'])
+    color = alt.Color('regressor:N', legend=None)
+
+    if timecourse:
+        line = alt.Chart(
+            dm,
+            title='Timecourse (shift-click legend to select)').mark_line(
+                clip=True).encode(
+            alt.X('scan_number',
+                axis=alt.Axis(
+                    title='Time (TRs)', values=time_labels, ticks=False),
+                scale=alt.Scale(domain=(0, int(dm.scan_number.max()))),
+                ),
+            y=alt.Y('value:Q',
+                    axis=alt.Axis(title='Amplitude (scaled for visualization)')),
+            color=color
+        ).transform_filter(
+            selection
+        ).properties(
+            width=650,
+            height=225
+        )
+
+    legend = alt.Chart(dm).mark_square(size=200).encode(
+        y=alt.Y('regressor:N', axis=alt.Axis(orient='right')),
+        color=color
+    ).add_selection(
+        selection
+    )
+
+    resolve = alt.Resolve(scale=alt.LegendResolveMap(color=alt.ResolveMode('independent')))
+    if timecourse:
+        plt = alt.vconcat(
+            heat,
+            line | legend,
+            resolve=resolve)
+    else:
+        plt = alt.vconcat(
+            heat,
+            legend,
+            resolve=resolve)
+
+    return plt.configure_scale(bandPaddingInner=0.0)
\ No newline at end of file
diff --git a/bids/modeling/statsmodels.py b/bids/modeling/statsmodels.py
index 9e5a2233..77028607 100644
--- a/bids/modeling/statsmodels.py
+++ b/bids/modeling/statsmodels.py
@@ -15,6 +15,7 @@ from bids.utils import matches_entities, convert_JSON, listify
 from bids.variables import (BIDSVariableCollection, merge_collections)
 from bids.modeling import transformations as tm
 from .model_spec import GLMMSpec, MetaAnalysisSpec
+from .report.utils import node_report
 import warnings
 
 
@@ -201,7 +202,7 @@ class BIDSStatsModelsGraph:
                                                       **node_kwargs)
             node.add_collections(collections)
 
-    def write_graph(self, dotfilename='graph.dot', format='png'):
+    def write_graph(self, dotfilename='graph.dot', format='png', pipe=False):
         """Generates a graphviz dot file and a png file
 
         Parameters
@@ -218,7 +219,7 @@ class BIDSStatsModelsGraph:
                 filename=dotfilename,
                 node_attr={'shape': 'record'},
                 comment=self.model['name'],
-                format=format
+                format=format, 
             )
 
         for node, nobj in self.nodes.items():
@@ -227,10 +228,47 @@ class BIDSStatsModelsGraph:
         for edge in self.edges:
             dot.edge(edge['source'], edge['destination'])
 
-        dot.render()
+        if pipe:
+            dot = dot.pipe(encoding='utf-8')
+        else:
+            dot.render()
 
         return dot
 
+    def run_graph(self, entities=None, **kwargs):
+        """Run the entire graph recursively.
+
+        Parameters
+        ----------
+        entities : dict (optional)
+            Optional dictionary of BIDS entities to use when loading collections.
+        kwargs : dict
+            Optional dictionary of keyword arguments to pass to
+            BIDSStatsModelsNode.run()
+        """
+        if entities is None:
+            entities = {}
+        
+        _run_node_recursive(self.root_node, filters=entities, **kwargs)
+
+def _run_node_recursive(node, inputs=None, filters=None, **kwargs):
+    """
+    Run a node recursively, storing outputs in place as
+    BIDSStatsModelsNode.outputs_.
+
+    """
+    if filters == None:
+        filters = {}
+        
+    # Run node
+    run_kwargs = {**filters, **kwargs}
+    node.outputs_ =  node.run(inputs, group_by=node.group_by, **run_kwargs)
+
+    # Inputs to next node
+    contrasts = list(itertools.chain(*[s.contrasts for s in node.outputs_]))
+
+    for edge in node.children:
+        _run_node_recursive(edge.destination, contrasts, filters=edge.filter, **kwargs)
 
 class BIDSStatsModelsNode:
     """Represents a single node in a BIDS-StatsModel graph.
@@ -398,8 +436,8 @@ class BIDSStatsModelsNode:
         return groups
 
     def run(self, inputs=None, group_by=None, force_dense=True,
-              sampling_rate='TR', invalid_contrasts='drop',
-              transformation_history=False, **filters):
+              sampling_rate='TR', invalid_contrasts='drop', 
+              transformation_history=False, node_reports=False, **filters):
         """Execute node with provided inputs.
 
         Parameters
@@ -439,6 +477,8 @@ class BIDSStatsModelsNode:
         transformation_history: bool
             If True, the returned ModelSpec instances will include a history of
             variable collections after each transformation.
+        node_reports: bool
+            If True, report measures and plots will be included in the output
         filters: dict
             Optional keyword arguments used to constrain the subset of the data
             that's processed. E.g., passing subject='01' will process and
@@ -478,7 +518,7 @@ class BIDSStatsModelsNode:
                 node=self, entities=dict(grp_ents), collections=grp_colls,
                 inputs=grp_inputs, force_dense=force_dense,
                 sampling_rate=sampling_rate, invalid_contrasts=invalid_contrasts,
-                transformation_history=transformation_history)
+                transformation_history=transformation_history, node_reports=node_reports)
             results.append(node_output)
 
         return results
@@ -582,10 +622,12 @@ class BIDSStatsModelsNodeOutput:
     transformation_history: bool
         If True, the returned ModelSpec instances will include a history of
         variable collections after each transformation.
+    node_reports: bool
+        If True, a report will be generated for each node output.
     """
     def __init__(self, node, entities={}, collections=None, inputs=None,
                  force_dense=True, sampling_rate='TR', invalid_contrasts='drop',
-                 *, transformation_history=False):
+                 *, transformation_history=False, node_reports=False):
         """Initialize a new BIDSStatsModelsNodeOutput instance.
         Applies the node's model to the specified collections and inputs, including
         applying transformations and generating final model specs and design matrices (X).
@@ -653,6 +695,8 @@ class BIDSStatsModelsNodeOutput:
         self.model_spec = SpecCls.from_df(self.data, node.model, self.metadata)
         self.contrasts = self._build_contrasts(unique_in_contrast)
 
+        self.report_ = node_report(self) if node_reports else None
+
     def _collections_to_dfs(self, collections, *, collection_history=False):
         """Merges collections and converts them to a pandas DataFrame."""
         if not collections:
@@ -802,4 +846,4 @@ class BIDSStatsModelsNodeOutput:
         return self.model_spec.X
 
     def __repr__(self):
-        return f"<{self.__class__.__name__}(level={self.node.name}, entities={self.entities})>"
+        return f"<{self.__class__.__name__}(name={self.node.name}, entities={self.entities})>"
diff --git a/pyproject.toml b/pyproject.toml
index 1c3d4e17..9cfaba63 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -50,6 +50,11 @@ test = [
   "pytest-cov",
   "bsmschema",
   "coverage[toml]",
+  "altair"
+]
+model_reports = [
+  "jinja2",
+  "altair"
 ]
 tutorial = [
   "nbconvert",


Test patch:
diff --git a/bids/modeling/tests/test_statsmodels.py b/bids/modeling/tests/test_statsmodels.py
index db4c6fce..c66ffcf1 100644
--- a/bids/modeling/tests/test_statsmodels.py
+++ b/bids/modeling/tests/test_statsmodels.py
@@ -54,7 +54,7 @@ def test_repr(graph):
     assert graph.__repr__() == "<BIDSStatsModelsGraph[{name='test_model', description='simple test model', ... }]>"
     node = graph.nodes['run']
     assert node.__repr__() == "<BIDSStatsModelsNode(level=run, name=run)>"
-    assert node.run()[0].__repr__() == "<BIDSStatsModelsNodeOutput(level=run, entities={'run': 1, 'subject': '01'})>"
+    assert node.run()[0].__repr__() == "<BIDSStatsModelsNodeOutput(name=run, entities={'run': 1, 'subject': '01'})>"
 
 def test_manual_intercept(graph_intercept):
     # Test that a automatic intercept (1) is correct
diff --git a/bids/tests/test_cli.py b/bids/tests/test_cli.py
index ae9fec1f..58179aff 100644
--- a/bids/tests/test_cli.py
+++ b/bids/tests/test_cli.py
@@ -65,4 +65,4 @@ def test_layout(runner, tmp_path):
         ],
         catch_exceptions=False,
     )
-    assert is_success(res)
+    assert is_success(res)
\ No newline at end of file