Repository: mne-tools/mne-nirs (language: python)

Issue:
Change function run_GLM to run_glm
According to https://www.python.org/dev/peps/pep-0008/#function-and-variable-names 

"Function names should be lowercase, with words separated by underscores as necessary to improve readability."


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/NumpyDoc_PR.yml b/.github/workflows/NumpyDoc_PR.yml
index fbc9f91cd..8e65939b9 100644
--- a/.github/workflows/NumpyDoc_PR.yml
+++ b/.github/workflows/NumpyDoc_PR.yml
@@ -5,8 +5,8 @@ on:
 
 jobs:
   # PIP + non-default stim channel + log level info
-  job:
-    name: 'build'
+  numpydoc:
+    name: 'numpydoc'
     runs-on: ubuntu-20.04
     defaults:
       run:
diff --git a/.github/workflows/PyLint_PR.yml b/.github/workflows/PyLint_PR.yml
index a3ae235b1..180136ef8 100644
--- a/.github/workflows/PyLint_PR.yml
+++ b/.github/workflows/PyLint_PR.yml
@@ -11,7 +11,7 @@ on:
 jobs:
 
   lint:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-20.04
 
     steps:
       - uses: actions/checkout@v1
diff --git a/.github/workflows/SphinxBuild.yml b/.github/workflows/SphinxBuild.yml
index 720649eac..8ea2ec4ca 100644
--- a/.github/workflows/SphinxBuild.yml
+++ b/.github/workflows/SphinxBuild.yml
@@ -8,7 +8,7 @@ on:
 jobs:
   # PIP + non-default stim channel + log level info
   job:
-    name: 'build and push'
+    name: 'docs and push'
     runs-on: ubuntu-20.04
     defaults:
       run:
diff --git a/.github/workflows/SphinxBuild_PR.yml b/.github/workflows/SphinxBuild_PR.yml
index 1377858a1..59cf6f15a 100644
--- a/.github/workflows/SphinxBuild_PR.yml
+++ b/.github/workflows/SphinxBuild_PR.yml
@@ -6,7 +6,7 @@ on:
 jobs:
   # PIP + non-default stim channel + log level info
   job:
-    name: 'build'
+    name: 'docs'
     runs-on: ubuntu-20.04
     defaults:
       run:
diff --git a/doc/api.rst b/doc/api.rst
index e8bd0f92a..e9eaeabf2 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -118,7 +118,6 @@ Data quality evaluation.
    scalp_coupling_index_windowed
 
 
-
 Experimental Design
 -------------------
 
@@ -135,22 +134,57 @@ Experimental Design
    create_boxcar
 
 
-
 Statistics
 ----------
 
+First level analysis
+^^^^^^^^^^^^^^^^^^^^
+
+Individual (first) level analysis functions.
+
+.. currentmodule:: mne_nirs.statistics
+
+.. automodule:: mne_nirs.statistics
+   :no-members:
+   :no-inherited-members:
+
+.. autosummary::
+   :toctree: generated/
+
+   run_glm
+
+
+Individual (first) level result classes.
+
+.. currentmodule:: mne_nirs.statistics
+
+.. automodule:: mne_nirs.statistics
+   :no-members:
+   :no-inherited-members:
+   :noindex:
+
+.. autosummary::
+   :toctree: generated/
+
+   RegressionResults
+   ContrastResults
+
+
+Second level analysis
+^^^^^^^^^^^^^^^^^^^^^
+
+Group (second) level analysis functions.
+
 .. currentmodule:: mne_nirs.statistics
 
 .. automodule:: mne_nirs.statistics
    :no-members:
    :no-inherited-members:
+   :noindex:
 
 .. autosummary::
    :toctree: generated/
 
-   run_GLM
-   compute_contrast
-   glm_region_of_interest
    statsmodels_to_results
 
 
@@ -253,6 +287,3 @@ General helper functions.
    :toctree: generated/
 
    glm_to_tidy
-
-
-
diff --git a/doc/changelog.rst b/doc/changelog.rst
index e3b42fff0..74001349a 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -15,6 +15,29 @@ Patch version changes indicate backward compatible bug fixes.
 v0.1.0
 ------
 
+API changes
+
+* Add :class:`~mne_nirs.statistics.RegressionResults` and :class:`~mne_nirs.statistics.ContrastResults` classes to store GLM results. By `Robert Luke`_.
+
+Adding a class simplifies user code and common use cases.
+To generate results in the new format, use the function run_glm rather than run_GLM.
+This will return a RegressionResults type that contains all relevant information.
+All previous functionality still exists with this new type,
+but is now accessible as more succinct methods that handle the relevant information,
+this results in less arguments being passed around by the user.
+For example, to access the previous glm_to_tidy(results) functionality use the new results.to_dataframe().
+A full list of replacement methods is provided below.
+
+Enhancements
+
+* Add :meth:`mne_nirs.statistics.RegressionResults.compute_contrast` to replace compute_contrast.
+* Add :meth:`mne_nirs.statistics.RegressionResults.plot_topo` to replace plot_glm_topo.
+* Add :meth:`mne_nirs.statistics.RegressionResults.to_dataframe` to replace glm_to_tidy.
+* Add :meth:`mne_nirs.statistics.RegressionResults.to_dataframe_region_of_interest` to replace glm_region_of_interest.
+* Add :meth:`mne_nirs.statistics.RegressionResults.scatter` to display GLM results as a scatter plot.
+* Add :meth:`mne_nirs.statistics.RegressionResults.surface_projection` to display GLM results on a cortical surface.
+
+
 Infrastructure
 
 * Change the git repository to use `main`, rather than `master` branch.
diff --git a/doc/conf.py b/doc/conf.py
index b568e20ec..b66834e80 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -26,7 +26,6 @@
 from mne.tests.test_docstring_parameters import error_ignores
 
 
-
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
@@ -65,6 +64,7 @@
 
 # generate autosummary even if no references.
 autosummary_generate = True
+autodoc_default_options = {'inherited-members': None}
 
 # The suffix of source filenames.
 source_suffix = '.rst'
diff --git a/examples/general/plot_10_hrf_simulation.py b/examples/general/plot_10_hrf_simulation.py
index bc620ad41..06eb463f1 100644
--- a/examples/general/plot_10_hrf_simulation.py
+++ b/examples/general/plot_10_hrf_simulation.py
@@ -26,7 +26,7 @@
 import matplotlib.pylab as plt
 import numpy as np
 from mne_nirs.experimental_design import make_first_level_design_matrix
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm
 from nilearn.plotting import plot_design_matrix
 np.random.seed(1)
 
@@ -75,11 +75,17 @@
 # which was the amplitude we used in the simulation.
 # We also see that the mean square error of the model fit is close to zero.
 
-glm_est = run_GLM(raw, design_matrix)
+glm_est = run_glm(raw, design_matrix)
+
+
+def print_results(glm_est, truth):
+    """Function to print the results of GLM estimate"""
+    print("Estimate:", glm_est.theta()[0][0],
+          "  MSE:", glm_est.MSE()[0],
+          "  Error (uM):", 1e6*(glm_est.theta()[0][0] - truth * 1e-6))
+
+print_results(glm_est, amp)
 
-print("Estimate:", glm_est['Simulated'].theta[0],
-      "  MSE:", glm_est['Simulated'].MSE,
-      "  Error (uM):", 1e6*(glm_est['Simulated'].theta[0] - amp*1e-6))
 
 
 ###############################################################################
@@ -98,18 +104,16 @@
 raw_noise_free = raw.copy()
 
 raw._data += np.random.normal(0, np.sqrt(1e-11), raw._data.shape)
-glm_est = run_GLM(raw, design_matrix)
+glm_est = run_glm(raw, design_matrix)
 
 plt.plot(raw.times, raw_noise_free.get_data().T * 1e6)
 plt.plot(raw.times, raw.get_data().T * 1e6, alpha=0.3)
-plt.plot(raw.times, glm_est['Simulated'].theta[0] * design_matrix["A"].values * 1e6)
+plt.plot(raw.times, glm_est.theta()[0][0] * design_matrix["A"].values * 1e6)
 plt.xlabel("Time (s)")
 plt.ylabel("Haemoglobin (uM)")
 plt.legend(["Clean Data", "Noisy Data", "GLM Estimate"])
 
-print("Estimate:", glm_est['Simulated'].theta[0],
-      "  MSE:", glm_est['Simulated'].MSE,
-      "  Error (uM):", 1e6*(glm_est['Simulated'].theta[0] - amp*1e-6))
+print_results(glm_est, amp)
 
 
 ###############################################################################
@@ -132,18 +136,16 @@
 design_matrix = make_first_level_design_matrix(raw, stim_dur=5.0,
                                                drift_order=1,
                                                drift_model='polynomial')
-glm_est = run_GLM(raw, design_matrix)
+glm_est = run_glm(raw, design_matrix)
 
 plt.plot(raw.times, raw_noise_free.get_data().T * 1e6)
 plt.plot(raw.times, raw.get_data().T * 1e6, alpha=0.3)
-plt.plot(raw.times, glm_est['Simulated'].theta[0] * design_matrix["A"].values * 1e6)
+plt.plot(raw.times, glm_est.theta()[0][0] * design_matrix["A"].values * 1e6)
 plt.xlabel("Time (s)")
 plt.ylabel("Haemoglobin (uM)")
 plt.legend(["Clean Data", "Noisy Data", "GLM Estimate"])
 
-print("Estimate:", glm_est['Simulated'].theta[0],
-      "  MSE:", glm_est['Simulated'].MSE,
-      "  Error (uM):", 1e6*(glm_est['Simulated'].theta[0] - amp*1e-6))
+print_results(glm_est, amp)
 
 
 ###############################################################################
@@ -167,17 +169,15 @@
 design_matrix = make_first_level_design_matrix(raw, stim_dur=5.0,
                                                drift_order=1,
                                                drift_model='polynomial')
-glm_est = run_GLM(raw, design_matrix)
+glm_est = run_glm(raw, design_matrix)
 
 plt.plot(raw.times, raw.get_data().T * 1e6, alpha=0.3)
-plt.plot(raw.times, glm_est['Simulated'].theta[0] * design_matrix["A"].values * 1e6)
+plt.plot(raw.times, glm_est.theta()[0][0] * design_matrix["A"].values * 1e6)
 plt.xlabel("Time (s)")
 plt.ylabel("Haemoglobin (uM)")
 plt.legend(["Noisy Data", "GLM Estimate"])
 
-print("Estimate:", glm_est['Simulated'].theta[0],
-      "  MSE:", glm_est['Simulated'].MSE,
-      "  Error (uM):", 1e6*(glm_est['Simulated'].theta[0] - amp*1e-6))
+print_results(glm_est, amp)
 
 
 ###############################################################################
@@ -191,12 +191,12 @@
 # properties was extracted from the data and if this
 # improved the response estimate.
 
