[ENH] add residuals to LSS/LSA interfaces (#294)

* bump nistats version

* add residuals to outputs

* test existance/shape of residuals

* add more comprehensive testing of beta maps

Author: jdkent

setup.cfg

@@ -31,7 +31,7 @@ install_requires =
     nipype ~= 1.4.2
     pybids ~= 0.9.3
     nibabel ~= 2.4.0
-    nistats == 0.0.1b1
+    nistats == 0.0.1b2
     nilearn ~= 0.4.2
     pandas ~= 0.24.0
     numpy

src/nibetaseries/interfaces/nistats.py

@@ -34,12 +34,14 @@ class LSSBetaSeriesInputSpec(BaseInterfaceInputSpec):
                                      desc="full wide half max smoothing kernel")
     high_pass = traits.Float(0.0078125, desc="the high pass filter (Hz)")
     fir_delays = traits.Either(None, traits.List(traits.Int),
-                               desc="FIR delays (in scans)")
+                               desc="FIR delays (in scans)",
+                               default=None, usedefault=True)
 
 
 class LSSBetaSeriesOutputSpec(TraitedSpec):
     beta_maps = OutputMultiPath(File)
     design_matrices = traits.Dict()
+    residual = traits.File(exists=True)
 
 
 class LSSBetaSeries(NistatsBaseInterface, SimpleInterface):
@@ -73,10 +75,12 @@ def _run_interface(self, runtime):
             drift_model='cosine',
             verbose=1,
             fir_delays=self.inputs.fir_delays,
+            minimize_memory=False,
         )
 
         # initialize dictionary to contain trial estimates (betas)
         beta_maps = {}
+        residuals = None
         design_matrix_collector = {}
         for target_trial_df, trial_type, trial_idx in \
                 _lss_events_iterator(self.inputs.events_file):
@@ -108,6 +112,21 @@ def _run_interface(self, runtime):
                 else:
                     beta_maps[trial_type] = [beta_map]
 
+            # add up all the residuals (to be divided later)
+            if residuals is None:
+                residuals = model.residuals[0].get_fdata()
+            else:
+                residuals += model.residuals[0].get_fdata()
+
+        # make an average residual
+        ave_residual = residuals / (trial_idx + 1)
+        # make residual nifti image
+        residual_file = os.path.join(runtime.cwd, 'desc-residuals_bold.nii.gz')
+        nib.Nifti1Image(
+            ave_residual,
+            model.residuals[0].affine,
+            model.residuals[0].header,
+        ).to_filename(residual_file)
         # make a beta series from each beta map list
         beta_series_template = os.path.join(runtime.cwd,
                                             'desc-{trial_type}_betaseries.nii.gz')
@@ -120,6 +139,7 @@ def _run_interface(self, runtime):
 
         self._results['beta_maps'] = beta_series_lst
         self._results['design_matrices'] = design_matrix_collector
+        self._results['residual'] = residual_file
         return runtime
 
 
@@ -147,6 +167,7 @@ class LSABetaSeriesInputSpec(BaseInterfaceInputSpec):
 class LSABetaSeriesOutputSpec(TraitedSpec):
     beta_maps = OutputMultiPath(File)
     design_matrices = traits.List()
+    residual = traits.File(exists=True)
 
 
 class LSABetaSeries(NistatsBaseInterface, SimpleInterface):
@@ -178,7 +199,8 @@ def _run_interface(self, runtime):
             signal_scaling=0,
             high_pass=self.inputs.high_pass,
             drift_model='cosine',
-            verbose=1
+            verbose=1,
+            minimize_memory=False,
         )
 
         # initialize dictionary to contain trial estimates (betas)
@@ -199,6 +221,9 @@ def _run_interface(self, runtime):
             else:
                 beta_maps[t_type] = [beta_map]
 
+        # calculate the residual
+        residual_file = os.path.join(runtime.cwd, 'desc-residuals_bold.nii.gz')
+        model.residuals[0].to_filename(residual_file)
         # make a beta series from each beta map list
         beta_series_template = os.path.join(runtime.cwd,
                                             'desc-{trial_type}_betaseries.nii.gz')
@@ -211,6 +236,7 @@ def _run_interface(self, runtime):
 
         self._results['beta_maps'] = beta_series_lst
         self._results['design_matrices'] = [design_matrix]
+        self._results['residual'] = residual_file
         return runtime
 
 

src/nibetaseries/interfaces/tests/test_nistats.py

@@ -1,7 +1,11 @@
 ''' Testing module for nibetaseries.interfaces.nistats '''
 import os
 import json
+import re
 
+import nibabel as nib
+from nilearn.image import load_img
+import pandas as pd
 import pytest
 
 
