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

Issue:
Common MNE annotation names do not work with nilearn
#### Whats the problem?

Nilearn recently changed the requirements for design matrix column names in https://github.com/nilearn/nilearn/pull/3025 which disallowed strings containing only numbers and strings with `/` as column names. 

Getting rid of the numbers is no problem, I prefer meaningful annotation names anyway.

However, I think its common in MNE to use names like `tapping/left` or `response/right`. See https://mne.tools/dev/auto_tutorials/epochs/40_autogenerate_metadata.html?highlight=hed#keeping-only-the-first-events-of-a-group

But nilearn now requires that the column names pass the following test:

```python
'response/right'.isidentifier()
# False

'response_right'.isidentifier()
# True
```

See lines: https://github.com/nilearn/nilearn/blob/ab9c68e502c4b92695cd658f5839a8207fd5e07e/nilearn/glm/first_level/hemodynamic_models.py#L428-L436

#### Solution
See #403 for an initial attempt at a fix. Other suggestions are welcome
Automatically sanitse annotation names that arent compatible with nilearn
#### Whats the problem?

Nilearn recently changed the requirements for design matrix column names in https://github.com/nilearn/nilearn/pull/3025 which disallowed strings containing only numbers and strings with `/` as column names. 

Getting rid of the numbers is no problem, I prefer meaningful annotation names anyway.

However, I think its common in MNE to use names like `tapping/left` or `response/right`. See https://mne.tools/dev/auto_tutorials/epochs/40_autogenerate_metadata.html?highlight=hed#keeping-only-the-first-events-of-a-group

But nilearn now requires that the column names pass the following test:

```python
'response/right'.isidentifier()
# False

'response_right'.isidentifier()
# True
```

See lines: https://github.com/nilearn/nilearn/blob/ab9c68e502c4b92695cd658f5839a8207fd5e07e/nilearn/glm/first_level/hemodynamic_models.py#L428-L436

#### Whats the fix?

This PR checks that our design matrix column names are compatible with nilearn, and if they are not then prepends a `t_` to force them to be valid.

#### The problem

This PR also changes `/` to `_` in each design matrix column name, but this is not then harmonious with the rest of MNE. Further, it means that in some places of my script I use `tapping/left` and other places I will need to use `tapping_left`. This solution doesn't seem great, what is a better alternative?


Discussion / hints:

@agramfort and @larsoner I am a bit stumped about what to do here. Nilearn has removed the ability to use `/` in design matrix columns. These columns represent experimental conditions and so I have previously just used the annotation name. However, annotation names commonly include `/` in them. Any fresh ideas of how I can tackle this problem?
I will pin nilearn to previous version rather than main until I know how to fix this.
To me this seems like a bug / regression in nilearn. It's entirely possible they haven't considered spaces or `/` or number-first names to be widespread use cases, but clearly they are at least for us. I would open an issue and see if it's possible for them to fix -- to me it seems like they should sanitize the strings when doing whatever processing they need to do, and then restore them afterward...

If they don't like this solution, it's what I'd probably implement at the MNE-NIRS end.
We now know that we can not use the `/` symbol in the design matrix. So I will try and write some sort of sanitisation approach as suggested by larsoner and will use `___` to equal `/`.

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 4e1483c4c..ae71d2ebc 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -10,8 +10,10 @@ _xvfb: &xvfb
 
 jobs:
     build_docs:
-      machine:
-        image: ubuntu-2004:202111-01
+      docker:
+        # Use 18.04 rather than 20.04 because MESA 20.0.8 on 18.04 has working
+        # transparency but 21.0.3 on 20.04 does not!
+        - image: cimg/base:stable-18.04
       steps:
         - checkout
         - run:
@@ -32,14 +34,21 @@ jobs:
             name: Set BASH_ENV
             command: |
               set -e
-              python3 -m venv ~/python_env --upgrade-deps
+              sudo apt update -qq
+              sudo apt install -qq libosmesa6 libglx-mesa0 libopengl0 libglx0 libdbus-1-3 \
+                libxkbcommon-x11-0 libxcb-icccm4 libxcb-image0 libxcb-keysyms1 libxcb-randr0 \
+                libxcb-render-util0 libxcb-shape0 libxcb-xfixes0 libxcb-xinerama0 \
+                graphviz optipng \
+                python3.8-venv python3-venv \
+                xvfb libxft2 ffmpeg
+              python3.8 -m venv ~/python_env
               echo "set -e" >> $BASH_ENV
-              echo "export DISPLAY=:99" >> $BASH_ENV
               echo "export OPENBLAS_NUM_THREADS=4" >> $BASH_ENV
               echo "export XDG_RUNTIME_DIR=/tmp/runtime-circleci" >> $BASH_ENV
-              source tools/get_minimal_commands.sh
+              echo "export MNE_FULL_DATE=true" >> $BASH_ENV
               echo "export MNE_3D_OPTION_ANTIALIAS=false" >> $BASH_ENV
-              echo "export MNE_3D_BACKEND=pyvista" >> $BASH_ENV
+              source tools/get_minimal_commands.sh
+              echo "export MNE_3D_BACKEND=pyvistaqt" >> $BASH_ENV
               echo "export PATH=~/.local/bin/:$PATH" >> $BASH_ENV
               echo "source ~/python_env/bin/activate" >> $BASH_ENV
               mkdir -p ~/.local/bin
@@ -47,16 +56,11 @@ jobs:
               echo "BASH_ENV:"
               cat $BASH_ENV
               mkdir -p ~/mne_data
-              touch pattern.txt;
-
+              touch pattern.txt
         - run:
-            name: Install 3D rendering libraries \ PyQt5 dependencies \ graphviz \ optipng (for optimized images)
+            name: check neuromag2ft
             command: |
-              sudo apt update
-              sudo apt install libosmesa6 libglx-mesa0 libopengl0 libglx0 libdbus-1-3 \
-                  libxkbcommon-x11-0 libxcb-icccm4 libxcb-image0 libxcb-keysyms1 libxcb-randr0 libxcb-render-util0 libxcb-shape0 libxcb-xfixes0 libxcb-xinerama0 \
-                  graphviz \
-                  optipng
+              neuromag2ft --version
 
         - run:
             <<: *xvfb
@@ -83,7 +87,7 @@ jobs:
         - run:
             name: Get Python running
             command: |
-              python -m pip install --upgrade --progress-bar off pip setuptools
+              python -m pip install --upgrade pip setuptools wheel
               python -m pip install --upgrade --progress-bar off --pre sphinx
               python -m pip install --upgrade --progress-bar off -r requirements.txt -r requirements_testing.txt -r requirements_doc.txt
               python -m pip install -e .
diff --git a/.github/workflows/linux_pip.yml b/.github/workflows/linux_pip.yml
index 01722bd98..5ccd93293 100644
--- a/.github/workflows/linux_pip.yml
+++ b/.github/workflows/linux_pip.yml
@@ -15,7 +15,7 @@ concurrency:
 jobs:
   # PIP + non-default stim channel + log level info
   job:
-    name: 'py3.8'
+    name: 'py3.10'
     runs-on: ubuntu-20.04
     if: "!contains(github.event.head_commit.message, '[skip tests]')"
     defaults:
diff --git a/environment.yml b/environment.yml
index be54413bb..44f2b5b6a 100644
--- a/environment.yml
+++ b/environment.yml
@@ -7,12 +7,17 @@ dependencies:
 - numpy
 - scipy
 - matplotlib
+- tqdm
+- pooch>=1.5
+- decorator
+- h5io
+- packaging
 - numba
-- pandas==1.3.2
+- pandas
 - xlrd
 - scikit-learn
 - h5py
-- h5io
+- jinja2
 - pillow
 - statsmodels
 - jupyter
@@ -20,25 +25,23 @@ dependencies:
 - psutil
 - numexpr
 - imageio
-- tqdm
 - spyder-kernels>=1.10.0
 - imageio-ffmpeg>=0.4.1
 - vtk>=9.0.1
-- pyvista>=0.30
+- traitlets
+- pyvista>=0.32
 - pyvistaqt>=0.4
 - qdarkstyle
 - darkdetect
 - dipy
 - nibabel
+- nilearn
 - python-picard
 - pyqt!=5.15.3
 - mffpy>=0.5.7
 - ipywidgets
-- lxml
-- pytables
-- pooch
-- nilearn
 - ipyvtklink
 - mne-qt-browser
+- pymatreader
 - pip:
   - https://github.com/mne-tools/mne-python/archive/main.zip
diff --git a/examples/general/plot_11_hrf_measured.py b/examples/general/plot_11_hrf_measured.py
index d369d676d..da0bd1167 100644
--- a/examples/general/plot_11_hrf_measured.py
+++ b/examples/general/plot_11_hrf_measured.py
@@ -7,10 +7,12 @@
 In this example we analyse data from a real multichannel
 functional near-infrared spectroscopy (fNIRS)
 experiment (see :ref:`tut-fnirs-hrf-sim` for a simplified simulated
-analysis). The experiment consists of three conditions
-1) tapping with the left hand,
-2) tapping with the right hand,
-3) a control condition where the participant does nothing.
+analysis). The experiment consists of three conditions:
+
+1. tapping with the left hand,
+2. tapping with the right hand,
+3. a control condition where the participant does nothing.
+
 We use a GLM analysis to examine the neural activity associated with
 the different tapping conditions.
 An alternative epoching style analysis on the same data can be
@@ -27,12 +29,12 @@
    :local:
    :depth: 2
 
-.. note:: Parts of this tutorial require the latest development version of MNE-Python. See these instructions for  
+.. note:: Parts of this tutorial require the latest development version of MNE-Python. See these instructions for
           `how to upgrade <https://mne.tools/dev/install/updating.html>`__.
           But basically boils down to running
           ``pip install -U --no-deps https://github.com/mne-tools/mne-python/archive/main.zip``.
           Sections of the code that require this version will be noted below.