-glm_est = run_GLM(raw, design_matrix, noise_model='ar5')
+glm_est = run_glm(raw, design_matrix, noise_model='ar5')
 
 fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(15, 6))
 plt.plot([-0.58853134, -0.29575669, -0.52246482, 0.38735476, 0.02428681],
          axes=axes)  # actual values from model above
-plt.plot(glm_est['Simulated'].model.rho * -1.0,  axes=axes)  # estimates
+plt.plot(glm_est.model()[0].rho * -1.0,  axes=axes)  # estimates
 plt.legend(["Simulation AR coefficients", "Estimated AR coefficients"])
 plt.xlabel("Coefficient")
 
@@ -205,9 +205,7 @@
 # We can see that the estimates from the GLM AR model are quite accurate,
 # but how does this affect the accuracy of the response estimate?
 
-print("Estimate:", glm_est['Simulated'].theta[0],
-      "  MSE:", glm_est['Simulated'].MSE,
-      "  Error (uM):", 1e6*(glm_est['Simulated'].theta[0] - amp*1e-6))
+print_results(glm_est, amp)
 
 ###############################################################################
 # The response estimate using the AR(5) model is more accurate than the
diff --git a/examples/general/plot_11_hrf_measured.py b/examples/general/plot_11_hrf_measured.py
index 60604f836..3d52d60cb 100644
--- a/examples/general/plot_11_hrf_measured.py
+++ b/examples/general/plot_11_hrf_measured.py
@@ -43,7 +43,7 @@
 import mne_nirs
 
 from mne_nirs.experimental_design import make_first_level_design_matrix
-from mne_nirs.statistics import run_GLM, glm_region_of_interest
+from mne_nirs.statistics import run_glm
 from mne_nirs.visualisation import plot_glm_topo
 from mne_nirs.channels import (get_long_channels, get_short_channels,
                                picks_pair_to_idx)
@@ -92,8 +92,7 @@
 raw_intensity.annotations.rename({'1.0': 'Control',
                                   '2.0': 'Tapping/Left',
                                   '3.0': 'Tapping/Right'})
-raw_intensity.annotations.delete(
-    np.where([d == '15.0' for d in raw_intensity.annotations.description]))
+raw_intensity.annotations.delete(raw_intensity.annotations.description == '15.0')
 raw_intensity.annotations.set_durations(5)
 
 
@@ -133,10 +132,8 @@
 # We observe that the order of conditions was randomised and the time between
 # events is also randomised.
 
-events, _ = mne.events_from_annotations(raw_haemo, verbose=False)
-event_dict = {'Control': 1, 'Tapping/Left': 2, 'Tapping/Right': 3}
-mne.viz.plot_events(events, event_id=event_dict,
-                    sfreq=raw_haemo.info['sfreq'])
+events, event_dict = mne.events_from_annotations(raw_haemo, verbose=False)
+mne.viz.plot_events(events, event_id=event_dict, sfreq=raw_haemo.info['sfreq'])
 
 
 ###############################################################################
@@ -247,42 +244,58 @@
 # of a single source detector pair.
 
 data_subset = raw_haemo.copy().pick(picks=range(2))
-glm_est = run_GLM(data_subset, design_matrix)
+glm_est = run_glm(data_subset, design_matrix)
 
 ###############################################################################
 #
 # This returns a GLM regression estimate for each channel.
-# This data is stored in a dictionary, which can be addressed using the
-# desired channel name. You can view an overview of the estimates as:
+# This data is stored in a dedicated type.
+# You can view an overview of the estimates by addressing the variable:
 
 glm_est
 
+###############################################################################
+#
+# As with other MNE types you can use the `pick` function.
+# To query the mean square error of a single channel you would call.
+#
+# Note: as we wish to retain both channels for further the analysis below,
+# we operate on a copy to demonstrate this channel picking functionality.
+
+glm_est.copy().pick('S1_D1 hbr')
 
 ###############################################################################
 #
-# Or you can view the estimate for a single channel by indexing it by name.
-# You can see that for each channel a standard
+# Underlying the data for each channel is a standard
 # `Nilearn RegressionResults object <https://nilearn.github.io/modules/generated/nilearn.glm.RegressionResults.html>`_
-# is returned. These objects are rich with information that can be requested
+# object. These objects are rich with information that can be requested
 # from the object, for example to determine the mean square error of the
-# estimate you would call
+# estimates for two channels you would call:
+
+glm_est.MSE()
 
-glm_est['S1_D1 hbo'].MSE
+###############################################################################
+#
+# And we can chain the methods to quickly access required details.
+# For example, to determine the MSE for channel `S1` `D1` for the hbr type
+# you would call:
+
+glm_est.copy().pick('S1_D1 hbr').MSE()
 
 
 ###############################################################################
 #
-# Due to the richness of the object we provide a function `glm_to_tidy` to
+# Due to the richness of the objects we provide a function to
 # extract commonly used information and put it in a convenient dataframe/table.
 # Below this is demonstrated and then we just display the first 9 rows of the
 # table which correspond to the 9 components of the design matrix for the
 # first channel.
 
-glm_to_tidy(data_subset, glm_est, design_matrix).head(9)
+glm_est.to_dataframe().head(9)
 
 ###############################################################################
 #
-# We then display the results. In this example we address the objects directly.
+# We then display the results using the scatter plot function.
 # Note that the control condition sits
 # around zero
 # and that the HbO is positive and larger than the HbR, this is to be expected.
@@ -290,12 +303,7 @@
 # right hand is larger than the left. And the values are similar to what
 # is seen in the epoching tutorial.
 
-plt.scatter(design_matrix.columns[:3], glm_est['S1_D1 hbo'].theta[:3] * 1e6)
-plt.scatter(design_matrix.columns[:3], glm_est['S1_D1 hbr'].theta[:3] * 1e6)
-plt.xlabel("Experiment Condition")
-plt.ylabel("Haemoglobin (μM)")
-plt.legend(["Oxyhaemoglobin", "Deoxyhaemoglobin"])
-plt.hlines([0.0], 0, 2)
+glm_est.scatter()
 
 
 ###############################################################################
@@ -309,10 +317,8 @@
 # contralateral to the tapping side. Also note that HbR tends to be the
 # negative of HbO as expected.
 
-glm_est = run_GLM(raw_haemo, design_matrix)
-plot_glm_topo(raw_haemo, glm_est, design_matrix,
-              requested_conditions=['Tapping/Left',
-                                    'Tapping/Right'])
+glm_est = run_glm(raw_haemo, design_matrix)
+glm_est.plot_topo(conditions=['Tapping/Left', 'Tapping/Right'])
 
 
 ###############################################################################
@@ -337,23 +343,28 @@
 # apparent that the data does not indicate that activity spreads across
 # the center of the head.
 
-fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 6),
-                         gridspec_kw=dict(width_ratios=[0.92, 1]))
+fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 6), gridspec_kw=dict(width_ratios=[0.92, 1]))
+
+glm_hbo = glm_est.copy().pick(picks="hbo")
+conditions = ['Tapping/Right']
+
+glm_hbo.plot_topo(axes=axes[0], colorbar=False, conditions=conditions)
 
-plot_glm_topo(raw_haemo.copy().pick(picks="hbo"), glm_est, design_matrix,
-              axes=axes[0], colorbar=False,
-              requested_conditions=['Tapping/Right'])
-plot_glm_topo(raw_haemo.copy().pick(picks="hbo").pick(picks=range(10)),
-              glm_est, design_matrix, vmin=-16, vmax=16, axes=axes[1],
-              requested_conditions=['Tapping/Right'], colorbar=False)
-plot_glm_topo(raw_haemo.copy().pick(picks="hbo").pick(picks=range(10, 20)),
-              glm_est, design_matrix, axes=axes[1], vmin=-16, vmax=16,
-              requested_conditions=['Tapping/Right'])
+glm_hbo.copy().pick(picks=range(10)).plot_topo(conditions=conditions, axes=axes[1], colorbar=False, vmin=-16, vmax=16)
+glm_hbo.copy().pick(picks=range(10, 20)).plot_topo(conditions=conditions, axes=axes[1], colorbar=False, vmin=-16, vmax=16)
 
 axes[0].set_title("Smoothed across hemispheres")
 axes[1].set_title("Hemispheres plotted independently")
 
 
+###############################################################################
+#
+# Another way to view the data is to project the GLM estimates to the nearest
+# cortical surface
+
+glm_est.copy().surface_projection(condition="Tapping/Right", view="dorsal", chroma="hbo")
+
+
 ###############################################################################
 # Analyse regions of interest
 # ---------------------------
@@ -383,9 +394,10 @@
 groups = dict(Left_ROI=picks_pair_to_idx(raw_haemo, left),
               Right_ROI=picks_pair_to_idx(raw_haemo, right))
 
-df = pd.DataFrame()
-for idx, col in enumerate(design_matrix.columns[:3]):
-    df = df.append(glm_region_of_interest(glm_est, groups, idx, col))
+conditions = ['Control', 'Tapping/Left', 'Tapping/Right']
+
+df = glm_est.to_dataframe_region_of_interest(groups, conditions)
+
 
 ###############################################################################
 # As with the single channel results above, this is placed in a tidy dataframe
@@ -394,6 +406,7 @@
 
 df
 
+
 ###############################################################################
 #
 # Compute contrasts
@@ -409,8 +422,9 @@
 basic_conts = dict([(column, contrast_matrix[i])
                    for i, column in enumerate(design_matrix.columns)])
 contrast_LvR = basic_conts['Tapping/Left'] - basic_conts['Tapping/Right']
-contrast = mne_nirs.statistics.compute_contrast(glm_est, contrast_LvR)
-mne_nirs.visualisation.plot_glm_contrast_topo(raw_haemo, contrast)
+
+contrast = glm_est.compute_contrast(contrast_LvR)
+contrast.plot_topo()
 
 
 ###############################################################################
@@ -429,7 +443,7 @@
 # if a significant response (p<0.05) was observed, which sensor and detector is
 # in the channel, which chroma, etc.
 
-df = glm_to_tidy(raw_haemo, glm_est, design_matrix)
+df = glm_est.to_dataframe()
 
 
 ###############################################################################
diff --git a/examples/general/plot_12_group_glm.py b/examples/general/plot_12_group_glm.py
index 91412c634..8d22232b5 100644
--- a/examples/general/plot_12_group_glm.py
+++ b/examples/general/plot_12_group_glm.py
@@ -80,9 +80,9 @@
 from mne.preprocessing.nirs import optical_density, beer_lambert_law
 
 # Import MNE-NIRS processing
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm
 from mne_nirs.experimental_design import make_first_level_design_matrix
-from mne_nirs.statistics import glm_region_of_interest, statsmodels_to_results
+from mne_nirs.statistics import statsmodels_to_results
 from mne_nirs.statistics import compute_contrast
 from mne_nirs.channels import get_short_channels, get_long_channels
 from mne_nirs.channels import picks_pair_to_idx