@@ -16,7 +20,6 @@ def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
     """Test lss interface with nibabel nifti images
     """
     if use_nibabel:
-        import nibabel as nib
         bold_file = nib.load(str(preproc_file))
         mask_file = nib.load(str(brainmask_file))
     else:
@@ -39,16 +42,35 @@ def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
                                 high_pass=0.008)
     res = beta_series.run()
 
+    events_df = pd.read_csv(str(sub_events), sep='\t')
+    trial_types = events_df['trial_type'].unique()
+
+    # check for the correct number of beta maps
+    assert len(trial_types) == len(res.outputs.beta_maps)
+
+    input_img_dim = load_img(bold_file).shape[:3]
     for beta_map in res.outputs.beta_maps:
+        trial_type = re.search(r'desc-([A-Za-z0-9]+)_', beta_map).groups()[0]
+        # check if correct trial type is made
+        assert trial_type in trial_types
+        # check if output is actually a file
         assert os.path.isfile(beta_map)
+        # check if image dimensions are correct
+        assert input_img_dim == load_img(beta_map).shape[:3]
+        # check if number of trials are correct
+        assert (events_df['trial_type'] == trial_type).sum() == load_img(beta_map).shape[-1]
+        # clean up
         os.remove(beta_map)
 
+    # check residual image
+    assert load_img(bold_file).shape == load_img(res.outputs.residual).shape
+    os.remove(res.outputs.residual)
+
 
 @pytest.mark.parametrize("use_nibabel", [(True), (False)])
 def test_fs_beta_series(sub_metadata, preproc_file, sub_events,
                         confounds_file, brainmask_file, use_nibabel):
     if use_nibabel:
-        import nibabel as nib
         bold_file = nib.load(str(preproc_file))
         mask_file = nib.load(str(brainmask_file))
     else:
@@ -73,16 +95,35 @@ def test_fs_beta_series(sub_metadata, preproc_file, sub_events,
                                 high_pass=0.008)
     res = beta_series.run()
 
+    events_df = pd.read_csv(str(sub_events), sep='\t')
+    trial_types = events_df['trial_type'].unique()
+
+    # check for the correct number of beta maps
+    assert len(trial_types) * len(fir_delays) == len(res.outputs.beta_maps)
+
+    input_img_dim = load_img(bold_file).shape[:3]
     for beta_map in res.outputs.beta_maps:
+        trial_type = re.search(r'desc-([A-Za-z0-9]+)Delay[0-9]+Vol_', beta_map).groups()[0]
+        # check if correct trial type is made
+        assert trial_type in trial_types
+        # check if output is actually a file
         assert os.path.isfile(beta_map)
+        # check if image dimensions are correct
+        assert input_img_dim == load_img(beta_map).shape[:3]
+        # check if number of trials are correct
+        assert (events_df['trial_type'] == trial_type).sum() == load_img(beta_map).shape[-1]
+        # clean up
         os.remove(beta_map)
 
+    # check residual image
+    assert load_img(bold_file).shape == load_img(res.outputs.residual).shape
+    os.remove(res.outputs.residual)
+
 
 @pytest.mark.parametrize("use_nibabel", [(True), (False)])
 def test_lsa_beta_series(sub_metadata, preproc_file, sub_events,
                          confounds_file, brainmask_file, use_nibabel):
     if use_nibabel:
-        import nibabel as nib
         bold_file = nib.load(str(preproc_file))
         mask_file = nib.load(str(brainmask_file))
     else:
@@ -105,10 +146,30 @@ def test_lsa_beta_series(sub_metadata, preproc_file, sub_events,
                                 high_pass=0.008)
     res = beta_series.run()
 
+    events_df = pd.read_csv(str(sub_events), sep='\t')
+    trial_types = events_df['trial_type'].unique()
+
+    # check for the correct number of beta maps
+    assert len(trial_types) == len(res.outputs.beta_maps)
+
+    input_img_dim = load_img(bold_file).shape[:3]
     for beta_map in res.outputs.beta_maps:
+        trial_type = re.search(r'desc-([A-Za-z0-9]+)_', beta_map).groups()[0]
+        # check if correct trial type is made
+        assert trial_type in trial_types
+        # check if output is actually a file
         assert os.path.isfile(beta_map)
+        # check if image dimensions are correct
+        assert input_img_dim == load_img(beta_map).shape[:3]
+        # check if number of trials are correct
+        assert (events_df['trial_type'] == trial_type).sum() == load_img(beta_map).shape[-1]
+        # clean up
         os.remove(beta_map)
 
+    # check residual image
+    assert load_img(bold_file).shape == load_img(res.outputs.residual).shape
+    os.remove(res.outputs.residual)
+
 
 def test_lss_events_iterator(sub_events):
     # all but the first instance of waffle