-          
+
 """
 # sphinx_gallery_thumbnail_number = 9
 
@@ -98,10 +100,12 @@
 # Next we update the annotations by assigning names to each trigger ID.
 # Then we crop the recording to the section containing our
 # experimental conditions.
-
+#
+# Because of limitations with ``nilearn``, we use ``'_'`` to separate conditions
+# rather than the standard ``'/'``.
 raw_intensity.annotations.rename({'1.0': 'Control',
-                                  '2.0': 'Tapping/Left',
-                                  '3.0': 'Tapping/Right'})
+                                  '2.0': 'Tapping_Left',
+                                  '3.0': 'Tapping_Right'})
 raw_intensity.annotations.delete(raw_intensity.annotations.description == '15.0')
 raw_intensity.annotations.set_durations(5)
 
@@ -166,7 +170,7 @@
 # .. sidebar:: Relevant literature
 #
 #    For further discussion on design matrices see
-#    the Nilearn examples. Specifically the 
+#    the Nilearn examples. Specifically the
 #    `first level model example <http://nilearn.github.io/auto_examples/04_glm_first_level/plot_first_level_details.html>`_.
 #
 # Next we create a model to fit our data to.
@@ -221,7 +225,7 @@
 # Examine expected response
 # -------------------------
 #
-# The matrices above can be a bit abstract as they encompase multiple 
+# The matrices above can be a bit abstract as they encompase multiple
 # conditions and regressors.
 # Instead we can examine a single condition.
 # Here we observe the boxcar function for a single condition,
@@ -239,7 +243,7 @@
 
 s = mne_nirs.experimental_design.create_boxcar(raw_intensity, stim_dur=5.0)
 plt.plot(raw_intensity.times, s[:, 1])
-plt.plot(design_matrix['Tapping/Left'])
+plt.plot(design_matrix['Tapping_Left'])
 plt.xlim(180, 300)
 plt.legend(["Stimulus", "Expected Response"])
 plt.xlabel("Time (s)")
@@ -336,7 +340,7 @@
 # negative of HbO as expected.
 
 glm_est = run_glm(raw_haemo, design_matrix)
-glm_est.plot_topo(conditions=['Tapping/Left', 'Tapping/Right'])
+glm_est.plot_topo(conditions=['Tapping_Left', 'Tapping_Right'])
 
 
 # %%
@@ -364,7 +368,7 @@
 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']
+conditions = ['Tapping_Right']
 
 glm_hbo.plot_topo(axes=axes[0], colorbar=False, conditions=conditions)
 
@@ -380,7 +384,7 @@
 # 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")
+glm_est.copy().surface_projection(condition="Tapping_Right", view="dorsal", chroma="hbo")
 
 
 # %%
@@ -412,7 +416,7 @@
 groups = dict(Left_ROI=picks_pair_to_idx(raw_haemo, left),
               Right_ROI=picks_pair_to_idx(raw_haemo, right))
 
-conditions = ['Control', 'Tapping/Left', 'Tapping/Right']
+conditions = ['Control', 'Tapping_Left', 'Tapping_Right']
 
 df = glm_est.to_dataframe_region_of_interest(groups, conditions)
 
@@ -439,7 +443,7 @@
 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_LvR = basic_conts['Tapping_Left'] - basic_conts['Tapping_Right']
 
 contrast = glm_est.compute_contrast(contrast_LvR)
 contrast.plot_topo()
@@ -474,8 +478,8 @@
 # the tapping, but we do expect 5% or less for the false positive rate.
 
 (df
- .query('Condition in ["Control", "Tapping/Left", "Tapping/Right"]')
- .groupby(['Condition', 'Chroma'])
+ .query('Condition in ["Control", "Tapping_Left", "Tapping_Right"]')
+ .drop(['df', 'mse', 'p_value', 't'], axis=1)
+ .groupby(['Condition', 'Chroma', 'ch_name'])
  .agg(['mean'])
- .drop(['df', 'mse', 'p_value', 't'], 1)
  )
diff --git a/examples/general/plot_12_group_glm.py b/examples/general/plot_12_group_glm.py
index 903953eea..35cad5a85 100644
--- a/examples/general/plot_12_group_glm.py
+++ b/examples/general/plot_12_group_glm.py
@@ -172,11 +172,23 @@
 # Here we also resample to a 0.3 Hz sample rate just to speed up the example
 # and use less memory, resampling to 0.6 Hz is a better choice for full
 # analyses.
+#
+# .. note::
+#
+#    The nilearn library does not allow backslash characters in the condition
+#    name. So we must replace the backslash with an underscore to ensure the
+#    GLM computation is successful. Hopefully future versions of MNE-NIRS will
+#    automatically handle these characters, see https://github.com/mne-tools/mne-nirs/issues/420
+#    for more information. In the meantime use the following code to replace the
+#    illegal characters.
 
 
 def individual_analysis(bids_path, ID):
 
     raw_intensity = read_raw_bids(bids_path=bids_path, verbose=False)
+    # sanitize event names
+    raw_intensity.annotations.description[:] = [
+        d.replace('/', '_') for d in raw_intensity.annotations.description]
 
     # Convert signal to haemoglobin and resample
     raw_od = optical_density(raw_intensity)
@@ -218,7 +230,7 @@ def individual_analysis(bids_path, ID):
     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_LvR = basic_conts['Tapping_Left'] - basic_conts['Tapping_Right']
 
     # Compute defined contrast
     contrast = glm_est.compute_contrast(contrast_LvR)
@@ -257,9 +269,9 @@ def individual_analysis(bids_path, ID):
     raw_haemo, roi, channel, con = individual_analysis(bids_path, sub)
 
     # Append individual results to all participants
-    df_roi = df_roi.append(roi)
-    df_cha = df_cha.append(channel)
-    df_con = df_con.append(con)
+    df_roi = pd.concat([df_roi, roi], ignore_index=True)
+    df_cha = pd.concat([df_cha, channel], ignore_index=True)
+    df_con = pd.concat([df_con, con], ignore_index=True)
 
 
 # %%
@@ -273,7 +285,7 @@ def individual_analysis(bids_path, ID):
 # We can already see that the control condition is always near zero,
 # and that the responses look to be contralateral to the tapping hand.
 
-grp_results = df_roi.query("Condition in ['Control', 'Tapping/Left', 'Tapping/Right']")
+grp_results = df_roi.query("Condition in ['Control', 'Tapping_Left', 'Tapping_Right']")
 grp_results = grp_results.query("Chroma in ['hbo']")
 
 ggplot(grp_results, aes(x='Condition', y='theta', color='ROI', shape='ROI')) \
@@ -316,7 +328,7 @@ def individual_analysis(bids_path, ID):
 # model by using the code
 # `roi_model = rlm('theta ~ -1 + ROI:Condition:Chroma', grp_results).fit()`.
 
-grp_results = df_roi.query("Condition in ['Control','Tapping/Left', 'Tapping/Right']")
+grp_results = df_roi.query("Condition in ['Control','Tapping_Left', 'Tapping_Right']")
 
 roi_model = smf.mixedlm("theta ~ -1 + ROI:Condition:Chroma",
                         grp_results, groups=grp_results["ID"]).fit(method='nm')
@@ -346,7 +358,7 @@ def individual_analysis(bids_path, ID):
 # contralaterally dominant responses) and there is no significant
 # effect of gender.
 
-grp_results = df_roi.query("Condition in ['Tapping/Left', 'Tapping/Right']")
+grp_results = df_roi.query("Condition in ['Tapping_Left', 'Tapping_Right']")
 grp_results = grp_results.query("Chroma in ['hbo']")
 grp_results = grp_results.query("ROI in ['Right_Hemisphere']")
 
@@ -368,7 +380,7 @@ def individual_analysis(bids_path, ID):
 # Filled symbols represent HbO, unfilled symbols represent HbR.
 
 # Regenerate the results from the original group model above
-grp_results = df_roi.query("Condition in ['Control','Tapping/Left', 'Tapping/Right']")
+grp_results = df_roi.query("Condition in ['Control','Tapping_Left', 'Tapping_Right']")
 roi_model = smf.mixedlm("theta ~ -1 + ROI:Condition:Chroma",
                         grp_results, groups=grp_results["ID"]).fit(method='nm')
 
@@ -401,7 +413,7 @@ def individual_analysis(bids_path, ID):
                          gridspec_kw=dict(width_ratios=[1, 1]))
 
 # Cut down the dataframe just to the conditions we are interested in
-ch_summary = df_cha.query("Condition in ['Tapping/Left', 'Tapping/Right']")
+ch_summary = df_cha.query("Condition in ['Tapping_Left', 'Tapping_Right']")
 ch_summary = ch_summary.query("Chroma in ['hbo']")
 
 # Run group level model and convert to dataframe
@@ -411,17 +423,17 @@ def individual_analysis(bids_path, ID):
 
 # Plot the two conditions
 plot_glm_group_topo(raw_haemo.copy().pick(picks="hbo"),
-                    ch_model_df.query("Condition in ['Tapping/Left']"),
+                    ch_model_df.query("Condition in ['Tapping_Left']"),
                     colorbar=False, axes=axes[0, 0],
                     vmin=0, vmax=20, cmap=mpl.cm.Oranges)
 
 plot_glm_group_topo(raw_haemo.copy().pick(picks="hbo"),
-                    ch_model_df.query("Condition in ['Tapping/Right']"),
+                    ch_model_df.query("Condition in ['Tapping_Right']"),
                     colorbar=True, axes=axes[0, 1],
                     vmin=0, vmax=20, cmap=mpl.cm.Oranges)
 
 # Cut down the dataframe just to the conditions we are interested in
-ch_summary = df_cha.query("Condition in ['Tapping/Left', 'Tapping/Right']")
+ch_summary = df_cha.query("Condition in ['Tapping_Left', 'Tapping_Right']")
 ch_summary = ch_summary.query("Chroma in ['hbr']")
 
 # Run group level model and convert to dataframe
@@ -431,11 +443,11 @@ def individual_analysis(bids_path, ID):
 
 # Plot the two conditions
 plot_glm_group_topo(raw_haemo.copy().pick(picks="hbr"),
-                    ch_model_df.query("Condition in ['Tapping/Left']"),
+                    ch_model_df.query("Condition in ['Tapping_Left']"),
                     colorbar=False, axes=axes[1, 0],
                     vmin=-10, vmax=0, cmap=mpl.cm.Blues_r)
 plot_glm_group_topo(raw_haemo.copy().pick(picks="hbr"),
-                    ch_model_df.query("Condition in ['Tapping/Right']"),
+                    ch_model_df.query("Condition in ['Tapping_Right']"),
                     colorbar=True, axes=axes[1, 1],
                     vmin=-10, vmax=0, cmap=mpl.cm.Blues_r)
 
@@ -498,7 +510,7 @@ def individual_analysis(bids_path, ID):
 
 # Run model code as above
 clim = dict(kind='value', pos_lims=(0, 11.5, 17))
-for idx, cond in enumerate(['Tapping/Left', 'Tapping/Right']):
+for idx, cond in enumerate(['Tapping_Left', 'Tapping_Right']):
 
     # Run same model as explained in the sections above
     ch_summary = df_cha.query("Condition in [@cond]")
@@ -521,7 +533,7 @@ def individual_analysis(bids_path, ID):
 # Sometimes a reviewer wants a long table of results per channel.
 # This can be generated from the statistics dataframe.
 
-ch_summary = df_cha.query("Condition in ['Tapping/Left', 'Tapping/Right']")
+ch_summary = df_cha.query("Condition in ['Tapping_Left', 'Tapping_Right']")
 ch_summary = ch_summary.query("Chroma in ['hbo']")
 
 # Run group level model and convert to dataframe
diff --git a/examples/general/plot_13_fir_glm.py b/examples/general/plot_13_fir_glm.py
index f8ad382fe..42830e740 100644
--- a/examples/general/plot_13_fir_glm.py
+++ b/examples/general/plot_13_fir_glm.py
@@ -105,6 +105,9 @@
 def analysis(fname, ID):
 
     raw_intensity = read_raw_bids(bids_path=fname, verbose=False)
+    # sanitize event names
+    raw_intensity.annotations.description[:] = [
+        d.replace('/', '_') for d in raw_intensity.annotations.description]
 
     # Convert signal to haemoglobin and just keep hbo
     raw_od = optical_density(raw_intensity)
@@ -163,7 +166,7 @@ def analysis(fname, ID):
 
     df_individual, raw, dm = analysis(bids_path, ID)
 
-    df = df.append(df_individual)
+    df = pd.concat([df, df_individual])
 
 
 # %%
@@ -177,14 +180,14 @@ def analysis(fname, ID):
 # Keep only tapping and FIR delay information in the dataframe
 # I.e., for this example we are not interest in the drift coefficients,
 # short channel information, or control conditions.
-df["isTapping"] = ["Tapping/Right" in n for n in df["Condition"]]
+df["isTapping"] = ["Tapping_Right" in n for n in df["Condition"]]
 df["isDelay"] = ["delay" in n for n in df["Condition"]]
 df = df.query("isDelay in [True]")
 df = df.query("isTapping in [True]")
 # Make a new column that stores the condition name for tidier model below
 df.loc[df["isTapping"] == True, "TidyCond"] = "Tapping"
 # Finally, extract the FIR delay in to its own column in data frame
-df.loc[:, "delay"] = [n.split('_')[2] for n in df.Condition]
+df.loc[:, "delay"] = [n.split('_')[-1] for n in df.Condition]
 
 # To simplify this example we will only look at the right hand tapping
 # condition so we now remove the left tapping conditions from the
@@ -267,10 +270,11 @@ def analysis(fname, ID):
 # Extract the time scale for plotting.
 # Set time zero to be the onset of the finger tapping.
 index_values = dm_cond_scaled_hbo.index - np.ceil(raw.annotations.onset[0])
+index_values = np.asarray(index_values)
 
 # Plot the result
-axes[0].plot(index_values, dm_cond)
-axes[1].plot(index_values, dm_cond_scaled_hbo)
+axes[0].plot(index_values, np.asarray(dm_cond))
+axes[1].plot(index_values,np.asarray(dm_cond_scaled_hbo))
 axes[2].plot(index_values, np.sum(dm_cond_scaled_hbo, axis=1), 'r')
 axes[2].plot(index_values, np.sum(dm_cond_scaled_hbr, axis=1), 'b')
 
diff --git a/examples/general/plot_16_waveform_group.py b/examples/general/plot_16_waveform_group.py
index 19b72e8c4..885ea6476 100644
--- a/examples/general/plot_16_waveform_group.py
+++ b/examples/general/plot_16_waveform_group.py
@@ -317,17 +317,23 @@ def individual_analysis(bids_path):
 df = pd.DataFrame(columns=['ID', 'ROI', 'Chroma', 'Condition', 'Value'])
 
 for idx, evoked in enumerate(all_evokeds):
+    subj_id = 0
     for subj_data in all_evokeds[evoked]:
+        subj_id += 1
         for roi in rois:
             for chroma in ["hbo", "hbr"]:
-                subj_id = subj_data.info["subject_info"]['first_name']
+                assert subj_data.info["subject_info"]['first_name'] == \
+                    'mne_anonymize'  # have been anonymized
                 data = deepcopy(subj_data).pick(picks=rois[roi]).pick(chroma)
                 value = data.crop(tmin=5.0, tmax=7.0).data.mean() * 1.0e6
 
                 # Append metadata and extracted feature to dataframe
-                df = df.append({'ID': subj_id, 'ROI': roi, 'Chroma': chroma,
-                                'Condition': evoked, 'Value': value},
-                               ignore_index=True)
+                this_df = pd.DataFrame(
+                    {'ID': subj_id, 'ROI': roi, 'Chroma': chroma,
+                     'Condition': evoked, 'Value': value}, index=[0])
+                df = pd.concat([df, this_df], ignore_index=True)
+df.reset_index(inplace=True, drop=True)
+df['Value'] = pd.to_numeric(df['Value'])  # some Pandas have this as object
 
 # You can export the dataframe for analysis in your favorite stats program
 df.to_csv("stats-export.csv")
@@ -395,21 +401,22 @@ def individual_analysis(bids_path):
 # conditions and the model is applied.
 
 # Encode the ROIs as ipsi- or contralateral to the hand that is tapping.
-df["Hemishphere"] = "Unknown"
+df["Hemisphere"] = "Unknown"
 df.loc[(df["Condition"] == "Tapping/Right") &
-       (df["ROI"] == "Right_Hemisphere"), "Hemishphere"] = "Ipsilateral"
+       (df["ROI"] == "Right_Hemisphere"), "Hemisphere"] = "Ipsilateral"
 df.loc[(df["Condition"] == "Tapping/Right") &
-       (df["ROI"] == "Left_Hemisphere"), "Hemishphere"] = "Contralateral"
+       (df["ROI"] == "Left_Hemisphere"), "Hemisphere"] = "Contralateral"
 df.loc[(df["Condition"] == "Tapping/Left") &
-       (df["ROI"] == "Left_Hemisphere"), "Hemishphere"] = "Ipsilateral"
+       (df["ROI"] == "Left_Hemisphere"), "Hemisphere"] = "Ipsilateral"
 df.loc[(df["Condition"] == "Tapping/Left") &
-       (df["ROI"] == "Right_Hemisphere"), "Hemishphere"] = "Contralateral"
+       (df["ROI"] == "Right_Hemisphere"), "Hemisphere"] = "Contralateral"
 
 # Subset the data for example model
 input_data = df.query("Condition in ['Tapping/Right', 'Tapping/Left']")
 input_data = input_data.query("Chroma in ['hbo']")
+assert len(input_data)
 
-roi_model = smf.mixedlm("Value ~ Hemishphere", input_data,
+roi_model = smf.mixedlm("Value ~ Hemisphere", input_data,
                         groups=input_data["ID"]).fit()
 roi_model.summary()
 
diff --git a/examples/general/plot_19_snirf.py b/examples/general/plot_19_snirf.py
index 03b2391ca..506c280de 100644
--- a/examples/general/plot_19_snirf.py
+++ b/examples/general/plot_19_snirf.py
@@ -91,6 +91,6 @@
 # https://github.com/BUNPC/pysnirf2. Below we demonstrate that the files created
 # by MNE-NIRS are compliant with the specification.
 
-valid, result = pysnirf2.validateSnirf('test_raw.snirf')
-assert valid
-print(result)
+result = pysnirf2.validateSnirf('test_raw.snirf')
+assert result.is_valid()
+result.display()
diff --git a/examples/general/plot_30_frequency.py b/examples/general/plot_30_frequency.py
index 43850f4b0..9d861d95a 100644
--- a/examples/general/plot_30_frequency.py
+++ b/examples/general/plot_30_frequency.py
@@ -51,8 +51,8 @@
                                      verbose=True).load_data()
 new_des = [des for des in raw_intensity.annotations.description]
 new_des = ['Control' if x == "1.0" else x for x in new_des]
-new_des = ['Tapping/Left' if x == "2.0" else x for x in new_des]
-new_des = ['Tapping/Right' if x == "3.0" else x for x in new_des]
+new_des = ['Tapping_Left' if x == "2.0" else x for x in new_des]
+new_des = ['Tapping_Right' if x == "3.0" else x for x in new_des]
 annot = mne.Annotations(raw_intensity.annotations.onset,
                         raw_intensity.annotations.duration * 5., new_des)
 raw_intensity.set_annotations(annot)
@@ -90,8 +90,8 @@
 # Overwrite the first NIRS channel with the expected response.
 # Rescale to be in expected units of uM.
 hrf = raw_haemo.copy().pick(picks=[0])
-hrf._data[0] = 1e-6 * (design_matrix['Tapping/Left'] +
-                       design_matrix['Tapping/Right']).T
+hrf._data[0] = 1e-6 * (design_matrix['Tapping_Left'] +
+                       design_matrix['Tapping_Right']).T
 hrf.pick(picks='hbo').plot_psd(average=True, fmax=2, xscale='log',
                                color='r', show=False)
 
@@ -129,7 +129,7 @@
 # is still visible.
 
 events, _ = mne.events_from_annotations(raw_haemo)
-event_dict = {'Tapping/Left': 1, 'Tapping/Right': 2}
+event_dict = {'Tapping_Left': 1, 'Tapping_Right': 2}
 reject_criteria = dict(hbo=120e-6)
 tmin, tmax = -5, 15
 epochs = mne.Epochs(raw_haemo, events, event_id=event_dict,
diff --git a/examples/general/plot_70_visualise_brain.py b/examples/general/plot_70_visualise_brain.py
index 9faeae12c..70f8a5730 100644
--- a/examples/general/plot_70_visualise_brain.py
+++ b/examples/general/plot_70_visualise_brain.py
@@ -237,6 +237,9 @@
 def individual_analysis(bids_path, ID):
 
     raw_intensity = read_raw_bids(bids_path=bids_path, verbose=False)
+     # sanitize event names
+    raw_intensity.annotations.description[:] = [
+        d.replace('/', '_') for d in raw_intensity.annotations.description]
 
     # Convert signal to haemoglobin and resample
     raw_od = optical_density(raw_intensity)
@@ -285,9 +288,10 @@ def individual_analysis(bids_path, ID):
     raw_haemo, channel = individual_analysis(bids_path, sub)
 
     # Append individual results to all participants
-    df_cha = df_cha.append(channel)
+    df_cha = pd.concat([df_cha, channel], ignore_index=True)
 
-ch_summary = df_cha.query("Condition in 'Tapping/Right'")
+ch_summary = df_cha.query("Condition in ['Tapping_Right']")
+assert len(ch_summary)
 ch_summary = ch_summary.query("Chroma in ['hbo']")
 ch_model = smf.mixedlm("theta ~ -1 + ch_name", ch_summary,
                        groups=ch_summary["ID"]).fit(method='nm')
diff --git a/examples/general/plot_80_save_read_glm.py b/examples/general/plot_80_save_read_glm.py
index a9cde83b3..0c2533935 100644
--- a/examples/general/plot_80_save_read_glm.py
+++ b/examples/general/plot_80_save_read_glm.py
@@ -172,7 +172,7 @@ def individual_analysis(bids_path):
     individual_results["ID"] = sub
 
     # Append individual results to larger dataframe
-    df = df.append(individual_results)
+    df = pd.concat([df, individual_results], ignore_index=True)
 
 
 # %%
diff --git a/mne_nirs/datasets/audio_or_visual_speech/_audio_or_visual_speech.py b/mne_nirs/datasets/audio_or_visual_speech/_audio_or_visual_speech.py
index abad56c6c..42ace770c 100644
--- a/mne_nirs/datasets/audio_or_visual_speech/_audio_or_visual_speech.py
+++ b/mne_nirs/datasets/audio_or_visual_speech/_audio_or_visual_speech.py
@@ -4,7 +4,6 @@
 # This downloads my fnirs audio or video speech dataset.
 
 import os
-import tempfile
 import shutil
 import pooch
 from functools import partial
@@ -74,9 +73,8 @@ def data_path(path=None, force_update=False, update_path=True, download=True,
     bad_name = os.path.join(dpath, '2021-fNIRS-Audio-visual-speech-'
                                    'Broad-vs-restricted-regions')
     if os.path.isdir(bad_name):
-        tmppath = tempfile.mkdtemp()
-        os.rename(bad_name, tmppath)
+        os.rename(bad_name, dpath + '.true')
         shutil.rmtree(dpath)
-        os.rename(tmppath, dpath)
+        os.rename(dpath + '.true', dpath)
 
     return dpath
diff --git a/mne_nirs/datasets/block_speech_noise/_block_speech_noise.py b/mne_nirs/datasets/block_speech_noise/_block_speech_noise.py
index 175736bde..5434bb739 100644
--- a/mne_nirs/datasets/block_speech_noise/_block_speech_noise.py
+++ b/mne_nirs/datasets/block_speech_noise/_block_speech_noise.py
@@ -4,7 +4,6 @@
 # This downloads my fnirs speech and noise dataset.
 
 import os
-import tempfile
 import shutil
 import pooch
 from functools import partial
@@ -72,9 +71,8 @@ def data_path(path=None, force_update=False, update_path=True, download=True,
     bad_name = os.path.join(dpath, '2021-fNIRS-Analysis-Methods-'
                                    'Passive-Auditory')
     if os.path.isdir(bad_name):
-        tmppath = tempfile.mkdtemp()
-        os.rename(bad_name, tmppath)
+        os.rename(bad_name, dpath + '.true')
         shutil.rmtree(dpath)
-        os.rename(tmppath, dpath)
+        os.rename(dpath + '.true', dpath)
 
     return dpath
diff --git a/mne_nirs/datasets/fnirs_motor_group/fnirs_motor_group.py b/mne_nirs/datasets/fnirs_motor_group/fnirs_motor_group.py
index 865413036..43b7a2e94 100644
--- a/mne_nirs/datasets/fnirs_motor_group/fnirs_motor_group.py
+++ b/mne_nirs/datasets/fnirs_motor_group/fnirs_motor_group.py
@@ -5,7 +5,6 @@
 # and the dataset complies I will add this to the MNE-Python library.
 
 import os
-import tempfile
 import shutil
 import pooch
 from functools import partial
@@ -57,7 +56,7 @@ def data_path(path=None, force_update=False, update_path=True, download=True,
 
     dataset_params = dict(
         archive_name='BIDS-NIRS-Tapping-master.zip',
-        hash='md5:3ed36f9cd56ae7f0c54826fcf5b5eefe',
+        hash='md5:d40529b85d4b7ebdf7243168c0e8c390',
         url='https://github.com/rob-luke/BIDS-NIRS-Tapping/archive/master.zip',
         folder_name='fNIRS-motor-group',
         dataset_name='fnirs_motor_group',
@@ -72,9 +71,8 @@ def data_path(path=None, force_update=False, update_path=True, download=True,
     # Do some wrangling to deal with nested directories
     bad_name = os.path.join(dpath, 'BIDS-NIRS-Tapping-master')
     if os.path.isdir(bad_name):
-        tmppath = tempfile.mkdtemp()
-        os.rename(bad_name, tmppath)
+        os.rename(bad_name, dpath + '.true')
         shutil.rmtree(dpath)
-        os.rename(tmppath, dpath)
+        os.rename(dpath + '.true', dpath)
 
     return dpath
diff --git a/mne_nirs/io/fold/_fold.py b/mne_nirs/io/fold/_fold.py
index 69d413ad4..c2677b61d 100644
--- a/mne_nirs/io/fold/_fold.py
+++ b/mne_nirs/io/fold/_fold.py
@@ -144,7 +144,8 @@ def fold_landmark_specificity(raw, landmark, fold_files=[None],
 
     fold_tbl = pd.DataFrame()
     for fname in fold_files:
-        fold_tbl = fold_tbl.append(_read_fold_xls(fname, atlas=atlas))
+        fold_tbl = pd.concat([fold_tbl, _read_fold_xls(fname, atlas=atlas)],
+                             ignore_index=True)
 
     specificity = np.zeros(len(raw.ch_names))
     for cidx in range(len(raw.ch_names)):
@@ -205,7 +206,8 @@ def fold_channel_specificity(raw, fold_files=[None], atlas="Juelich"):
 
     fold_tbl = pd.DataFrame()
     for fname in fold_files:
-        fold_tbl = fold_tbl.append(_read_fold_xls(fname, atlas=atlas))
+        fold_tbl = pd.concat([fold_tbl, _read_fold_xls(fname, atlas=atlas)],
+                             ignore_index=True)
 
     chan_spec = list()
     for cidx in range(len(raw.ch_names)):
diff --git a/mne_nirs/preprocessing/_mayer.py b/mne_nirs/preprocessing/_mayer.py
index 71d86a9dd..c8c39651a 100644
--- a/mne_nirs/preprocessing/_mayer.py
+++ b/mne_nirs/preprocessing/_mayer.py
@@ -103,7 +103,8 @@ def quantify_mayer_fooof(raw, num_oscillations=1, centre_frequency=0.01,
             data["Chromaphore"] = chroma
             data = {**data, **extra_df_fields}
 
-            df = df.append(pd.DataFrame(data, index=[0]), ignore_index=True)
+            df = pd.concat([df, pd.DataFrame(data, index=[0])],
+                           ignore_index=True)
 
     return df
 
diff --git a/mne_nirs/statistics/_glm_level_first.py b/mne_nirs/statistics/_glm_level_first.py
index 65b2a1408..d3bd6c936 100644
--- a/mne_nirs/statistics/_glm_level_first.py
+++ b/mne_nirs/statistics/_glm_level_first.py
@@ -18,7 +18,7 @@
 
 from mne.channels.channels import ContainsMixin
 from mne.utils import fill_doc, warn, verbose, check_fname, _validate_type
-from mne.io.pick import _picks_to_idx
+from mne.io.pick import _picks_to_idx, pick_info
 from mne.io.constants import FIFF
 from mne import Info
 
@@ -267,8 +267,7 @@ def pick(self, picks, exclude=()):
         """
         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)