@@ -152,7 +152,7 @@ def individual_analysis(bids_path, ID):
     design_matrix["ShortHbR"] = np.mean(sht_chans.copy().pick(picks="hbr").get_data(), axis=0)
 
     # Run GLM
-    glm_est = run_GLM(raw_haemo, design_matrix)
+    glm_est = run_glm(raw_haemo, design_matrix)
 
     # Define channels in each region of interest
     # List the channel pairs manually
@@ -164,23 +164,23 @@ def individual_analysis(bids_path, ID):
         Right_Hemisphere=picks_pair_to_idx(raw_haemo, right, on_missing='ignore'))
 
     # Extract channel metrics
-    cha = glm_to_tidy(raw_haemo, glm_est, design_matrix)
-    cha["ID"] = ID  # Add the participant ID to the dataframe
+    cha = glm_est.to_dataframe()
 
     # Compute region of interest results from channel data
-    roi = pd.DataFrame()
-    for idx, col in enumerate(design_matrix.columns):
-        roi = roi.append(glm_region_of_interest(glm_est, groups, idx, col))
-    roi["ID"] = ID  # Add the participant ID to the dataframe
+    roi = glm_est.to_dataframe_region_of_interest(groups, design_matrix.columns)
 
-    # Contrast left vs right tapping
+    # Define left vs right tapping contrast
     contrast_matrix = np.eye(design_matrix.shape[1])
     basic_conts = dict([(column, contrast_matrix[i])
                         for i, column in enumerate(design_matrix.columns)])
     contrast_LvR = basic_conts['Tapping/Left'] - basic_conts['Tapping/Right']
-    contrast = compute_contrast(glm_est, contrast_LvR)
-    con = glm_to_tidy(raw_haemo, contrast, design_matrix)
-    con["ID"] = ID  # Add the participant ID to the dataframe
+
+    # Compute defined contrast
+    contrast = glm_est.compute_contrast(contrast_LvR)
+    con = contrast.to_dataframe()
+
+    # Add the participant ID to the dataframes
+    roi["ID"] = cha["ID"] = con["ID"] = ID
 
     # Convert to uM for nicer plotting below.
     cha["theta"] = [t * 1.e6 for t in cha["theta"]]
diff --git a/examples/general/plot_13_fir_glm.py b/examples/general/plot_13_fir_glm.py
index 1885a41e9..f95fb5314 100644
--- a/examples/general/plot_13_fir_glm.py
+++ b/examples/general/plot_13_fir_glm.py
@@ -60,7 +60,7 @@
 from mne.preprocessing.nirs import optical_density, beer_lambert_law
 
 # Import MNE-NIRS processing
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm
 from mne_nirs.experimental_design import make_first_level_design_matrix
 from mne_nirs.statistics import glm_region_of_interest, statsmodels_to_results
 from mne_nirs.datasets import fnirs_motor_group
@@ -113,14 +113,14 @@ def analysis(fname, ID):
         design_matrix[f"short_{chan}"] = short_chans.get_data(chan).T
 
     # Run GLM
-    glm_est = run_GLM(raw_haemo, design_matrix)
+    glm_est = run_glm(raw_haemo, design_matrix)
 
     # Create a single ROI that includes all channels for example
     rois = dict(AllChannels=range(len(raw_haemo.ch_names)))
+    # Calculate ROI for all conditions
+    conditions = design_matrix.columns
     # Compute output metrics by ROI
-    df_ind = pd.DataFrame()
-    for idx, col in enumerate(design_matrix.columns):
-        df_ind = df_ind.append(glm_region_of_interest(glm_est, rois, idx, col))
+    df_ind = glm_est.to_dataframe_region_of_interest(rois, conditions)
 
     df_ind["ID"] = ID
     df_ind["theta"] = [t * 1.e6 for t in df_ind["theta"]]
diff --git a/mne_nirs/statistics/__init__.py b/mne_nirs/statistics/__init__.py
index d3584bb95..ed0955a2f 100644
--- a/mne_nirs/statistics/__init__.py
+++ b/mne_nirs/statistics/__init__.py
@@ -1,9 +1,8 @@
-"""Statistical analysis."""
-
 # Authors: Robert Luke <mail@robertluke.net>
 #
 # License: BSD (3-clause)
 
-from ._run_GLM import run_GLM, compute_contrast
 from ._roi import glm_region_of_interest
 from ._statsmodels import statsmodels_to_results
