[ENH,FIX] optimally combine temporal + dispersion derivatives for beta estimation (#296)

* add function _calc_beta_map

- the function implements calhoun 2004 beta estimate combination with temporal derivatives
- also minor switch from Nift1Image to Nifti2image

* test various hrf_models on the _calc_beta_map function

Author: jdkent

src/nibetaseries/interfaces/nistats.py

@@ -104,7 +104,7 @@ def _run_interface(self, runtime):
                         beta_maps[new_delay_ttype] = [beta_map]
             else:
                 # calculate the beta map
-                beta_map = model.compute_contrast(trial_type, output_type='effect_size')
+                beta_map = _calc_beta_map(model, trial_type, self.inputs.hrf_model)
                 design_matrix_collector[trial_idx] = model.design_matrices_[0]
                 # assign beta map to appropriate list
                 if trial_type in beta_maps:
@@ -122,7 +122,7 @@ def _run_interface(self, runtime):
         ave_residual = residuals / (trial_idx + 1)
         # make residual nifti image
         residual_file = os.path.join(runtime.cwd, 'desc-residuals_bold.nii.gz')
-        nib.Nifti1Image(
+        nib.Nifti2Image(
             ave_residual,
             model.residuals[0].affine,
             model.residuals[0].header,
@@ -213,7 +213,7 @@ def _run_interface(self, runtime):
             t_type = lsa_df.loc[i_trial, 'original_trial_type']
 
             # calculate the beta map
-            beta_map = model.compute_contrast(t_name, output_type='effect_size')
+            beta_map = _calc_beta_map(model, t_name, self.inputs.hrf_model)
 
             # assign beta map to appropriate list
             if t_type in beta_maps:
@@ -344,3 +344,51 @@ def _select_confounds(confounds_file, selected_confounds):
         msg = "The selected confounds contain nans: {conf}".format(conf=expanded_confounds)
         raise ValueError(msg)
     return desired_confounds
+
+
+def _calc_beta_map(model, trial_type, hrf_model):
+    """
+    Calculates the beta estimates for every voxel from
+    a nistats model
+
+    Parameters
+    ----------
+    model : nistats.first_level_model.FirstLevelModel
+        a fit model of the first level results
+
+    Returns
+    -------
+    beta_map : nibabel.nifti2.Nifti2Image
+        nifti image containing voxelwise beta estimates
+    """
+    import numpy as np
+
+    # calculate the beta map
+    beta_map_list = []
+    beta_map_base = model.compute_contrast(trial_type, output_type='effect_size')
+    beta_map_list.append(beta_map_base.get_fdata())
+    beta_sign = np.where(beta_map_list[0] < 0, -1, 1)
+    if 'derivative' in hrf_model:
+        td_contrast = '_'.join([trial_type, 'derivative'])
+        beta_map_list.append(
+            model.compute_contrast(
+                td_contrast, output_type='effect_size').get_fdata())
+    if 'dispersion' in hrf_model:
+        dd_contrast = '_'.join([trial_type, 'dispersion'])
+        beta_map_list.append(
+            model.compute_contrast(
+                dd_contrast, output_type='effect_size').get_fdata())
+
+    if len(beta_map_list) == 1:
+        beta_map = beta_map_base
+    else:
+        beta_map_array = beta_sign * \
+            np.sqrt(
+                np.sum(
+                    np.array([np.power(c, 2) for c in beta_map_list]), axis=0))
+        beta_map = nib.Nifti2Image(
+            beta_map_array,
+            beta_map_base.affine,
+            beta_map_base.header)
+
+    return beta_map

src/nibetaseries/interfaces/tests/test_nistats.py

@@ -6,17 +6,27 @@
 import nibabel as nib
 from nilearn.image import load_img
 import pandas as pd
+from nistats import first_level_model
 import pytest
 
 
 from ..nistats import (LSSBetaSeries, LSABetaSeries,
                        _lss_events_iterator, _lsa_events_converter,
-                       _select_confounds)
+                       _select_confounds,
+                       _calc_beta_map)
 
 
-@pytest.mark.parametrize("use_nibabel", [(True), (False)])
+@pytest.mark.parametrize(
+    "use_nibabel,hrf_model",
+    [
+        (True, 'spm'),
+        (False, 'spm + derivative'),
+        (False, 'spm + derivative + dispersion')
+    ]
+)
 def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
-                         confounds_file, brainmask_file, use_nibabel):
+                         confounds_file, brainmask_file, use_nibabel,
+                         hrf_model):
     """Test lss interface with nibabel nifti images
     """
     if use_nibabel:
@@ -27,7 +37,6 @@ def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
         mask_file = str(brainmask_file)
 
     selected_confounds = ['white_matter', 'csf']
-    hrf_model = 'spm'
     with open(str(sub_metadata), 'r') as md:
         bold_metadata = json.load(md)
 
@@ -274,3 +283,39 @@ def test_select_confounds(confounds_file, selected_confounds, nan_confounds,
             vals = confounds_df[nan_c].values
             expected_result = np.nanmean(vals[vals != 0])
             assert res_df[nan_c][0] == expected_result
+
+
+@pytest.mark.parametrize(
+    "hrf_model",
+    [
+        ("glover"),
+        ("glover + derivative"),
+        ("glover + derivative + dispersion"),
+    ],
+)
+def test_calc_beta_map(sub_metadata, preproc_file, sub_events,
+                       confounds_file, brainmask_file, hrf_model):
+
+    model = first_level_model.FirstLevelModel(
+            t_r=2,
+            slice_time_ref=0,
+            hrf_model=hrf_model,
+            mask=str(brainmask_file),
+            smoothing_fwhm=0,
+            signal_scaling=False,
+            high_pass=0.008,
+            drift_model='cosine',
+            verbose=1,
+            minimize_memory=False,
+        )
+
+    lsa_df = _lsa_events_converter(str(sub_events))
+
+    model.fit(str(preproc_file), events=lsa_df)
+
+    i_trial = lsa_df.index[0]
+    t_name = lsa_df.loc[i_trial, 'trial_type']
+
+    beta_map = _calc_beta_map(model, t_name, hrf_model)
+
+    assert beta_map.shape == nib.load(str(brainmask_file)).shape