+        pick_info(self.info, picks, copy=False)
         self._data = {key: self._data[key] for key in self.info.ch_names}
         return self
 
@@ -414,10 +413,13 @@ def to_dataframe_region_of_interest(self, group_by, condition,
         tidy = pd.DataFrame()
         for cond in condition:
             cond_idx = where([c == cond for c in self.design.columns])[0]
+            if not len(cond_idx):
+                raise KeyError(f'condition {repr(cond)} not found in '
+                               f'self.design.columns: {self.design.columns}')
 
             roi = _glm_region_of_interest(self._data, group_by,
                                           cond_idx, cond, weighted)
-            tidy = tidy.append(roi)
+            tidy = pd.concat([tidy, roi])
 
         if weighted is True:
             tidy["Weighted"] = "Inverse standard error"
@@ -520,8 +522,13 @@ def surface_projection(self, chroma="hbo", condition=None,
         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_use = df.query("Condition in @condition")
+        if len(df_use) == 0:
+            raise KeyError(
+                f'condition={repr(condition)} not found in conditions: '
+                f'{sorted(set(df["Condition"]))}')
+        df = df_use
+        df = df.query("Chroma in @chroma").copy()
         df["theta"] = df["theta"] * 1e6
 
         info = self.copy().pick(chroma).info
diff --git a/mne_nirs/statistics/_roi.py b/mne_nirs/statistics/_roi.py
index e30941ad2..6720bb004 100644
--- a/mne_nirs/statistics/_roi.py
+++ b/mne_nirs/statistics/_roi.py
@@ -135,15 +135,18 @@ def _glm_region_of_interest(stats, group_by, cond_idx,
 
             p = 2 * ss.t.cdf(-1.0 * np.abs(t), df=dfe)
 
-            df = df.append({'ROI': roi_name,
-                            'Condition': cond_name,
-                            'Chroma': chroma,
-                            'theta': theta / 1.0e6,
-                            'se': s,
-                            't': t,
-                            'dfe': dfe,
-                            'p': p, },
-                           ignore_index=True)
+            this_df = pd.DataFrame(
+                {'ROI': roi_name,
+                 'Condition': cond_name,
+                 'Chroma': chroma,
+                 'theta': theta / 1.0e6,
+                 'se': s,
+                 't': t,
+                 'dfe': dfe,
+                 'p': p, }, index=[0])
+            df = pd.concat([df, this_df], ignore_index=True)
+
+    df.reset_index(inplace=True, drop=True)
 
     return df
 
diff --git a/mne_nirs/statistics/_statsmodels.py b/mne_nirs/statistics/_statsmodels.py
index 1ecba9b9a..0ba05ed49 100644
--- a/mne_nirs/statistics/_statsmodels.py
+++ b/mne_nirs/statistics/_statsmodels.py
@@ -31,10 +31,10 @@ def expand_summary_dataframe(summary):
         col_names.append(col_name)
 
     # Fill in values
-    indices = summary.copy(deep=True).index
     if 'Group Var' in summary.index:
         summary = summary[:-1]
-        indices = summary.index
+    summary = summary.copy(deep=True)
+    indices = summary.index
     for row_idx, row in enumerate(indices):
         col_vals = row.split(':')
         for col_idx, col in enumerate(col_names):
diff --git a/mne_nirs/utils/_io.py b/mne_nirs/utils/_io.py
index 68913121e..fe2a010ba 100644
--- a/mne_nirs/utils/_io.py
+++ b/mne_nirs/utils/_io.py
@@ -66,26 +66,29 @@ def glm_to_tidy(info, statistic, design_matrix, wide=True, order=None):
 def _tidy_Contrast(data, glm_est, design_matrix):
     df = pd.DataFrame()
     for idx, ch in enumerate(data.ch_names):
-        df = df.append({'ch_name': ch,
-                        'ContrastType': glm_est.contrast_type,
-                        'variable': "effect",
-                        'value': glm_est.effect[0][idx]},
-                       ignore_index=True)
-        df = df.append({'ch_name': ch,
-                        'ContrastType': glm_est.contrast_type,
-                        'variable': "p_value",
-                        'value': glm_est.p_value()[idx]},
-                       ignore_index=True)
-        df = df.append({'ch_name': ch,
-                        'ContrastType': glm_est.contrast_type,
-                        'variable': "stat",
-                        'value': glm_est.stat()[idx]},
-                       ignore_index=True)
-        df = df.append({'ch_name': ch,
-                        'ContrastType': glm_est.contrast_type,
-                        'variable': "z_score",
-                        'value': glm_est.z_score()[idx]},
-                       ignore_index=True)
+        df = pd.concat([
+            df,
+            pd.DataFrame({'ch_name': ch,
+                          'ContrastType': glm_est.contrast_type,
+                          'variable': "effect",
+                          'value': glm_est.effect[0][idx]},
+                         index=[0]),
+            pd.DataFrame({'ch_name': ch,
+                          'ContrastType': glm_est.contrast_type,
+                          'variable': "p_value",
+                          'value': glm_est.p_value()[idx]},
+                         index=[1]),
+            pd.DataFrame({'ch_name': ch,
+                          'ContrastType': glm_est.contrast_type,
+                          'variable': "stat",
+                          'value': glm_est.stat()[idx]},
+                         index=[2]),
+            pd.DataFrame({'ch_name': ch,
+                          'ContrastType': glm_est.contrast_type,
+                          'variable': "z_score",
+                          'value': glm_est.z_score()[idx]},
+                         index=[3]),
+        ], ignore_index=True)
     return df
 
 
diff --git a/mne_nirs/visualisation/_plot_GLM_topo.py b/mne_nirs/visualisation/_plot_GLM_topo.py
index 28e15aa02..40d4c2e9f 100644
--- a/mne_nirs/visualisation/_plot_GLM_topo.py
+++ b/mne_nirs/visualisation/_plot_GLM_topo.py
@@ -13,7 +13,7 @@
 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, _get_channel_types
+from mne.io.pick import _picks_to_idx, _get_channel_types, pick_info
 
 
 def plot_glm_topo(inst, glm_estimates, design_matrix,
@@ -383,8 +383,7 @@ def plot_glm_group_topo(inst, statsmodel_df,
 def _handle_overlaps(info, t, sphere, estimates):
     """Prepare for topomap including merging channels"""
     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)
+    info_subset = pick_info(info, picks)
     _, pos, merge_channels, ch_names, ch_type, sphere, clip_origin = \
         mne.viz.topomap._prepare_topomap_plot(info_subset, t, sphere=sphere)
     estmrg, ch_names = _merge_ch_data(estimates.copy()[picks], t, ch_names)
diff --git a/requirements.txt b/requirements.txt
index 82f65925b..fc4f55002 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,16 +2,24 @@
 numpy>=1.15.4
 scipy>=1.1.0
 matplotlib
+tqdm
+pooch>=1.5
+decorator
+h5io
+packaging
+pymatreader
 pyqt5>=5.10,<5.14; platform_system == "Darwin"
 pyqt5>=5.10,!=5.15.2,!=5.15.3; platform_system == "Linux"
-pyqt5>=5.10,!=5.15.3; platform_system != "Linux" and platform_system != "Darwin"
+pyqt5>=5.10,!=5.15.3; platform_system == "Windows"
 pyqt5-sip
+pyobjc-framework-Cocoa>=5.2.0; platform_system=="Darwin"
 sip
 scikit-learn
 nibabel
 numba
 h5py
-pandas==1.3.2
+jinja2
+pandas
 numexpr
 jupyter
 python-picard
@@ -19,20 +27,17 @@ statsmodels
 joblib
 psutil
 dipy
-vtk; platform_system == "Darwin"
-vtk; platform_system == "Linux"
-vtk<=9.0.1; platform_system == "Windows"
+vtk
 nilearn
 xlrd
 imageio>=2.6.1
 imageio-ffmpeg>=0.4.1
-pyvista>=0.30
+traitlets
+pyvista>=0.32
 pyvistaqt>=0.4
-tqdm
 mffpy>=0.5.7
 ipywidgets
 ipyvtklink
-pooch
 
 # mne-nirs specific requirements
 https://github.com/mne-tools/mne-python/archive/main.zip
diff --git a/requirements_doc.txt b/requirements_doc.txt
index 4a1d8ae1a..96a9195fc 100644
--- a/requirements_doc.txt
+++ b/requirements_doc.txt
@@ -5,11 +5,8 @@ pydata-sphinx-theme>=0.6.3
 https://github.com/sphinx-gallery/sphinx-gallery/archive/master.zip
 sphinxcontrib-bibtex>=2.3
 memory_profiler
-neo
 seaborn
 sphinx_copybutton
-# https://github.com/mne-tools/mne-bids/archive/main.zip
-pyxdf
 
 # mne-nirs specific requirements
 sphinx-multiversion
@@ -26,5 +23,5 @@ seaborn
 fooof
 pysnirf2>=0.4.2
 statsmodels
-lets-plot
+lets-plot; python_version < '3.10'
 dabest
diff --git a/setup.cfg b/setup.cfg
index d688c0384..9948d1b02 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -9,4 +9,9 @@ exclude = __init__.py,*externals*,constants.py,fixes.py
 ignore = W504,I100,I101,I201,N806
 
 [tool:pytest]
-addopts = --doctest-modules
+addopts =
+    --durations=20 --doctest-modules -ra --cov-report= --tb=short
+    --doctest-ignore-import-errors --junit-xml=junit-results.xml
+    --ignore=doc --ignore=examples --ignore=tools
+    --color=yes --capture=sys
+junit_family = xunit2
diff --git a/tools/github_actions_dependencies.sh b/tools/github_actions_dependencies.sh
index 6591a2513..03029455d 100755
--- a/tools/github_actions_dependencies.sh
+++ b/tools/github_actions_dependencies.sh
@@ -1,37 +1,48 @@
 #!/bin/bash -ef
 
+STD_ARGS="--progress-bar off --upgrade"
+EXTRA_ARGS=""
 if [ ! -z "$CONDA_ENV" ]; then
 	pip uninstall -yq mne
 elif [ ! -z "$CONDA_DEPENDENCIES" ]; then
 	conda install -y $CONDA_DEPENDENCIES
 else
 	# Changes here should also go in the interactive_test CircleCI job
-	python -m pip install --progress-bar off --upgrade pip setuptools wheel
+	python -m pip install $STD_ARGS pip setuptools wheel
 	echo "Numpy"
 	pip uninstall -yq numpy
 	echo "Date utils"
 	# https://pip.pypa.io/en/latest/user_guide/#possible-ways-to-reduce-backtracking-occurring
-	pip install --progress-bar off --upgrade --pre --only-binary ":all:" python-dateutil pytz joblib threadpoolctl six
+	pip install $STD_ARGS --pre --only-binary ":all:" python-dateutil pytz joblib threadpoolctl six
 	echo "PyQt5"
-	pip install --progress-bar off --upgrade --pre --only-binary ":all:" --no-deps --extra-index-url https://www.riverbankcomputing.com/pypi/simple PyQt5 PyQt5-sip PyQt5-Qt5
-	echo "NumPy etc."
-	pip install --progress-bar off --upgrade --pre --only-binary ":all:" --no-deps -i "https://pypi.anaconda.org/scipy-wheels-nightly/simple" numpy "scikit-learn>=0.24.2" dipy
+	pip install $STD_ARGS --pre --only-binary ":all:" --no-deps --extra-index-url https://www.riverbankcomputing.com/pypi/simple PyQt5 PyQt5-sip PyQt5-Qt5
+	echo "NumPy/SciPy/pandas etc."
+	# TODO: Currently missing dipy for 3.10 https://github.com/dipy/dipy/issues/2489
+	pip install $STD_ARGS --pre --only-binary ":all:" --no-deps -i "https://pypi.anaconda.org/scipy-wheels-nightly/simple" numpy scipy pandas "scikit-learn>=0.24.2" statsmodels
 	echo "H5py, pillow, matplotlib"
-	pip install --progress-bar off --upgrade --pre --only-binary ":all:" --no-deps -f "https://7933911d6844c6c53a7d-47bd50c35cd79bd838daf386af554a83.ssl.cf2.rackcdn.com" h5py pillow matplotlib
+	pip install $STD_ARGS --pre --only-binary ":all:" --no-deps -f "https://7933911d6844c6c53a7d-47bd50c35cd79bd838daf386af554a83.ssl.cf2.rackcdn.com" h5py pillow matplotlib
 	echo "Numba, nilearn"
-	pip install --progress-bar off --upgrade --pre --only-binary ":all:" numba llvmlite nilearn
+	pip install $STD_ARGS --pre --only-binary ":all:" numba llvmlite https://github.com/nilearn/nilearn/zipball/main
 	echo "VTK"
-	# built using vtk master branch on an Ubuntu 18.04.5 VM and uploaded to OSF,
-	# can't use VTK's pre wheel because it breaks Mayavi (probably event processing)
-	wget -q https://osf.io/kej3v/download -O vtk-9.0.20201117-cp39-cp39-linux_x86_64.whl
-	pip install --progress-bar off vtk-9.0.20201117-cp39-cp39-linux_x86_64.whl
+	# Have to use our own version until VTK releases a 3.10 build
+	if `python -c "import sys; assert sys.version_info >= (3, 10), sys.version_info"`; then
+		wget -q https://osf.io/hjyvx/download -O vtk-9.1.20211213.dev0-cp310-cp310-linux_x86_64.whl
+		pip install $STD_ARGS --pre --only-binary ":all:" vtk-9.1.20211213.dev0-cp310-cp310-linux_x86_64.whl
+	else
+		pip install $STD_ARGS --pre --only-binary ":all:" vtk
+	fi;
 	echo "PyVista"
 	pip install --progress-bar off https://github.com/pyvista/pyvista/zipball/main
 	echo "pyvistaqt"
 	pip install --progress-bar off https://github.com/pyvista/pyvistaqt/zipball/main
 	echo "imageio-ffmpeg, xlrd, mffpy"
 	pip install --progress-bar off --pre imageio-ffmpeg xlrd mffpy
+	if [ "$OSTYPE" == "darwin"* ]; then
+	  echo "pyobjc-framework-Cocoa"
+	  pip install --progress-bar off pyobjc-framework-Cocoa>=5.2.0
+	fi
+	EXTRA_ARGS="--pre"
 fi
 
-pip install --progress-bar off -r requirements.txt
-pip install --progress-bar off --upgrade -r requirements_testing.txt
+pip install $STD_ARGS $EXTRA_ARGS -r requirements.txt
+pip install $STD_ARGS $EXTRA_ARGS -r requirements_testing.txt


Test patch:
diff --git a/mne_nirs/conftest.py b/mne_nirs/conftest.py
new file mode 100644
index 000000000..825c1d538
--- /dev/null
+++ b/mne_nirs/conftest.py
@@ -0,0 +1,79 @@
+# -*- coding: utf-8 -*-
+# Author: Eric Larson <larson.eric.d@gmail.com>
+#
+# License: BSD-3-Clause
+
+import warnings
+
+import pytest
+
+
+# most of this adapted from MNE-Python
+# TODO: Add to mne's allowed config keys:
+# MNE_DATASETS_FNIRSMOTORGROUP_PATH
+# MNE_DATASETS_FNIRSSPEECHNOISE_PATH
+# MNE_DATASETS_FNIRSAUDIOVISUALSPEECH_PATH
+
+def pytest_configure(config):
+    """Configure pytest options."""
+    # Markers
+    for marker in ('examples',):
+        config.addinivalue_line('markers', marker)
+    for fixture in ('matplotlib_config', 'close_all'):
+        config.addinivalue_line('usefixtures', fixture)
+
+    warning_lines = r"""
+    error::
+    ignore:.*np\.MachAr.*:DeprecationWarning
+    ignore:.*sysconfig module is deprecated.*:DeprecationWarning
+    ignore:.*nilearn.glm module is experimental.*:
+    ignore:.*Using or importing the ABCs from.*:
+    ignore:.*distutils Version classes are deprecated.*:
+    ignore:.*LUTSIZE was deprecated in Matplotlib.*:
+    ignore:.*pandas\.Int64Index is deprecated and will be removed.*:
+    ignore:.*Setting non-standard config type.*:
+    ignore:.*The MLE may be on the boundary.*:
+    ignore:.*The Hessian matrix at the estimated parameter values.*:
+    always:`np\..*is a deprecated alias for the builtin.*:DeprecationWarning
+    always::ResourceWarning
+    """  # noqa: E501
+    for warning_line in warning_lines.split('\n'):
+        warning_line = warning_line.strip()
+        if warning_line and not warning_line.startswith('#'):
+            config.addinivalue_line('filterwarnings', warning_line)
+
+
+@pytest.fixture(scope='session')
+def matplotlib_config():
+    """Configure matplotlib for viz tests."""
+    from matplotlib import cbook, use
+    want = 'agg'  # don't pop up windows
+    with warnings.catch_warnings(record=True):  # ignore warning
+        warnings.filterwarnings('ignore')
+        use(want, force=True)
+    import matplotlib.pyplot as plt
+    assert plt.get_backend() == want
+    # overwrite some params that can horribly slow down tests that
+    # users might have changed locally (but should not otherwise affect
+    # functionality)
+    plt.ioff()
+    plt.rcParams['figure.dpi'] = 100
+    plt.rcParams['figure.max_open_warning'] = 100
+
+    # Make sure that we always reraise exceptions in handlers
+    orig = cbook.CallbackRegistry
+
+    class CallbackRegistryReraise(orig):
+        def __init__(self, exception_handler=None):
+            super(CallbackRegistryReraise, self).__init__(exception_handler)
+
+    cbook.CallbackRegistry = CallbackRegistryReraise
+
+
+@pytest.fixture(autouse=True)
+def close_all():
+    """Close all matplotlib plots, regardless of test status."""
+    # This adds < 1 µS in local testing, and we have ~2500 tests, so ~2 ms max
+    import matplotlib.pyplot as plt
+    yield
+    plt.close('all')
diff --git a/mne_nirs/experimental_design/tests/test_experimental_design.py b/mne_nirs/experimental_design/tests/test_experimental_design.py
index 23552d659..9a0bde5a7 100644
--- a/mne_nirs/experimental_design/tests/test_experimental_design.py
+++ b/mne_nirs/experimental_design/tests/test_experimental_design.py
@@ -3,7 +3,7 @@
 # License: BSD (3-clause)
 
 import os
-import pytest
+
 import mne
 import mne_nirs
 import numpy as np
@@ -60,7 +60,6 @@ def test_create_boxcar():
     assert np.max(np.mean(bc, axis=1)) * bc.shape[1] == 1
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_create_design():
     raw_intensity = _load_dataset()
     raw_intensity.crop(450, 600)  # Keep the test fast
diff --git a/mne_nirs/io/fold/tests/test_fold.py b/mne_nirs/io/fold/tests/test_fold.py
index 04ba8c056..7a9ae9d23 100644
--- a/mne_nirs/io/fold/tests/test_fold.py
+++ b/mne_nirs/io/fold/tests/test_fold.py
@@ -5,11 +5,14 @@
 
 import os.path as op
 import pathlib
+
 import numpy as np
 import pandas as pd
+import pytest
 
 from mne.datasets.testing import data_path, requires_testing_data
 from mne.io import read_raw_nirx
+from mne.utils import check_version
 
 from mne_nirs.io.fold._fold import _generate_montage_locations,\
     _find_closest_standard_location, _read_fold_xls
@@ -23,7 +26,11 @@
 fname_nirx_15_3_short = op.join(data_path(download=False),
                                 'NIRx', 'nirscout', 'nirx_15_3_recording')
 
+requires_xlrd = pytest.mark.skipif(
+    not check_version('xlrd', '1.0'), reason='Requires xlrd >= 1.0')
+
 
+@requires_xlrd
 def test_channel_specificity():
     raw = read_raw_nirx(fname_nirx_15_3_short, preload=True)
     raw.pick(range(2))
@@ -32,6 +39,7 @@ def test_channel_specificity():
     assert res[0].shape == (6, 10)
 
 
+@requires_xlrd
 def test_landmark_specificity():
     raw = read_raw_nirx(fname_nirx_15_3_short, preload=True)
     res = fold_landmark_specificity(raw, "L Superior Frontal Gyrus",
@@ -41,6 +49,7 @@ def test_landmark_specificity():
     assert np.min(res) >= 0
 
 
+@requires_xlrd
 def test_fold_workflow():
     # Read raw data
     raw = read_raw_nirx(fname_nirx_15_3_short, preload=True)
@@ -68,6 +77,7 @@ def test_fold_workflow():
     assert specificity.values == 12.34
 
 
+@requires_xlrd
 def test_fold_reader():
     tbl = _read_fold_xls(foldfile, atlas="Juelich")
     assert isinstance(tbl, pd.DataFrame)
diff --git a/mne_nirs/io/snirf/tests/test_snirf.py b/mne_nirs/io/snirf/tests/test_snirf.py
index 94b8046d5..8e0358e8f 100644
--- a/mne_nirs/io/snirf/tests/test_snirf.py
+++ b/mne_nirs/io/snirf/tests/test_snirf.py
@@ -38,10 +38,10 @@ def test_snirf_write(fname, tmpdir):
     write_raw_snirf(raw_orig, test_file)
     raw = read_raw_snirf(test_file)
 
-    valid, result = validateSnirf(str(test_file))
-    if not valid:
+    result = validateSnirf(str(test_file))
+    if result.is_valid():
         result.display()
-    assert valid
+    assert result.is_valid()
 
     # Check annotations are the same
     assert_allclose(raw.annotations.onset, raw_orig.annotations.onset)
diff --git a/mne_nirs/simulation/tests/test_simulation.py b/mne_nirs/simulation/tests/test_simulation.py
index a7ea2f6f7..38dd84508 100644
--- a/mne_nirs/simulation/tests/test_simulation.py
+++ b/mne_nirs/simulation/tests/test_simulation.py
@@ -8,7 +8,6 @@
 from mne_nirs.experimental_design import make_first_level_design_matrix
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_simulate_NIRS_single_channel():
 
     raw = simulate_nirs_raw(sfreq=3., amplitude=1., sig_dur=300., stim_dur=5.,
@@ -27,7 +26,6 @@ def test_simulate_NIRS_single_channel():
                               stim_dur=5., isi_min=15., isi_max=45.)
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_simulate_NIRS_multi_channel():
 
     raw = simulate_nirs_raw(sfreq=3.,
diff --git a/mne_nirs/statistics/tests/test_glm_type.py b/mne_nirs/statistics/tests/test_glm_type.py
index 2c88b8a4a..b31be27c9 100644
--- a/mne_nirs/statistics/tests/test_glm_type.py
+++ b/mne_nirs/statistics/tests/test_glm_type.py
@@ -25,6 +25,8 @@ def _get_minimal_haemo_data(tmin=0, tmax=60):
     raw = mne.preprocessing.nirs.optical_density(raw)
     raw = mne.preprocessing.nirs.beer_lambert_law(raw, ppf=0.1)
     raw.resample(0.3)
+    raw.annotations.description[:] = [
+        'e' + d.replace('.', 'p') for d in raw.annotations.description]
     return raw
 
 
@@ -47,7 +49,8 @@ def _get_glm_contrast_result(tmin=60, tmax=400):
     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']
+    assert 'e1p' in basic_conts, sorted(basic_conts)
+    contrast_LvR = basic_conts['e1p'] - basic_conts['e2p']
 
     return glm_est.compute_contrast(contrast_LvR)
 
@@ -104,8 +107,8 @@ def test_results_glm_properties():
     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
+    with pytest.raises(RuntimeWarning, match='could not be picked'):
+        assert len(res.copy().pick(picks=["S1_D1 hbr", "S1_D1 XXX"])) == 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
@@ -132,8 +135,10 @@ def test_glm_scatter():
 
 def test_glm_surface_projection():
 
-    _get_glm_result(tmax=2974, tmin=0).surface_projection(condition="3.0",
-                                                          view="dorsal")
+    res = _get_glm_result(tmax=2974, tmin=0)
+    res.surface_projection(condition="e3p0", view="dorsal")
+    with pytest.raises(KeyError, match='not found in conditions'):
+        res.surface_projection(condition='foo')
 
 
 def test_results_glm_export_dataframe():
@@ -155,58 +160,60 @@ def test_results_glm_export_dataframe_region_of_interest():
     rois["A"] = [0, 2, 4]
 
     # Single ROI, single condition
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    with pytest.raises(KeyError, match=r'not found in self\.design\.col'):
+        res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
     assert df.shape == (1, 9)
-    assert df.Condition[0] == "1.0"
+    assert df.Condition[0] == "e1p0"
     assert df.ROI[0] == "A"
     assert df.Chroma[0] == "hbo"
 
     # Single ROI, multiple conditions
-    df = res.to_dataframe_region_of_interest(rois, ["1.0", "3.0", "drift_1"])
+    df = res.to_dataframe_region_of_interest(rois, ["e1p0", "e3p0", "drift_1"])
     assert df.shape == (3, 9)
-    assert (df.Condition == ["1.0", "3.0", "drift_1"]).all()
+    assert (df.Condition == ["e1p0", "e3p0", "drift_1"]).all()
     assert (df.ROI == ["A", "A", "A"]).all()
 
     # Single ROI, all conditions (default)
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
 
     # HbR only ROI
     rois["B"] = [1, 3, 5]
 
     # Multiple ROI, single condition
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
     assert df.shape == (2, 9)
-    assert df.Condition[0] == "1.0"
-    assert df.Condition[1] == "1.0"
+    assert df.Condition[0] == "e1p0"
+    assert df.Condition[1] == "e1p0"
     assert df.ROI[0] == "A"
     assert df.ROI[1] == "B"
     assert df.Chroma[0] == "hbo"
     assert df.Chroma[1] == "hbr"
 
     # Multiple ROI, multiple conditions
-    df = res.to_dataframe_region_of_interest(rois, ["1.0", "3.0", "drift_1"])
+    df = res.to_dataframe_region_of_interest(rois, ["e1p0", "e3p0", "drift_1"])
     assert df.shape == (6, 9)
 
     # Multiple ROI, all conditions (default)
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
     assert df.shape == (2, 9)
 
     # HbO and HbR ROI
     rois["C"] = [6, 7, 8]
 
     # Multiple ROI, single condition
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
     assert df.shape == (4, 9)
 
     # Multiple ROI, multiple conditions
-    df = res.to_dataframe_region_of_interest(rois, ["1.0", "3.0", "drift_1"])
+    df = res.to_dataframe_region_of_interest(rois, ["e1p0", "e3p0", "drift_1"])
     assert df.shape == (12, 9)
 
     # Multiple ROI, all conditions (default)
-    df = res.to_dataframe_region_of_interest(rois, "1.0")
+    df = res.to_dataframe_region_of_interest(rois, "e1p0")
 
     # With demographic information
-    df = res.to_dataframe_region_of_interest(rois, ["1.0", "3.0", "drift_1"],
+    df = res.to_dataframe_region_of_interest(rois, ["e1p0", "e3p0", "drift_1"],
                                              demographic_info=True)
     assert df.shape == (12, 10)
     assert "Sex" in df.columns
@@ -215,7 +222,8 @@ def test_results_glm_export_dataframe_region_of_interest():
 def test_results_glm_export_dataframe_region_of_interest_weighted():
 
     res = _get_glm_result(tmax=400)
-    res_df = res.to_dataframe().query("Condition == '1.0'")
+    res_df = res.to_dataframe().query("Condition == 'e1p0'")
+    assert len(res_df)
     res_df["theta_uM"] = res_df["theta"] * 1e6
 
     # Create ROIs
@@ -224,7 +232,7 @@ def test_results_glm_export_dataframe_region_of_interest_weighted():
     rois["B"] = [1, 3, 5, 7, 9, 11]
     rois["C"] = [6, 7, 8, 9]
 
-    df_uw = res.to_dataframe_region_of_interest(rois, "1.0", weighted=False)
+    df_uw = res.to_dataframe_region_of_interest(rois, "e1p0", weighted=False)
     assert df_uw.shape == (4, 9)
     assert df_uw.Weighted[0] == "Equal"
     thetas = np.array(res_df.theta)
@@ -232,7 +240,7 @@ def test_results_glm_export_dataframe_region_of_interest_weighted():
     assert np.isclose(df_uw.query("ROI == 'A'").theta,
                       thetas[rois["A"]].mean())
 
-    df_w = res.to_dataframe_region_of_interest(rois, "1.0", weighted=True)
+    df_w = res.to_dataframe_region_of_interest(rois, "e1p0", weighted=True)
     assert df_w.shape == (4, 9)
     assert df_w.Weighted[0] == "Inverse standard error"
     # weighted option should result in larger response
@@ -245,7 +253,7 @@ def test_results_glm_export_dataframe_region_of_interest_weighted():
     weights["B"] = [0.6, 0.1, 0.3, 10000.4, 10000.4, 0.4]
     weights["C"] = [16, 7, 8, 9]
 
-    df = res.to_dataframe_region_of_interest(rois, "1.0", weighted=weights)
+    df = res.to_dataframe_region_of_interest(rois, "e1p0", weighted=weights)
     assert df.shape == (4, 9)
     assert df.Weighted[0] == "Custom"
     assert np.isclose(thetas[0], df.theta[0], atol=0.1e-6)
@@ -255,18 +263,18 @@ def test_results_glm_export_dataframe_region_of_interest_weighted():
 
     with pytest.raises(ValueError, match='must be positive'):
         weights["C"] = [16, 7, -8, 9]
-        _ = res.to_dataframe_region_of_interest(rois, "1.0",
+        _ = res.to_dataframe_region_of_interest(rois, "e1p0",
                                                 weighted=weights)
 
     with pytest.raises(ValueError, match='length of the keys'):
         weights["C"] = [16, 7]
-        _ = res.to_dataframe_region_of_interest(rois, "1.0",
+        _ = res.to_dataframe_region_of_interest(rois, "e1p0",
                                                 weighted=weights)
 
     with pytest.raises(KeyError, match='Keys of group_by and weighted'):
         bad_weights = dict()
         bad_weights["Z"] = [0, 2, 4]
-        _ = res.to_dataframe_region_of_interest(rois, "1.0",
+        _ = res.to_dataframe_region_of_interest(rois, "e1p0",
                                                 weighted=bad_weights)
 
 
diff --git a/mne_nirs/statistics/tests/test_statistics.py b/mne_nirs/statistics/tests/test_statistics.py
index 77b87a157..590412d24 100644
--- a/mne_nirs/statistics/tests/test_statistics.py
+++ b/mne_nirs/statistics/tests/test_statistics.py
@@ -16,7 +16,6 @@
 iir_filter = [1., -0.58853134, -0.29575669, -0.52246482, 0.38735476, 0.024286]
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_run_GLM():
     raw = simulate_nirs_raw(sig_dur=200, stim_dur=5.)
     design_matrix = make_first_level_design_matrix(raw, stim_dur=5.,
@@ -25,7 +24,8 @@ def test_run_GLM():
     glm_estimates = run_glm(raw, design_matrix)
 
     # Test backwards compatibility
-    old_res = run_GLM(raw, design_matrix)
+    with pytest.deprecated_call(match='more comprehensive'):
+        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()
diff --git a/mne_nirs/statistics/tests/test_statsmodels.py b/mne_nirs/statistics/tests/test_statsmodels.py
index 16d73e6ad..d67db43ef 100644
--- a/mne_nirs/statistics/tests/test_statsmodels.py
+++ b/mne_nirs/statistics/tests/test_statsmodels.py
@@ -9,13 +9,17 @@
 import pandas as pd
 import statsmodels.formula.api as smf
 
-from ...simulation import simulate_nirs_raw
-from ...experimental_design import make_first_level_design_matrix
-from ...statistics import run_glm, statsmodels_to_results
+from mne.utils import check_version
 
+from mne_nirs.simulation import simulate_nirs_raw
+from mne_nirs.experimental_design import make_first_level_design_matrix
+from mne_nirs.statistics import run_glm, statsmodels_to_results
 
+
+@pytest.mark.skipif(not check_version('lxml'), reason='Requires lxml')
 @pytest.mark.parametrize('func', ('mixedlm', 'ols', 'rlm'))
 @pytest.mark.filterwarnings('ignore:.*optimization.*:')
+@pytest.mark.filterwarnings('ignore:.*unknown kwargs.*:')
 @pytest.mark.filterwarnings('ignore:.*on the boundary.*:')
 @pytest.mark.filterwarnings('ignore:.*The Hessian matrix at the estimated.*:')
 def test_statsmodel_to_df(func):
@@ -36,7 +40,7 @@ def test_statsmodel_to_df(func):
         with pytest.warns(RuntimeWarning, match='Non standard source detect'):
             cha = glm_est.to_dataframe()
         cha["ID"] = '%02d' % n
-        df_cha = df_cha.append(cha)
+        df_cha = pd.concat([df_cha, cha], ignore_index=True)
     df_cha["theta"] = df_cha["theta"] * 1.0e6
     roi_model = func("theta ~ -1 + Condition", df_cha,
                      groups=df_cha["ID"]).fit()
diff --git a/mne_nirs/tests/test_examples.py b/mne_nirs/tests/test_examples.py
index 84f2b0e81..3f4b1f0ad 100644
--- a/mne_nirs/tests/test_examples.py
+++ b/mne_nirs/tests/test_examples.py
@@ -7,6 +7,8 @@
 import os
 import pytest
 
+from mne.utils import check_version
+
 
 def examples_path():
 
@@ -23,29 +25,34 @@ def examples_path():
 
 
 def run_script_and_check(test_file_path):
-    import matplotlib as mpl
-    mpl.use("Agg")  # Useful when testing locally as not to block progress
-    return exec(open(test_file_path).read(), locals(), locals())
+    with open(test_file_path) as fid:
+        return exec(fid.read(), locals(), locals())
+
+
+requires_lets_plot = pytest.mark.skipif(
+    not check_version('lets_plot'), reason='Requires lets-plot')
 
 
+@pytest.mark.filterwarnings('ignore:No bad channels to interpolate.*:')
 @pytest.mark.examples
-@pytest.mark.parametrize('fname', (["plot_01_data_io.py",
-                                    "plot_05_datasets.py",
-                                    "plot_10_hrf_simulation.py",
-                                    "plot_11_hrf_measured.py",
-                                    "plot_12_group_glm.py",
-                                    "plot_13_fir_glm.py",
-                                    "plot_14_glm_components.py",
-                                    "plot_15_waveform.py",
-                                    "plot_16_waveform_group.py",
-                                    "plot_19_snirf.py",
-                                    "plot_20_enhance.py",
-                                    "plot_21_artifacts.py",
-                                    "plot_22_quality.py",
-                                    "plot_30_frequency.py",
-                                    "plot_40_mayer.py",
-                                    "plot_80_save_read_glm.py",
-                                    "plot_99_bad.py"]))
+@pytest.mark.parametrize('fname', ([
+    "plot_01_data_io.py",
+    "plot_05_datasets.py",
+    "plot_10_hrf_simulation.py",
+    "plot_11_hrf_measured.py",
+    pytest.param("plot_12_group_glm.py", marks=requires_lets_plot),
+    "plot_13_fir_glm.py",
+    "plot_14_glm_components.py",
+    "plot_15_waveform.py",
+    pytest.param("plot_16_waveform_group.py", marks=requires_lets_plot),
+    "plot_19_snirf.py",
+    "plot_20_enhance.py",
+    "plot_21_artifacts.py",
+    "plot_22_quality.py",
+    "plot_30_frequency.py",
+    "plot_40_mayer.py",
+    "plot_80_save_read_glm.py",
+    "plot_99_bad.py"]))
 def test_examples(fname):
     test_file_path = examples_path() + fname
     run_script_and_check(test_file_path)
diff --git a/mne_nirs/utils/tests/test_io.py b/mne_nirs/utils/tests/test_io.py
index fd0b5ebd9..786f1b770 100644
--- a/mne_nirs/utils/tests/test_io.py
+++ b/mne_nirs/utils/tests/test_io.py
@@ -14,13 +14,16 @@
 from mne_nirs.statistics import run_glm
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
+@pytest.mark.filterwarnings('ignore:.*more comprehensive.*:')
 def test_io():
     num_chans = 6
     fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
     fnirs_raw_dir = os.path.join(fnirs_data_folder, 'Participant-1')
     raw_intensity = mne.io.read_raw_nirx(fnirs_raw_dir).load_data()
     raw_intensity.resample(0.2)
+    raw_intensity.annotations.description[:] = [
+        'e' + d.replace('.', 'p')
+        for d in raw_intensity.annotations.description]
     raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
     raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)
     raw_haemo = mne_nirs.channels.get_long_channels(raw_haemo)
@@ -52,7 +55,7 @@ def test_io():
     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['2.0'] - basic_conts['3.0']
+    contrast_LvR = basic_conts['e2p0'] - basic_conts['e3p0']
 
     contrast = mne_nirs.statistics.compute_contrast(glm_est.data, contrast_LvR)
     df = glm_to_tidy(raw_haemo, contrast, design_matrix)
diff --git a/mne_nirs/visualisation/tests/test_visualisation.py b/mne_nirs/visualisation/tests/test_visualisation.py
index 62830128a..1929a66ad 100644
--- a/mne_nirs/visualisation/tests/test_visualisation.py
+++ b/mne_nirs/visualisation/tests/test_visualisation.py
@@ -49,7 +49,7 @@ def test_plot_nirs_source_detector_pyvista():
         verbose=True)
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
+@pytest.mark.filterwarnings('ignore:"plot_glm_topo" has been deprecated.*:')
 def test_run_plot_GLM_topo():
     raw_intensity = _load_dataset()
     raw_intensity.crop(450, 600)  # Keep the test fast
@@ -102,7 +102,6 @@ def test_run_plot_GLM_topo():
         plot_glm_topo(raw_haemo, glm_estimates_subset, design_matrix)
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_run_plot_GLM_contrast_topo():
     raw_intensity = _load_dataset()
     raw_intensity.crop(450, 600)  # Keep the test fast
@@ -117,12 +116,15 @@ def test_run_plot_GLM_contrast_topo():
     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.data, contrast_LvR)
-    fig = mne_nirs.visualisation.plot_glm_contrast_topo(raw_haemo, contrast)
+    with pytest.deprecated_call(match='comprehensive GLM'):
+        contrast = mne_nirs.statistics.compute_contrast(
+            glm_est.data, contrast_LvR)
+    with pytest.deprecated_call(match='comprehensive GLM'):
+        fig = mne_nirs.visualisation.plot_glm_contrast_topo(
+            raw_haemo, contrast)
     assert len(fig.axes) == 3
 
 
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_run_plot_GLM_contrast_topo_single_chroma():
     raw_intensity = _load_dataset()
     raw_intensity.crop(450, 600)  # Keep the test fast
@@ -138,8 +140,12 @@ def test_run_plot_GLM_contrast_topo_single_chroma():
     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.data, contrast_LvR)
-    fig = mne_nirs.visualisation.plot_glm_contrast_topo(raw_haemo, contrast)
+    with pytest.deprecated_call(match='comprehensive GLM'):
+        contrast = mne_nirs.statistics.compute_contrast(
+            glm_est.data, contrast_LvR)
+    with pytest.deprecated_call(match='comprehensive GLM'):
+        fig = mne_nirs.visualisation.plot_glm_contrast_topo(
+            raw_haemo, contrast)
     assert len(fig.axes) == 2
 
 
@@ -150,7 +156,6 @@ def test_fig_from_axes():
 
 
 @requires_pyvista()
-@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
 def test_run_plot_GLM_projection():
     mne.viz.set_3d_backend('pyvista')
     raw_intensity = _load_dataset()