+from ._glm_level_first import (RegressionResults, ContrastResults,
+                               run_glm, compute_contrast, run_GLM)
diff --git a/mne_nirs/statistics/_glm_level_first.py b/mne_nirs/statistics/_glm_level_first.py
new file mode 100644
index 000000000..33340b1d6
--- /dev/null
+++ b/mne_nirs/statistics/_glm_level_first.py
@@ -0,0 +1,622 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+from copy import deepcopy
+
+import pandas as pd
+import numpy as np
+from numpy import array_equal, where
+
+import nilearn.glm
+from nilearn.glm.first_level import run_glm as nilearn_glm
+
+from mne.channels.channels import ContainsMixin
+from mne.utils import fill_doc, warn, verbose
+from mne.io.pick import _picks_to_idx
+
+from ..visualisation._plot_GLM_topo import _plot_glm_topo,\
+    _plot_glm_contrast_topo
+from ..visualisation import plot_glm_surface_projection
+from ..statistics._roi import _glm_region_of_interest
+
+
+@fill_doc
+class _BaseGLM(ContainsMixin):
+    """Base GLM results class."""
+
+    @property
+    def ch_names(self):
+        """Return the channel names.
+
+        Returns
+        -------
+        names : array
+            The channel names.
+        """
+        return self.info['ch_names']
+
+    def __str__(self):
+        return (f"GLM Results for {len(self.ch_names)} channels")
+
+    def __repr__(self):
+        return (f"GLM Results for {len(self.ch_names)} channels")
+
+    def __len__(self):
+        return len(self.info.ch_names)
+
+    def copy(self):
+        """Return a copy of the GLM results.
+
+        Returns
+        -------
+        inst : instance of ResultsGLM
+            A copy of the object.
+        """
+        return deepcopy(self)
+
+    def to_dataframe(self, order=None):
+        """Return a tidy dataframe representing the GLM results.
+
+        Parameters
+        ----------
+        order : list
+            Order in which the rows should be returned by channel name.
+
+        Returns
+        -------
+        tidy : pandas.DataFrame
+            Dataframe containing GLM results.
+        """
+        from ..utils import glm_to_tidy
+        return glm_to_tidy(self.info, self._data, self.design, order=order)
+
+    def scatter(self, conditions=[], exclude_no_interest=True, axes=None,
+                no_interest=None):
+        """Scatter plot of the GLM results.
+
+        Parameters
+        ----------
+        conditions : list
+            List of condition names to plot. By default plots all regressors
+            of interest.
+        exclude_no_interest : bool
+            Exclude regressors of no interest from the figure.
+        axes : Axes
+            Optional axes on which to plot the data.
+        no_interest : list
+            List of regressors that are of no interest. If none are specified
+            then conditions starting with
+            ["drift", "constant", "short", "Short"] will be excluded.
+
+        Returns
+        -------
+        plt : matplotlib.Figure
+            Scatter plot.
+        """
+        if no_interest is None:
+            no_interest = ["drift", "constant", "short", "Short"]
+        import matplotlib.pyplot as plt
+        df = self.to_dataframe()
+
+        x_column = "Condition"
+        y_column = "theta"
+        if "ContrastType" in df.columns:
+            x_column = "ch_name"
+            y_column = "effect"
+            if len(conditions) == 0:
+                conditions = ["t", "f"]
+            df = df.query('ContrastType in @conditions')
+
+        else:
+            if len(conditions) == 0:
+                conditions = self.design.columns
+            df = df.query('Condition in @conditions')
+
+        if exclude_no_interest:
+            for no_i in no_interest:
+                df = df[~df[x_column].astype(str).str.startswith(no_i)]
+
+        df_hbo = df.query('Chroma == "hbo"')
+        df_hbr = df.query('Chroma == "hbr"')
+
+        if axes is None:
+            fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 5))
+
+        axes.scatter(df_hbo[x_column], df_hbo[y_column] * 1e6, c="r")
+        axes.scatter(df_hbr[x_column], df_hbr[y_column] * 1e6, c="b")
+        axes.set_xlabel(x_column)
+        axes.set_ylabel(y_column.capitalize())
+        axes.legend(["Oxyhaemoglobin", "Deoxyhaemoglobin"])
+        axes.hlines([0.0], 0, len(np.unique(df[x_column])) - 1)
+        if len(np.unique(df[x_column])) > 8:
+            plt.xticks(rotation=45, ha='right')
+
+        return axes
+
+
+@fill_doc
+class RegressionResults(_BaseGLM):
+    """
+    Class containing GLM regression results.
+
+    Parameters
+    ----------
+    info : mne.Info
+        Info.
+    data : dict
+        Dictionary.
+    design : dataframe
+        Design matrix.
+
+    Returns
+    -------
+    glm_est : ResultsGLM,
+        Result class.
+    """
+
+    @property
+    def data(self):
+        return self._data
+
+    @data.setter
+    def data(self, data):
+        if type(data) is not dict:
+            raise TypeError("Data must be a dictionary type")
+        if not array_equal(list(data.keys()), self.info.ch_names):
+            raise TypeError("Dictionary keys must match ch_names")
+        for idx in range(len(self.info.ch_names)):
+            if list(data.keys())[idx] is not self.info.ch_names[idx]:
+                raise TypeError("Data names and channel names do not match")
+        for d in data:
+            if type(data[d]) is not nilearn.glm.regression.RegressionResults:
+                raise TypeError("Dictionary items must be"
+                                " nilearn RegressionResults")
+
+        self._data = data
+
+    def __init__(self, info, data, design):  # noqa: D102
+        self.info = info
+        self.data = data
+        self.design = design
+        self.preload = True
+
+    @fill_doc
+    def pick(self, picks, exclude=()):
+        """Pick a subset of channels.
+
+        Parameters
+        ----------
+        %(picks_all)s
+        exclude : list | str
+            Set of channels to exclude, only used when picking based on
+            types (e.g., exclude="bads" when picks="meg").
+
+        Returns
+        -------
+        inst : instance of ResultsGLM
+            The modified instance.
+        """
+        picks = _picks_to_idx(self.info, picks, 'all', exclude,
+                              allow_empty=False)
+        picks = [self.ch_names[p] for p in picks]
+        self.info.pick_channels(picks)
+        self._data = {key: self._data[key] for key in self.info.ch_names}
+        return self
+
+    def theta(self):
+        """Return the GLM theta results.
+
+        Returns
+        -------
+        thetas : array
+            Array of thetas. An array is provided per channel.
+        """
+        return [self._data[key].theta for key in self.info.ch_names]
+
+    def MSE(self):
+        """Return the GLM MSE.
+
+        Returns
+        -------
+        thetas : array
+            Array of MSEs. A value is provided per channel.
+        """
+        return [self._data[key].MSE[0] for key in self.info.ch_names]
+
+    def model(self):
+        """Return the GLM model.
+
+        Returns
+        -------
+        models : array
+            Array of models. An model is provided per channel.
+        """
+        return [self._data[key].model for key in self.info.ch_names]
+
+    def compute_contrast(self, contrast, contrast_type=None):
+        """
+        Compute contrasts on regression results.
+
+        This is a wrapper function for nilearn.stats.contrasts.
+
+        Parameters
+        ----------
+        contrast : numpy.ndarray of shape (p) or (q, p),
+            Where q = number of contrast vectors and p = number of regressors.
+        contrast_type : {None, ‘t’, ‘F’}, optional
+            Type of the contrast. If None, then defaults to ‘t’ for 1D con_val
+            and ‘F’ for 2D con_val.
+
+        Returns
+        -------
+        contrast : Contrast instance,
+            Yields the statistics of the contrast
+            (effects, variance, p-values).
+        """
+        cont = _compute_contrast(self._data, contrast,
+                                 contrast_type=contrast_type)
+        return ContrastResults(self.info, cont, self.design)
+
+    def plot_topo(self, conditions=None,
+                  axes=None, vmin=None, vmax=None, colorbar=True,
+                  figsize=(12, 7), sphere=None):
+        """Plot 2D topography of GLM data.
+
+        Parameters
+        ----------
+        conditions : array
+            Which conditions should be displayed.
+        axes : instance of Axes | None
+            The axes to plot to. If None, a new figure is used.
+        vmin : float | None
+            The value specifying the lower bound of the color range.
+            If None, and vmax is None, -vmax is used. Else np.min(data).
+            Defaults to None.
+        vmax : float | None
+            The value specifying the upper bound of the color range.
+            If None, the maximum absolute value is used. Defaults to None.
+        colorbar : Bool
+            Should a colorbar be plotted.
+        figsize : two values
+            Figure size.
+        sphere : As specified in MNE
+            Sphere parameter from mne.viz.topomap.plot_topomap.
+
+        Returns
+        -------
+        fig : figure
+            Figure of each design matrix componenent for hbo (top row)
+            and hbr (bottom row).
+        """
+        return _plot_glm_topo(self.info, self._data, self.design,
+                              requested_conditions=conditions,
+                              axes=axes, vmin=vmin, vmax=vmax,
+                              colorbar=colorbar,
+                              figsize=figsize, sphere=sphere)
+
+    def to_dataframe_region_of_interest(self, group_by, condition,
+                                        weighted=True):
+        """Calculate region of interest results and return as a dataframe.
+
+        Parameters
+        ----------
+        group_by : dict
+            Specifies which channels are aggregated into a single ROI.
+            The dict key will be used as the ROI label and the dict
+            values must be lists of picks (either channel names or
+            integer indices of ``epochs.ch_names``). For example::
+
+                group_by=dict(Left_ROI=[1, 2, 3, 4], Right_ROI=[5, 6, 7, 8])
+
+            Note that within a dict entry all channels must have the same type.
+        condition : str | list
+            Name to be used for condition.
+        weighted : Bool
+            Should channels be weighted by inverse of standard error (True).
+
+        Returns
+        -------
+        tidy : pandas.DataFrame
+            Dataframe containing GLM results.
+        """
+        if isinstance(condition, str):
+            condition = [condition]
+
+        tidy = pd.DataFrame()
+        for cond in condition:
+            cond_idx = where([c == cond for c in self.design.columns])[0]
+
+            roi = _glm_region_of_interest(self._data, group_by,
+                                          cond_idx, cond, weighted)
+            tidy = tidy.append(roi)
+
+        return tidy
+
+    @verbose
+    def surface_projection(self, chroma="hbo", condition=None,
+                           background='w', figure=None, clim='auto',
+                           mode='weighted', colormap='RdBu_r',
+                           surface='pial', hemi='both', size=800,
+                           view=None, colorbar=True, distance=0.03,
+                           subjects_dir=None, src=None, verbose=False):
+        """
+        Project GLM results on to the surface of the brain.
+
+        Note: This function provides a convenient wrapper around low level
+        MNE-Python functions. It is convenient if you wish to use a generic
+        head model.
+        If you have acquired fMRI images you may wish to use the underlying
+        lower level functions.
+
+        Note: This function does not conduct a forward model analysis with
+        photon migration etc. It simply projects the values from each channel
+        to the local cortical surface.
+        It is useful for visualisation, but users should
+        take this in to consideration when drawing conclusions from the
+        visualisation.
+
+        Parameters
+        ----------
+        chroma : str
+            Chromophore to plot`.
+        condition : str
+            Condition to plot`.
+        background : matplotlib color
+            Color of the background of the display window.
+        figure : mayavi.core.api.Scene, matplotlib.figure.Figure, list, None
+            If None, a new figure will be created. If multiple views or a
+            split view is requested, this must be a list of the appropriate
+            length. If int is provided it will be used to identify the Mayavi
+            figure by it's id or create a new figure with the given id. If an
+            instance of matplotlib figure, mpl backend is used for plotting.
+        %(clim)s
+        mode : str
+            Can be "sum" to do a linear sum of weights, "weighted" to make this
+            a weighted sum, "nearest" to
+            use only the weight of the nearest sensor, or "single" to
+            do a distance-weight of the nearest sensor. Default is "sum".
+        colormap : str
+            Colormap to use.
+        surface : str
+            The type of surface (inflated, white etc.).
+        hemi : str
+            Hemisphere id (ie 'lh', 'rh', 'both', or 'split'). In the case
+            of 'both', both hemispheres are shown in the same window.
+            In the case of 'split' hemispheres are displayed side-by-side
+            in different viewing panes.
+        size : float or tuple of float
+            The size of the window, in pixels. can be one number to specify
+            a square window, or the (width, height) of a rectangular window.
+            Has no effect with mpl backend.
+        view : str
+            View to set brain to.
+        colorbar : bool
+            If True, display colorbar on scene.
+        distance : float
+            Distance (m) defining the activation "ball" of the sensor.
+        %(subjects_dir)s
+        src : instance of SourceSpaces
+            The source space.
+        %(verbose)s
+
+        Returns
+        -------
+        figure : instance of mne.viz.Brain | matplotlib.figure.Figure
+            An instance of :class:`mne.viz.Brain` or matplotlib figure.
+        """
+
+        df = self.to_dataframe(order=self.ch_names)
+        if condition is None:
+            warn("You must provide a condition to plot", ValueError)
+        df = df.query("Condition in @condition")
+        df = df.query("Chroma in @chroma")
+        df["theta"] = df["theta"] * 1e6
+
+        info = self.copy().pick(chroma).info
+
+        return plot_glm_surface_projection(info, df, value="theta",
+                                           background=background,
+                                           figure=figure, clim=clim,
+                                           mode=mode, colormap=colormap,
+                                           surface=surface, hemi=hemi,
+                                           size=size,
+                                           view=view, colorbar=colorbar,
+                                           distance=distance,
+                                           subjects_dir=subjects_dir, src=src,
+                                           verbose=verbose
+                                           )
+
+
+@fill_doc
+class ContrastResults(_BaseGLM):
+    """
+    Class containing GLM contrast results.
+
+    Parameters
+    ----------
+    info : mne.Info
+        Info.
+    data : dict
+        Dictionary.
+    design : dataframe
+        Design matrix.
+
+    Returns
+    -------
+    glm_est : ResultsGLM,
+        Result class.
+    """
+
+    def __init__(self, info, data, design):  # noqa: D102
+        self.info = info
+        self.data = data
+        self.design = design
+        self.preload = True
+
+    @property
+    def data(self):
+        return self._data
+
+    @data.setter
+    def data(self, data):
+        if not isinstance(data, nilearn.glm.contrasts.Contrast):
+            raise TypeError("Data must be a nilearn glm contrast type")
+        if data.effect.size != len(self.info.ch_names):
+            raise TypeError("Data results must be the same length "
+                            "as the number of channels")
+
+        self._data = data
+
+    def plot_topo(self, figsize=(12, 7), sphere=None):
+        """
+        Plot topomap GLM contrast data.
+
+        Parameters
+        ----------
+        figsize : numbers
+            TODO: Remove this, how does MNE usually deal with this.
+        sphere : numbers
+            As specified in MNE.
+
+        Returns
+        -------
+        fig : figure
+            Figure of each design matrix componenent for hbo (top row)
+            and hbr (bottom row).
+        """
+        return _plot_glm_contrast_topo(self.info, self._data,
+                                       figsize=figsize, sphere=sphere)
+
+
+def run_GLM(raw, design_matrix, noise_model='ar1', bins=0,
+            n_jobs=1, verbose=0):
+    """
+    Run GLM on data using supplied design matrix.
+    This is a wrapper function for nilearn.stats.first_level_model.run_glm.
+    Parameters
+    ----------
+    raw : instance of Raw
+        The haemoglobin data.
+    design_matrix : as specified in Nilearn
+        The design matrix.
+    noise_model : {'ar1', 'ols', 'arN', 'auto'}, optional
+        The temporal variance model. Defaults to first order
+        auto regressive model 'ar1'.
+        The AR model can be set to any integer value by modifying the value
+        of N. E.g. use `ar5` for a fifth order model.
+        If the string `auto` is provided a model with order 4 times the sample
+        rate will be used.
+    bins : int, optional
+        Maximum number of discrete bins for the AR coef histogram/clustering.
+        By default the value is 0, which will set the number of bins to the
+        number of channels, effectively estimating the AR model for each
+        channel.
+    n_jobs : int, optional
+        The number of CPUs to use to do the computation. -1 means
+        'all CPUs'.
+    verbose : int, optional
+        The verbosity level. Default is 0.
+    Returns
+    -------
+    glm_estimates : dict
+        Keys correspond to the different labels values values are
+        RegressionResults instances corresponding to the voxels.
+    """
+    warn('"run_GLM" has been deprecated in favor of the more '
+         'comprehensive run_glm function, and will be removed in v1.0.0. '
+         'See the changelog for further details.',
+         DeprecationWarning)
+    res = run_glm(raw, design_matrix, noise_model=noise_model, bins=bins,
+                  n_jobs=n_jobs, verbose=verbose)
+    return res.data
+
+
+def run_glm(raw, design_matrix, noise_model='ar1', bins=0,
+            n_jobs=1, verbose=0):
+    """
+    GLM fit for an MNE structure containing fNIRS data.
+
+    This is a wrapper function for nilearn.stats.first_level_model.run_glm.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        The haemoglobin data.
+    design_matrix : as specified in Nilearn
+        The design matrix.
+    noise_model : {'ar1', 'ols', 'arN', 'auto'}, optional
+        The temporal variance model. Defaults to first order
+        auto regressive model 'ar1'.
+        The AR model can be set to any integer value by modifying the value
+        of N. E.g. use `ar5` for a fifth order model.
+        If the string `auto` is provided a model with order 4 times the sample
+        rate will be used.
+    bins : int, optional
+        Maximum number of discrete bins for the AR coef histogram/clustering.
+        By default the value is 0, which will set the number of bins to the
+        number of channels, effectively estimating the AR model for each
+        channel.
+    n_jobs : int, optional
+        The number of CPUs to use to do the computation. -1 means
+        'all CPUs'.
+    verbose : int, optional
+        The verbosity level. Default is 0.
+
+    Returns
+    -------
+    glm_estimates : RegressionResults
+        RegressionResults class which stores the GLM results.
+    """
+    picks = _picks_to_idx(raw.info, 'fnirs', exclude=[], allow_empty=True)
+    ch_names = raw.ch_names
+
+    if noise_model == 'auto':
+        noise_model = f"ar{int(np.round(raw.info['sfreq'] * 4))}"
+
+    if bins == 0:
+        bins = len(raw.ch_names)
+
+    results = dict()
+    for pick in picks:
+        labels, glm_estimates = nilearn_glm(raw.get_data(pick).T,
+                                            design_matrix.values,
+                                            noise_model=noise_model, bins=bins,
+                                            n_jobs=n_jobs, verbose=verbose)
+        results[ch_names[pick]] = glm_estimates[labels[0]]
+
+    return RegressionResults(raw.info, results, design_matrix)
+
+
+def compute_contrast(glm_est, contrast, contrast_type=None):
+    """
+    Compute contrasts on regression results.
+
+    This is a wrapper function for nilearn.stats.contrasts.
+
+    Parameters
+    ----------
+    glm_est : dict
+        Dictionary of nilearn regression results as returned by `run_glm`.
+    contrast : numpy.ndarray of shape (p) or (q, p),
+        Where q = number of contrast vectors and p = number of regressors.
+    contrast_type : {None, ‘t’, ‘F’}, optional
+        Type of the contrast. If None, then defaults to ‘t’ for 1D con_val
+        and ‘F’ for 2D con_val.
+
+    Returns
+    -------
+    contrast : Contrast instance,
+        Yields the statistics of the contrast (effects, variance, p-values).
+    """
+    warn('"compute_contrast" has been deprecated in favor of the more '
+         'comprehensive GLM class and will be removed in v1.0.0. '
+         'Use the ResultsGLM class "compute_contrast()" method instead.',
+         DeprecationWarning)
+    return _compute_contrast(glm_est, contrast, contrast_type=contrast_type)
+
+
+def _compute_contrast(glm_est, contrast, contrast_type=None):
+    from nilearn.glm.contrasts import compute_contrast as _cc
+    return _cc(np.array(list(glm_est.keys())), glm_est, contrast,
+               contrast_type=contrast_type)
diff --git a/mne_nirs/statistics/_roi.py b/mne_nirs/statistics/_roi.py
index 7c2576b8e..be343cf18 100644
--- a/mne_nirs/statistics/_roi.py
+++ b/mne_nirs/statistics/_roi.py
@@ -1,11 +1,57 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
 import numpy as np
 
+from mne.utils import warn
+
 
-def glm_region_of_interest(stats, group_by, cond_idx,
+def glm_region_of_interest(glm, group_by, cond_idx,
                            cond_name, weighted=True):
     """
     Calculate statistics for region of interest.
 
+    Parameters
+    ----------
+    glm : dict
+        Need to write.
+    group_by : dict
+        Specifies which channels are aggregated into a single ROI.
+        The dict key will be used as the ROI label and the dict
+        values must be lists of picks (either channel names or integer indices
+        of ``epochs.ch_names``). For example::
+
+            group_by=dict(Left_ROI=[1, 2, 3, 4], Right_ROI=[5, 6, 7, 8])
+
+        Note that within a dict entry all channels must have the same type.
+    cond_name : str
+        Name to be used for condition.
+    cond_idx : int
+        Index of condition of interest.
+    weighted : Bool
+        Should channels be weighted by inverse of standard error (True).
+
+    Returns
+    -------
+    stats : DataFrame
+        Statistics for each ROI.
+    """
+    warn('"glm_region_of_interest" has been deprecated in favor of the more '
+         'comprehensive GLM class and will be removed in v1.0.0. '
+         'Use the RegressionResults class "region_of_interest_dataframe()" '
+         'method instead.',
+         DeprecationWarning)
+
+    return _glm_region_of_interest(glm, group_by,
+                                   cond_idx, cond_name, weighted=weighted)
+
+
+def _glm_region_of_interest(stats, group_by, cond_idx,
+                            cond_name, weighted=True):
+    """
+    Calculate statistics for region of interest.
+
     Parameters
     ----------
     stats : dict
diff --git a/mne_nirs/statistics/_run_GLM.py b/mne_nirs/statistics/_run_GLM.py
deleted file mode 100644
index 1497cf99b..000000000
--- a/mne_nirs/statistics/_run_GLM.py
+++ /dev/null
@@ -1,92 +0,0 @@
-# Authors: Robert Luke <mail@robertluke.net>
-#
-# License: BSD (3-clause)
-
-import numpy as np
-from mne.io.pick import _picks_to_idx
-
-
-def run_GLM(raw, design_matrix, noise_model='ar1', bins=0,
-            n_jobs=1, verbose=0):
-    """
-    Run GLM on data using supplied design matrix.
-
-    This is a wrapper function for nilearn.stats.first_level_model.run_glm.
-
-    Parameters
-    ----------
-    raw : instance of Raw
-        The haemoglobin data.
-    design_matrix : as specified in Nilearn
-        The design matrix.
-    noise_model : {'ar1', 'ols', 'arN', 'auto'}, optional
-        The temporal variance model. Defaults to first order
-        auto regressive model 'ar1'.
-        The AR model can be set to any integer value by modifying the value
-        of N. E.g. use `ar5` for a fifth order model.
-        If the string `auto` is provided a model with order 4 times the sample
-        rate will be used.
-    bins : int, optional
-        Maximum number of discrete bins for the AR coef histogram/clustering.
-        By default the value is 0, which will set the number of bins to the
-        number of channels, effectively estimating the AR model for each
-        channel.
-    n_jobs : int, optional
-        The number of CPUs to use to do the computation. -1 means
-        'all CPUs'.
-    verbose : int, optional
-        The verbosity level. Default is 0.
-
-    Returns
-    -------
-    glm_estimates : dict
-        Keys correspond to the different labels values values are
-        RegressionResults instances corresponding to the voxels.
-    """
-    from nilearn.glm.first_level import run_glm
-
-    picks = _picks_to_idx(raw.info, 'fnirs', exclude=[], allow_empty=True)
-    ch_names = raw.ch_names
-
-    if noise_model == 'auto':
-        noise_model = f"ar{int(np.round(raw.info['sfreq'] * 4))}"
-
-    if bins == 0:
-        bins = len(raw.ch_names)
-
-    results = dict()
-    for pick in picks:
-        labels, glm_estimates = run_glm(raw.get_data(pick).T,
-                                        design_matrix.values,
-                                        noise_model=noise_model, bins=bins,
-                                        n_jobs=n_jobs, verbose=verbose)
-        results[ch_names[pick]] = glm_estimates[labels[0]]
-
-    return results
-
-
-def compute_contrast(glm_est, contrast, contrast_type=None):
-    """
-    Compute contrasts on regression results.
-
-    This is a wrapper function for nilearn.stats.contrasts.
-
-    Parameters
-    ----------
-    glm_est : dict
-        Dictionary of nilearn regression results as returned by `run_glm`.
-    contrast : numpy.ndarray of shape (p) or (q, p),
-        Where q = number of contrast vectors and p = number of regressors.
-    contrast_type : {None, ‘t’, ‘F’}, optional
-        Type of the contrast. If None, then defaults to ‘t’ for 1D con_val
-        and ‘F’ for 2D con_val.
-
-    Returns
-    -------
-    contrast : Contrast instance,
-        Yields the statistics of the contrast (effects, variance, p-values).
-    """
-
-    from nilearn.glm.contrasts import compute_contrast as _cc
-    return _cc(np.array(list(glm_est.keys())), glm_est, contrast,
-               contrast_type=contrast_type)
diff --git a/mne_nirs/utils/_io.py b/mne_nirs/utils/_io.py
index 839668be7..68913121e 100644
--- a/mne_nirs/utils/_io.py
+++ b/mne_nirs/utils/_io.py
@@ -6,26 +6,28 @@
 import nilearn
 
 
-def glm_to_tidy(raw, statistic, design_matrix, wide=True):
+def glm_to_tidy(info, statistic, design_matrix, wide=True, order=None):
     """
     Export GLM regression or contrast results in tidy format.
 
     Creates a long pandas data frame from regression results or contrast
-    as computed by run_GLM or compute_contrast.
+    as computed by run_glm or compute_contrast.
 
     Parameters
     ----------
-    raw : MNE.Raw
-        Instance of MNE raw.
+    info : MNE.Info
+        Instance of MNE info.
     statistic : nilearn data,
         Either dict of nilearn.stats.regression.RegressionResults as returned
-        by run_GLM, or nilearn.stats.contrasts.Contrast as returned by
+        by run_glm, or nilearn.stats.contrasts.Contrast as returned by
         compute_contrast.
     design_matrix : DataFrame
         As specified in Nilearn.
     wide : Bool
         Should the returned dataframe be in wide format. If False, then the
         returned data will be in long format.
+    order : list
+        Order that the channels should be returned with.
 
     Returns
     -------
@@ -36,10 +38,10 @@ def glm_to_tidy(raw, statistic, design_matrix, wide=True):
     if isinstance(statistic, dict) and \
             isinstance(statistic[list(statistic.keys())[0]],
                        nilearn.glm.regression.RegressionResults):
-        df = _tidy_RegressionResults(raw, statistic, design_matrix)
+        df = _tidy_RegressionResults(info, statistic, design_matrix)
 
     elif isinstance(statistic, nilearn.glm.contrasts.Contrast):
-        df = _tidy_Contrast(raw, statistic, design_matrix)
+        df = _tidy_Contrast(info, statistic, design_matrix)
 
     else:
         raise TypeError(
@@ -49,6 +51,15 @@ def glm_to_tidy(raw, statistic, design_matrix, wide=True):
     if wide:
         df = _tidy_long_to_wide(df, expand_output=True)
 
+    if order is not None:
+        df['old_index'] = df.index
+        df = df.set_index('ch_name')
+        df = df.loc[order, :]
+        df['ch_name'] = df.index
+        df.index = df['old_index']
+        df.drop(columns='old_index', inplace=True)
+        df.rename_axis(None, inplace=True)
+
     return df
 
 
diff --git a/mne_nirs/visualisation/_plot_GLM_surface_projection.py b/mne_nirs/visualisation/_plot_GLM_surface_projection.py
index 788a4dba8..f32f0bb9c 100644
--- a/mne_nirs/visualisation/_plot_GLM_surface_projection.py
+++ b/mne_nirs/visualisation/_plot_GLM_surface_projection.py
@@ -4,15 +4,16 @@
 
 import os
 import numpy as np
+from copy import deepcopy
 
-from mne import stc_near_sensors, EvokedArray, read_source_spaces
+from mne import stc_near_sensors, EvokedArray, read_source_spaces, Info
 from mne.datasets import sample
 from mne.io.constants import FIFF
 from mne import verbose
 
 
 @verbose
-def plot_glm_surface_projection(raw, statsmodel_df, picks="hbo",
+def plot_glm_surface_projection(inst, statsmodel_df, picks="hbo",
                                 value="Coef.",
                                 background='w', figure=None, clim='auto',
                                 mode='weighted', colormap='RdBu_r',
@@ -35,7 +36,7 @@ def plot_glm_surface_projection(raw, statsmodel_df, picks="hbo",
 
     Parameters
     ----------
-    raw : instance of Raw
+    inst : instance of Raw
         Haemoglobin data.
     statsmodel_df : dataframe
         As produced by produced by `statsmodels_to_results`.
@@ -85,13 +86,14 @@ def plot_glm_surface_projection(raw, statsmodel_df, picks="hbo",
     figure : instance of mne.viz.Brain | matplotlib.figure.Figure
         An instance of :class:`mne.viz.Brain` or matplotlib figure.
     """
-    if not (raw.ch_names == list(statsmodel_df['ch_name'].values)):
+    info = deepcopy(inst if isinstance(inst, Info) else inst.info)
+    if not (info.ch_names == list(statsmodel_df['ch_name'].values)):
         raise RuntimeError('MNE data structure does not match dataframe '
-                           f'results. Raw = {raw.ch_names}. '
+                           f'results.\nMNE = {info.ch_names}.\n'
                            f'GLM = {list(statsmodel_df["ch_name"].values)}')
 
     ea = EvokedArray(np.tile(statsmodel_df[value].values.T, (1, 1)).T,
-                     raw.copy().info)
+                     info.copy())
 
     # TODO: mimic behaviour of other MNE-NIRS glm plotting options
     if picks is not None:
diff --git a/mne_nirs/visualisation/_plot_GLM_topo.py b/mne_nirs/visualisation/_plot_GLM_topo.py
index f9fa3c39f..87b86b859 100644
--- a/mne_nirs/visualisation/_plot_GLM_topo.py
+++ b/mne_nirs/visualisation/_plot_GLM_topo.py
@@ -2,6 +2,7 @@
 #
 # License: BSD (3-clause)
 
+from copy import deepcopy
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -9,12 +10,13 @@
 from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
 
 import mne
+from mne import Info
 from mne.utils import warn
 from mne.channels.layout import _merge_ch_data
-from mne.io.pick import _picks_to_idx
+from mne.io.pick import _picks_to_idx, _get_channel_types
 
 
-def plot_glm_topo(raw, glm_estimates, design_matrix,
+def plot_glm_topo(inst, glm_estimates, design_matrix,
                   requested_conditions=None,
                   axes=None, vmin=None, vmax=None, colorbar=True,
                   figsize=(12, 7), sphere=None):
@@ -26,8 +28,8 @@ def plot_glm_topo(raw, glm_estimates, design_matrix,
 
     Parameters
     ----------
-    raw : instance of Raw
-        Haemoglobin data.
+    inst : instance of Info or Raw
+        Raw data or info structure used to generate the GLM results.
     glm_estimates : dict
         Keys correspond to the different labels values values are
         RegressionResults instances corresponding to the voxels.
@@ -57,14 +59,30 @@ def plot_glm_topo(raw, glm_estimates, design_matrix,
         Figure of each design matrix componenent for hbo (top row)
         and hbr (bottom row).
     """
+    warn('"plot_glm_topo" has been deprecated in favor of the more '
+         'comprehensive GLM class and will be removed in v1.0.0. '
+         'Use the RegressionResults class "plot_topo()" method instead.',
+         DeprecationWarning)
+    return _plot_glm_topo(inst, glm_estimates, design_matrix,
+                          requested_conditions=requested_conditions,
+                          axes=axes, vmin=vmin, vmax=vmax, colorbar=colorbar,
+                          figsize=figsize, sphere=sphere)
 
-    if not (raw.ch_names == list(glm_estimates.keys())):
-        if len(raw.ch_names) < len(list(glm_estimates.keys())):
+
+def _plot_glm_topo(inst, glm_estimates, design_matrix,
+                   requested_conditions=None,
+                   axes=None, vmin=None, vmax=None, colorbar=True,
+                   figsize=(12, 7), sphere=None):
+
+    info = deepcopy(inst if isinstance(inst, Info) else inst.info)
+
+    if not (info.ch_names == list(glm_estimates.keys())):
+        if len(info.ch_names) < len(list(glm_estimates.keys())):
             warn("Reducing GLM results to match MNE data")
-            glm_estimates = {a: glm_estimates[a] for a in raw.ch_names}
+            glm_estimates = {a: glm_estimates[a] for a in info.ch_names}
         else:
             raise RuntimeError('MNE data structure does not match regression '
-                               f'results. Raw = {len(raw.ch_names)}. '
+                               f'results. Raw = {len(info.ch_names)}. '
                                f'GLM = {len(list(glm_estimates.keys()))}')
 
     estimates = np.zeros((len(glm_estimates), len(design_matrix.columns)))
@@ -72,7 +90,7 @@ def plot_glm_topo(raw, glm_estimates, design_matrix,
     for idx, name in enumerate(glm_estimates.keys()):
         estimates[idx, :] = glm_estimates[name].theta.T
 
-    types = np.unique(raw.get_channel_types())
+    types = np.unique(_get_channel_types(info))
 
     if requested_conditions is None:
         requested_conditions = design_matrix.columns
@@ -98,7 +116,7 @@ def plot_glm_topo(raw, glm_estimates, design_matrix,
 
     for t_idx, t in enumerate(types):
 
-        estmrg, pos, chs, sphere = _handle_overlaps(raw, t, sphere, estimates)
+        estmrg, pos, chs, sphere = _handle_overlaps(info, t, sphere, estimates)
 
         for idx, label in enumerate(design_matrix.columns):
             if label in requested_conditions:
@@ -135,14 +153,14 @@ def plot_glm_topo(raw, glm_estimates, design_matrix,
     return _get_fig_from_axes(axes)
 
 
-def plot_glm_contrast_topo(raw, contrast, figsize=(12, 7), sphere=None):
+def plot_glm_contrast_topo(inst, contrast, figsize=(12, 7), sphere=None):
     """
     Plot topomap of NIRS GLM data.
 
     Parameters
     ----------
-    raw : instance of Raw
-        Haemoglobin data.
+    inst : instance of Info or Raw
+        Raw data or info structure used to generate the GLM results.
     contrast : dict
         As in nilearn.stats.compute_contrast.
     figsize : numbers
@@ -156,9 +174,20 @@ def plot_glm_contrast_topo(raw, contrast, figsize=(12, 7), sphere=None):
         Figure of each design matrix componenent for hbo (top row)
         and hbr (bottom row).
     """
+    warn('"plot_glm_contrast_topo" has been deprecated in favor of the more '
+         'comprehensive GLM class and will be removed in v1.0.0. '
+         'Use the ContrastResults class "plot_topo()" method instead.',
+         DeprecationWarning)
+    return _plot_glm_contrast_topo(inst, contrast, figsize=figsize,
+                                   sphere=sphere)
+
+
+def _plot_glm_contrast_topo(inst, contrast, figsize=(12, 7), sphere=None):
+
+    info = deepcopy(inst if isinstance(inst, Info) else inst.info)
 
     # Extract types. One subplot is created per type (hbo/hbr)
-    types = np.unique(raw.get_channel_types())
+    types = np.unique(_get_channel_types(info))
 
     # Extract values to plot and rescale to uM
     estimates = contrast.effect[0]
@@ -176,7 +205,7 @@ def plot_glm_contrast_topo(raw, contrast, figsize=(12, 7), sphere=None):
 
     for t_idx, t in enumerate(types):
 
-        estmrg, pos, chs, sphere = _handle_overlaps(raw, t, sphere, estimates)
+        estmrg, pos, chs, sphere = _handle_overlaps(info, t, sphere, estimates)
 
         # Deal with case when only a single chroma is available
         if len(types) == 1:
@@ -212,7 +241,7 @@ def plot_glm_contrast_topo(raw, contrast, figsize=(12, 7), sphere=None):
     return fig
 
 
-def plot_glm_group_topo(raw, statsmodel_df,
+def plot_glm_group_topo(inst, statsmodel_df,
                         value="Coef.",
                         axes=None,
                         threshold=False,
@@ -231,8 +260,8 @@ def plot_glm_group_topo(raw, statsmodel_df,
 
     Parameters
     ----------
-    raw : instance of Raw
-        Haemoglobin data.
+    inst : instance of Info or Raw
+        Raw data or info structure used to generate the GLM results.
     statsmodel_df : DataFrame
         Dataframe created from a statsmodel summary.
     value : String
@@ -274,18 +303,19 @@ def plot_glm_group_topo(raw, statsmodel_df,
     fig : figure
         Figure with topographic representation of statsmodel_df value.
     """
+    info = deepcopy(inst if isinstance(inst, Info) else inst.info)
 
     # Check that the channels in two inputs match
-    if not (raw.ch_names == list(statsmodel_df["ch_name"].values)):
-        if len(raw.ch_names) < len(list(statsmodel_df["ch_name"].values)):
+    if not (info.ch_names == list(statsmodel_df["ch_name"].values)):
+        if len(info.ch_names) < len(list(statsmodel_df["ch_name"].values)):
             print("Reducing GLM results to match MNE data")
-            statsmodel_df["Keep"] = [g in raw.ch_names
+            statsmodel_df["Keep"] = [g in info.ch_names
                                      for g in statsmodel_df["ch_name"]]
             statsmodel_df = statsmodel_df.query("Keep == True")
         else:
             warn("MNE data structure does not match regression results")
     statsmodel_df = statsmodel_df.set_index('ch_name')
-    statsmodel_df = statsmodel_df.reindex(raw.ch_names)
+    statsmodel_df = statsmodel_df.reindex(info.ch_names)
 
     # Extract estimate of interest to plot
     estimates = statsmodel_df[value].values
@@ -321,7 +351,7 @@ def plot_glm_group_topo(raw, statsmodel_df,
         cmap = mpl.cm.RdBu_r
     norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
 
-    estmrg, pos, chs, sphere = _handle_overlaps(raw, t, sphere, estimates)
+    estmrg, pos, chs, sphere = _handle_overlaps(info, t, sphere, estimates)
 
     mne.viz.topomap.plot_topomap(estmrg, pos,
                                  extrapolate=extrapolate,
@@ -348,12 +378,13 @@ def plot_glm_group_topo(raw, statsmodel_df,
     return axes
 
 
-def _handle_overlaps(raw, t, sphere, estimates):
+def _handle_overlaps(info, t, sphere, estimates):
     """Prepare for topomap including merging channels"""
-    picks = _picks_to_idx(raw.info, t, exclude=[], allow_empty=True)
-    raw_subset = raw.copy().pick(picks=picks)
+    picks = _picks_to_idx(info, t, exclude=[], allow_empty=True)
+    pick_names = [info.ch_names[p] for p in picks]
+    info_subset = info.copy().pick_channels(pick_names)
     _, pos, merge_channels, ch_names, ch_type, sphere, clip_origin = \
-        mne.viz.topomap._prepare_topomap_plot(raw_subset, t, sphere=sphere)
+        mne.viz.topomap._prepare_topomap_plot(info_subset, t, sphere=sphere)
     estmrg, ch_names = _merge_ch_data(estimates.copy()[picks], t, ch_names)
     return estmrg, pos, ch_names, sphere
 


Test patch:
diff --git a/mne_nirs/statistics/tests/test_glm_type.py b/mne_nirs/statistics/tests/test_glm_type.py
new file mode 100644
index 000000000..f5a2640a3
--- /dev/null
+++ b/mne_nirs/statistics/tests/test_glm_type.py
@@ -0,0 +1,163 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import os
+
+import pandas
+import pytest
+import numpy as np
+import matplotlib
+from matplotlib.pyplot import Axes
+
+import mne
+import nilearn
+
+from mne_nirs.statistics import RegressionResults
+from mne_nirs.experimental_design import make_first_level_design_matrix
+from mne_nirs.statistics import run_glm
+
+
+def _get_minimal_haemo_data(tmin=0, tmax=60):
+    raw = mne.io.read_raw_nirx(os.path.join(
+        mne.datasets.fnirs_motor.data_path(), 'Participant-1'), preload=False)
+    raw.crop(tmax=tmax, tmin=tmin)
+    raw = mne.preprocessing.nirs.optical_density(raw)
+    raw = mne.preprocessing.nirs.beer_lambert_law(raw)
+    return raw
+
+
+def _get_glm_result(tmax=60, tmin=0):
+    raw = _get_minimal_haemo_data(tmin=tmin, tmax=tmax)
+    design_matrix = make_first_level_design_matrix(raw, stim_dur=5.,
+                                                   drift_order=1,
+                                                   drift_model='polynomial')
+    return run_glm(raw, design_matrix)
+
+
+def _get_glm_contrast_result(tmin=60, tmax=400):
+    raw = _get_minimal_haemo_data(tmin=tmin, tmax=tmax)
+    design_matrix = make_first_level_design_matrix(raw, stim_dur=5.,
+                                                   drift_order=1,
+                                                   drift_model='polynomial')
+
+    glm_est = run_glm(raw, design_matrix)
+
+    contrast_matrix = np.eye(design_matrix.shape[1])
+    basic_conts = dict([(column, contrast_matrix[i])
+                        for i, column in enumerate(design_matrix.columns)])
+    contrast_LvR = basic_conts['1.0'] - basic_conts['2.0']
+
+    return glm_est.compute_contrast(contrast_LvR)
+
+
+def test_create_results_glm():
+
+    # Create a relevant info structure
+    raw = _get_minimal_haemo_data()
+
+    # Create a minimal structure
+    res = _get_glm_result()
+
+    # Get arguments for type so we can test errors below
+    info = raw.info
+    minimal_structure = res._data
+
+    # Test construction
+    with pytest.raises(TypeError, match='must be a dictionary'):
+        _ = RegressionResults(info, np.zeros((5, 2)), 1)
+    with pytest.raises(TypeError, match='must be a dictionary'):
+        _ = RegressionResults(info, 3.2, 1)
+    with pytest.raises(TypeError, match='must be a dictionary'):
+        _ = RegressionResults(info, [], 1)
+    with pytest.raises(TypeError, match='must be a dictionary'):
+        _ = RegressionResults(info, "hbo", 1)
+    with pytest.raises(TypeError, match='keys must match'):
+        _ = RegressionResults(info, _take(4, minimal_structure), 1)
+
+    onewrongname = _take(55, minimal_structure)
+    onewrongname["test"] = onewrongname["S1_D1 hbo"]
+    with pytest.raises(TypeError, match='must match ch_names'):
+        _ = RegressionResults(info, onewrongname, 1)
+
+    # Test properties
+    assert len(res) == len(raw.ch_names)
+
+
+def test_results_glm_properties():
+
+    n_channels = 56
+
+    res = _get_glm_result()
+
+    # Test ContainsMixin
+    assert 'hbo' in res
+    assert 'hbr' in res
+    assert 'meg' not in res
+
+    # Test copy
+    assert len(res) == len(res.copy())
+
+    # Test picks
+    assert len(res.copy().pick(picks=range(4))) == 4
+    assert len(res.copy().pick(picks="S1_D1 hbr")) == 1
+    assert len(res.copy().pick(picks=["S1_D1 hbr"])) == 1
+    assert len(res.copy().pick(picks=["S1_D1 hbr", "S1_D1 hbo"])) == 2
+    assert len(res.copy().pick(picks=["S1_D1 hbr", "S1_D1 XXX"])) == 1
+    assert len(res.copy().pick(picks=["S1_D1 hbr", "S1_D1 hbr"])) == 1
+    assert len(res.copy().pick(picks="fnirs")) == n_channels
+    assert len(res.copy().pick(picks="hbo")) == n_channels / 2
+    assert len(res.copy().pick(picks="hbr")) == n_channels / 2
+
+    # Test results
+    assert len(res.theta()) == n_channels
+    assert len(res.copy().pick(picks=range(4)).theta()) == 4
+    assert len(res.copy().pick(picks=3).theta()) == 1
+    assert res.copy().pick(picks=3).theta()[0].shape == (3, 1)
+
+    # Test models
+    assert len(res.model()) == n_channels
+    assert len(res.copy().pick(picks=range(8)).model()) == 8
+    assert type(res.model()[0]) is nilearn.glm.regression.ARModel
+
+    assert isinstance(res.to_dataframe(), pandas.DataFrame)
+
+
+def test_glm_scatter():
+
+    assert isinstance(_get_glm_result().scatter(), Axes)
+    assert isinstance(_get_glm_contrast_result().scatter(), Axes)
+    _get_glm_result(tmax=2974, tmin=0).surface_projection(condition="3.0",
+                                                          view="dorsal")
+
+
+def test_results_glm_export_dataframe():
+
+    n_channels = 56
+    res = _get_glm_result(tmax=400)
+    df = res.to_dataframe()
+
+    assert df.shape == (6 * n_channels, 12)
+
+
+def test_create_results_glm_contrast():
+
+    # Create a minimal structure
+    res = _get_glm_contrast_result()
+    assert isinstance(res._data, nilearn.glm.contrasts.Contrast)
+    assert isinstance(res.info, mne.Info)
+
+    # Test copy
+    assert len(res) == len(res.copy())
+
+    assert isinstance(res.plot_topo(), matplotlib.figure.Figure)
+
+    n_channels = 56
+    assert isinstance(res.to_dataframe(), pandas.DataFrame)
+    df = res.to_dataframe()
+    assert df.shape == (n_channels, 10)
+
+
+def _take(n, mydict):
+    """Return first n items of the iterable as a list"""
+    return {k: mydict[k] for k in list(mydict)[:n]}
diff --git a/mne_nirs/statistics/tests/test_statistics.py b/mne_nirs/statistics/tests/test_statistics.py
index c280461c6..77b87a157 100644
--- a/mne_nirs/statistics/tests/test_statistics.py
+++ b/mne_nirs/statistics/tests/test_statistics.py
@@ -10,7 +10,7 @@
 from mne.simulation import add_noise
 
 from mne_nirs.experimental_design import make_first_level_design_matrix
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm, run_GLM
 from mne_nirs.simulation import simulate_nirs_raw
 
 iir_filter = [1., -0.58853134, -0.29575669, -0.52246482, 0.38735476, 0.024286]
@@ -22,17 +22,23 @@ def test_run_GLM():
     design_matrix = make_first_level_design_matrix(raw, stim_dur=5.,
                                                    drift_order=1,
                                                    drift_model='polynomial')
-    glm_estimates = run_GLM(raw, design_matrix)
+    glm_estimates = run_glm(raw, design_matrix)
+
+    # Test backwards compatibility
+    old_res = run_GLM(raw, design_matrix)
+    assert old_res.keys() == glm_estimates.data.keys()
+    assert (old_res["Simulated"].theta ==
+            glm_estimates.data["Simulated"].theta).all()
 
     assert len(glm_estimates) == len(raw.ch_names)
 
     # Check the estimate is correct within 10% error
-    assert abs(glm_estimates["Simulated"].theta[0] - 1.e-6) < 0.1e-6
+    assert abs(glm_estimates.pick("Simulated").theta()[0][0] - 1.e-6) < 0.1e-6
 
     # ensure we return the same type as nilearn to encourage compatibility
     _, ni_est = nilearn.glm.first_level.run_glm(
         raw.get_data(0).T, design_matrix.values)
-    assert type(ni_est) == type(glm_estimates)
+    assert isinstance(glm_estimates._data, type(ni_est))
 
 
 def test_run_GLM_order():
@@ -42,22 +48,22 @@ def test_run_GLM_order():
                                                    drift_model='polynomial')
 
     # Default should be first order AR
-    glm_estimates = run_GLM(raw, design_matrix)
-    assert glm_estimates['Simulated'].model.order == 1
+    glm_estimates = run_glm(raw, design_matrix)
+    assert glm_estimates.pick("Simulated").model()[0].order == 1
 
     # Default should be first order AR
-    glm_estimates = run_GLM(raw, design_matrix, noise_model='ar2')
-    assert glm_estimates['Simulated'].model.order == 2
+    glm_estimates = run_glm(raw, design_matrix, noise_model='ar2')
+    assert glm_estimates.pick("Simulated").model()[0].order == 2
 
-    glm_estimates = run_GLM(raw, design_matrix, noise_model='ar7')
-    assert glm_estimates['Simulated'].model.order == 7
+    glm_estimates = run_glm(raw, design_matrix, noise_model='ar7')
+    assert glm_estimates.pick("Simulated").model()[0].order == 7
 
     # Auto should be 4 times sample rate
     cov = Covariance(np.ones(1) * 1e-11, raw.ch_names,
                      raw.info['bads'], raw.info['projs'], nfree=0)
     raw = add_noise(raw, cov, iir_filter=iir_filter)
-    glm_estimates = run_GLM(raw, design_matrix, noise_model='auto')
-    assert glm_estimates['Simulated'].model.order == 3 * 4
+    glm_estimates = run_glm(raw, design_matrix, noise_model='auto')
+    assert glm_estimates.pick("Simulated").model()[0].order == 3 * 4
 
     raw = simulate_nirs_raw(sig_dur=10, stim_dur=5., sfreq=2)
     cov = Covariance(np.ones(1) * 1e-11, raw.ch_names,
@@ -67,5 +73,5 @@ def test_run_GLM_order():
                                                    drift_order=1,
                                                    drift_model='polynomial')
     # Auto should be 4 times sample rate
-    glm_estimates = run_GLM(raw, design_matrix, noise_model='auto')
-    assert glm_estimates['Simulated'].model.order == 2 * 4
+    glm_estimates = run_glm(raw, design_matrix, noise_model='auto')
+    assert glm_estimates.pick("Simulated").model()[0].order == 2 * 4
diff --git a/mne_nirs/statistics/tests/test_statsmodels.py b/mne_nirs/statistics/tests/test_statsmodels.py
index bbed3a933..16d73e6ad 100644
--- a/mne_nirs/statistics/tests/test_statsmodels.py
+++ b/mne_nirs/statistics/tests/test_statsmodels.py
@@ -11,8 +11,7 @@
 
 from ...simulation import simulate_nirs_raw
 from ...experimental_design import make_first_level_design_matrix
-from ...statistics import run_GLM, statsmodels_to_results
-from ...utils._io import glm_to_tidy
+from ...statistics import run_glm, statsmodels_to_results
 
 
 @pytest.mark.parametrize('func', ('mixedlm', 'ols', 'rlm'))
@@ -33,9 +32,9 @@ def test_statsmodel_to_df(func):
                                 isi_min=15., isi_max=45.)
         raw._data += np.random.normal(0, np.sqrt(1e-12), raw._data.shape)
         design_matrix = make_first_level_design_matrix(raw, stim_dur=5.0)
-        glm_est = run_GLM(raw, design_matrix)
+        glm_est = run_glm(raw, design_matrix)
         with pytest.warns(RuntimeWarning, match='Non standard source detect'):
-            cha = glm_to_tidy(raw, glm_est, design_matrix)
+            cha = glm_est.to_dataframe()
         cha["ID"] = '%02d' % n
         df_cha = df_cha.append(cha)
     df_cha["theta"] = df_cha["theta"] * 1.0e6
diff --git a/mne_nirs/utils/tests/test_io.py b/mne_nirs/utils/tests/test_io.py
index e3b25956c..374bfc808 100644
--- a/mne_nirs/utils/tests/test_io.py
+++ b/mne_nirs/utils/tests/test_io.py
@@ -11,7 +11,7 @@
 
 from mne_nirs.experimental_design import make_first_level_design_matrix
 from mne_nirs.utils._io import glm_to_tidy, _tidy_long_to_wide
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm
 
 
 @pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
@@ -30,8 +30,8 @@ def test_io():
                                                    stim_dur=5.0,
                                                    drift_order=3,
                                                    drift_model='polynomial')
-    glm_est = run_GLM(raw_haemo, design_matrix)
-    df = glm_to_tidy(raw_haemo, glm_est, design_matrix)
+    glm_est = run_glm(raw_haemo, design_matrix)
+    df = glm_to_tidy(raw_haemo, glm_est.data, design_matrix)
     assert df.shape == (48, 12)
     assert set(df.columns) == {'ch_name', 'Condition', 'df', 'mse', 'p_value',
                                't', 'theta', 'Source', 'Detector', 'Chroma',
@@ -54,14 +54,14 @@ def test_io():
                         for i, column in enumerate(design_matrix.columns)])
     contrast_LvR = basic_conts['2.0'] - basic_conts['3.0']
 
-    contrast = mne_nirs.statistics.compute_contrast(glm_est, contrast_LvR)
+    contrast = mne_nirs.statistics.compute_contrast(glm_est.data, contrast_LvR)
     df = glm_to_tidy(raw_haemo, contrast, design_matrix)
     assert df.shape == (6, 10)
     assert set(df.columns) == {'ch_name', 'ContrastType', 'z_score', 'stat',
                                'p_value', 'effect', 'Source', 'Detector',
                                'Chroma', 'Significant'}
 
-    contrast = mne_nirs.statistics.compute_contrast(glm_est, contrast_LvR,
+    contrast = mne_nirs.statistics.compute_contrast(glm_est.data, contrast_LvR,
                                                     contrast_type='F')
     df = glm_to_tidy(raw_haemo, contrast, design_matrix, wide=False)
     df = _tidy_long_to_wide(df)
diff --git a/mne_nirs/visualisation/tests/test_visualisation.py b/mne_nirs/visualisation/tests/test_visualisation.py
index ef03c12e8..bd9ab4e43 100644
--- a/mne_nirs/visualisation/tests/test_visualisation.py
+++ b/mne_nirs/visualisation/tests/test_visualisation.py
@@ -10,7 +10,7 @@
 from mne_nirs.experimental_design.tests.test_experimental_design import \
     _load_dataset
 from mne_nirs.experimental_design import make_first_level_design_matrix
-from mne_nirs.statistics import run_GLM
+from mne_nirs.statistics import run_glm
 from mne_nirs.visualisation import plot_glm_topo, plot_glm_surface_projection
 from mne_nirs.utils import glm_to_tidy
 
@@ -74,44 +74,44 @@ def test_run_plot_GLM_topo():
                                                    drift_model='polynomial')
     raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
     raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
-    glm_estimates = run_GLM(raw_haemo, design_matrix)
-    fig = plot_glm_topo(raw_haemo, glm_estimates, design_matrix)
+    glm_estimates = run_glm(raw_haemo, design_matrix)
+    fig = plot_glm_topo(raw_haemo, glm_estimates.data, design_matrix)
     # 5 conditions (A,B,C,Drift,Constant) * two chroma + 2xcolorbar
     assert len(fig.axes) == 12
 
     # Two conditions * two chroma + 2 x colorbar
-    fig = plot_glm_topo(raw_haemo, glm_estimates, design_matrix,
+    fig = plot_glm_topo(raw_haemo, glm_estimates.data, design_matrix,
                         requested_conditions=['A', 'B'])
     assert len(fig.axes) == 6
 
     # Two conditions * one chroma + 1 x colorbar
     with pytest.warns(RuntimeWarning, match='Reducing GLM results'):
         fig = plot_glm_topo(raw_haemo.copy().pick(picks="hbo"),
-                            glm_estimates, design_matrix,
+                            glm_estimates.data, design_matrix,
                             requested_conditions=['A', 'B'])
     assert len(fig.axes) == 3
 
     # One conditions * two chroma + 2 x colorbar
-    fig = plot_glm_topo(raw_haemo, glm_estimates, design_matrix,
+    fig = plot_glm_topo(raw_haemo, glm_estimates.data, design_matrix,
                         requested_conditions=['A'])
     assert len(fig.axes) == 4
 
     # One conditions * one chroma + 1 x colorbar
     with pytest.warns(RuntimeWarning, match='Reducing GLM results'):
         fig = plot_glm_topo(raw_haemo.copy().pick(picks="hbo"),
-                            glm_estimates,
+                            glm_estimates.data,
                             design_matrix, requested_conditions=['A'])
     assert len(fig.axes) == 2
 
     # One conditions * one chroma + 0 x colorbar
     with pytest.warns(RuntimeWarning, match='Reducing GLM results'):
         fig = plot_glm_topo(raw_haemo.copy().pick(picks="hbo"),
-                            glm_estimates, design_matrix,
+                            glm_estimates.data, design_matrix,
                             colorbar=False, requested_conditions=['A'])
     assert len(fig.axes) == 1
 
     # Ensure warning thrown if glm estimates is missing channels from raw
-    glm_estimates_subset = {a: glm_estimates[a]
+    glm_estimates_subset = {a: glm_estimates.data[a]
                             for a in raw_haemo.ch_names[0:3]}
     with pytest.raises(RuntimeError, match="does not match regression"):
         plot_glm_topo(raw_haemo, glm_estimates_subset, design_matrix)
@@ -127,12 +127,12 @@ def test_run_plot_GLM_contrast_topo():
                                                    drift_model='polynomial')
     raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
     raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
-    glm_est = run_GLM(raw_haemo, design_matrix)
+    glm_est = run_glm(raw_haemo, design_matrix)
     contrast_matrix = np.eye(design_matrix.shape[1])
     basic_conts = dict([(column, contrast_matrix[i])
                         for i, column in enumerate(design_matrix.columns)])
     contrast_LvR = basic_conts['A'] - basic_conts['B']
-    contrast = mne_nirs.statistics.compute_contrast(glm_est, contrast_LvR)
+    contrast = mne_nirs.statistics.compute_contrast(glm_est.data, contrast_LvR)
     fig = mne_nirs.visualisation.plot_glm_contrast_topo(raw_haemo, contrast)
     assert len(fig.axes) == 3
 
@@ -148,12 +148,12 @@ def test_run_plot_GLM_contrast_topo_single_chroma():
     raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
     raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
     raw_haemo = raw_haemo.pick(picks='hbo')
-    glm_est = run_GLM(raw_haemo, design_matrix)
+    glm_est = run_glm(raw_haemo, design_matrix)
     contrast_matrix = np.eye(design_matrix.shape[1])
     basic_conts = dict([(column, contrast_matrix[i])
                         for i, column in enumerate(design_matrix.columns)])
     contrast_LvR = basic_conts['A'] - basic_conts['B']
-    contrast = mne_nirs.statistics.compute_contrast(glm_est, contrast_LvR)
+    contrast = mne_nirs.statistics.compute_contrast(glm_est.data, contrast_LvR)
     fig = mne_nirs.visualisation.plot_glm_contrast_topo(raw_haemo, contrast)
     assert len(fig.axes) == 2
 
@@ -177,8 +177,8 @@ def test_run_plot_GLM_projection():
                                                    drift_model='polynomial')
     raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
     raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
-    glm_estimates = run_GLM(raw_haemo, design_matrix)
-    df = glm_to_tidy(raw_haemo, glm_estimates, design_matrix)
+    glm_estimates = run_glm(raw_haemo, design_matrix)
+    df = glm_to_tidy(raw_haemo, glm_estimates.data, design_matrix)
     df = df.query("Chroma in 'hbo'")
     df = df.query("Condition in 'A'")