[ENH] add --normalize-betas option (#299)

* refactor _calc_beta_maps to allow for tstats output

* test tstat/norm_beta functionality

* propagate norm_betas through the workflow and add tests

Author: jdkent

docs/workflows.rst

@@ -23,6 +23,7 @@ Participant Workflow
         hrf_model='glover',
         high_pass=0.008,
         name='subtest',
+        norm_betas=False,
         output_dir='.',
         preproc_img_list=[''],
         selected_confounds=[''],
@@ -54,6 +55,7 @@ Least Squares- Separate (LSS)
         fir_delays=None,
         hrf_model='glover',
         high_pass=0.008,
+        norm_betas=False,
         smoothing_kernel=0.0,
         signal_scaling=0,
         selected_confounds=[''])
@@ -76,6 +78,7 @@ Finite BOLD Response- Separate (FS)
         fir_delays=[0, 1, 2, 3, 4],
         hrf_model='fir',
         high_pass=0.008,
+        norm_betas=False,
         smoothing_kernel=0.0,
         signal_scaling=0,
         selected_confounds=[''])
@@ -99,6 +102,7 @@ Least Squares- All (LSA)
         fir_delays=None,
         hrf_model='glover',
         high_pass=0.008,
+        norm_betas=False,
         smoothing_kernel=0.0,
         signal_scaling=0,
         selected_confounds=[''])

src/nibetaseries/cli/run.py

@@ -106,6 +106,9 @@ def get_parser():
                            help='do not scale every voxel with respect to time meaning'
                                 ' beta estimates will be in the same units as the bold'
                                 ' signal')
+    proc_opts.add_argument('--normalize-betas', action='store_true', default=False,
+                           help='beta estimates will be divided by the square root '
+                                'of their variance')
 
     # Image Selection options
     bids_opts = parser.add_argument_group('Options for selecting images')
@@ -260,6 +263,7 @@ def main():
             fir_delays=opts.fir_delays,
             hrf_model=opts.hrf_model,
             high_pass=opts.high_pass,
+            norm_betas=opts.normalize_betas,
             output_dir=output_dir,
             run_label=opts.run_label,
             signal_scaling=signal_scaling,

src/nibetaseries/cli/tests/test_run.py

@@ -36,15 +36,18 @@ def test_conditional_arguments(monkeypatch):
         get_parser().parse_args(no_img)
 
 
-@pytest.mark.parametrize("use_atlas,estimator,fir_delays,hrf_model,part_label,use_signal_scaling",
-                         [(True, 'lsa', None, 'spm', '01', True),
-                          (False, 'lss', None, 'spm', 'sub-01', False),
-                          (True, 'lss', [0, 1, 2, 3, 4], 'fir', None, False)])
+@pytest.mark.parametrize(
+    "use_atlas,estimator,fir_delays,hrf_model,part_label,use_signal_scaling,norm_betas",
+    [
+        (True, 'lsa', None, 'spm', '01', True, True),
+        (False, 'lss', None, 'spm', 'sub-01', False, False),
+        (True, 'lss', [0, 1, 2, 3, 4], 'fir', None, False, True)
+    ])
 def test_nibs(
         bids_dir, deriv_dir, sub_fmriprep, sub_metadata, bold_file, preproc_file,
         sub_events, confounds_file, brainmask_file, atlas_file, atlas_lut,
         estimator, fir_delays, hrf_model, monkeypatch, part_label, use_atlas,
-        use_signal_scaling):
+        use_signal_scaling, norm_betas):
     import sys
     bids_dir = str(bids_dir)
     out_dir = os.path.join(bids_dir, 'derivatives')
@@ -65,6 +68,8 @@ def test_nibs(
             out_dir,
             "participant",
             "-c", ".*derivative.*"]
+    if norm_betas:
+        args.extend(['--normalize-betas'])
     if use_signal_scaling:
         args.extend(["--no-signal-scaling"])
     if use_atlas:

src/nibetaseries/interfaces/nistats.py

@@ -39,6 +39,7 @@ class LSSBetaSeriesInputSpec(BaseInterfaceInputSpec):
     fir_delays = traits.Either(None, traits.List(traits.Int),
                                desc="FIR delays (in scans)",
                                default=None, usedefault=True)
+    return_tstat = traits.Bool(desc="use the T-statistic instead of the raw beta estimates")
 
 
 class LSSBetaSeriesOutputSpec(TraitedSpec):
@@ -99,15 +100,20 @@ def _run_interface(self, runtime):
                     delay_ttype = trial_type+'_delay_{}'.format(delay)
                     new_delay_ttype = delay_ttype.replace('_delay_{}'.format(delay),
                                                           'Delay{}Vol'.format(delay))
-                    beta_map = model.compute_contrast(
-                        delay_ttype, output_type='effect_size')
+                    beta_map = _calc_beta_map(model,
+                                              delay_ttype,
+                                              self.inputs.hrf_model,
+                                              self.inputs.return_tstat)
                     if new_delay_ttype in beta_maps:
                         beta_maps[new_delay_ttype].append(beta_map)
                     else:
                         beta_maps[new_delay_ttype] = [beta_map]
             else:
                 # calculate the beta map
-                beta_map = _calc_beta_map(model, trial_type, self.inputs.hrf_model)
+                beta_map = _calc_beta_map(model,
+                                          trial_type,
+                                          self.inputs.hrf_model,
+                                          self.inputs.return_tstat)
                 design_matrix_collector[trial_idx] = model.design_matrices_[0]
                 # assign beta map to appropriate list
                 if trial_type in beta_maps:
@@ -168,6 +174,7 @@ class LSABetaSeriesInputSpec(BaseInterfaceInputSpec):
     smoothing_kernel = traits.Either(None, traits.Float(),
                                      desc="full wide half max smoothing kernel")
     high_pass = traits.Float(0.0078125, desc="the high pass filter (Hz)")
+    return_tstat = traits.Bool(desc="use the T-statistic instead of the raw beta estimates")
 
 
 class LSABetaSeriesOutputSpec(TraitedSpec):
@@ -219,7 +226,10 @@ def _run_interface(self, runtime):
             t_type = lsa_df.loc[i_trial, 'original_trial_type']
 
             # calculate the beta map
-            beta_map = _calc_beta_map(model, t_name, self.inputs.hrf_model)
+            beta_map = _calc_beta_map(model,
+                                      t_name,
+                                      self.inputs.hrf_model,
+                                      self.inputs.return_tstat)
 
             # assign beta map to appropriate list
             if t_type in beta_maps:
@@ -352,7 +362,7 @@ def _select_confounds(confounds_file, selected_confounds):
     return desired_confounds
 
 
-def _calc_beta_map(model, trial_type, hrf_model):
+def _calc_beta_map(model, trial_type, hrf_model, tstat):
     """
     Calculates the beta estimates for every voxel from
     a nistats model
@@ -361,6 +371,12 @@ def _calc_beta_map(model, trial_type, hrf_model):
     ----------
     model : nistats.first_level_model.FirstLevelModel
         a fit model of the first level results
+    trial_type : str
+        the trial to create the beta estimate
+    hrf_model : str
+        the hemondynamic response function used to fit the model
+    tstat : bool
+        return the t-statistic for the betas instead of the raw estimates
 
     Returns
     -------
@@ -369,32 +385,68 @@ def _calc_beta_map(model, trial_type, hrf_model):
     """
     import numpy as np
 
+    # make it so we do not divide by zero
+    TINY = 1e-50
+    raw_beta_map = _estimate_map(model, trial_type, hrf_model, 'effect_size')
+    if tstat:
+        var_map = _estimate_map(model, trial_type, hrf_model, 'effect_variance')
+        tstat_array = raw_beta_map.get_fdata() / np.sqrt(np.maximum(var_map.get_fdata(), TINY))
+        return nib.Nifti2Image(tstat_array, raw_beta_map.affine, raw_beta_map.header)
+    else:
+        return raw_beta_map
+
+
+def _estimate_map(model, trial_type, hrf_model, output_type):
+    """
+    Calculates model output for every voxel from
+    a nistats model
+
+    Parameters
+    ----------
+    model : nistats.first_level_model.FirstLevelModel
+        a fit model of the first level results
+    trial_type : str
+        the trial to create the beta estimate
+    hrf_model : str
+        the hemondynamic response function used to fit the model
+    output_type : str
+        Type of the output map.
+        Can be ‘z_score’, ‘stat’, ‘p_value’, ‘effect_size’, or ‘effect_variance’
+
+    Returns
+    -------
+    map_img : nibabel.nifti2.Nifti2Image
+        nifti image containing voxelwise output_type 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)
+    map_list = []
+    map_base = model.compute_contrast(trial_type, output_type=output_type)
+    map_list.append(map_base.get_fdata())
+    sign = np.where(map_list[0] < 0, -1, 1)
     if 'derivative' in hrf_model:
         td_contrast = '_'.join([trial_type, 'derivative'])
-        beta_map_list.append(
+        map_list.append(
             model.compute_contrast(
-                td_contrast, output_type='effect_size').get_fdata())
+                td_contrast, output_type=output_type).get_fdata())
     if 'dispersion' in hrf_model:
         dd_contrast = '_'.join([trial_type, 'dispersion'])
-        beta_map_list.append(
+        map_list.append(
             model.compute_contrast(
-                dd_contrast, output_type='effect_size').get_fdata())
+                dd_contrast, output_type=output_type).get_fdata())
 
-    if len(beta_map_list) == 1:
-        beta_map = beta_map_base
+    if len(map_list) == 1:
+        map_img = map_base
     else:
-        beta_map_array = beta_sign * \
+        map_array = 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)
+                    np.array([np.power(c, 2) for c in map_list]), axis=0))
+
+        map_img = nib.Nifti2Image(
+            map_array,
+            map_base.affine,
+            map_base.header)
 
-    return beta_map
+    return map_img

src/nibetaseries/interfaces/tests/test_nistats.py

@@ -17,16 +17,16 @@
 
 
 @pytest.mark.parametrize(
-    "use_nibabel,hrf_model",
+    "use_nibabel,hrf_model,return_tstat",
     [
-        (True, 'spm'),
-        (False, 'spm + derivative'),
-        (False, 'spm + derivative + dispersion')
+        (True, 'spm', True),
+        (False, 'spm + derivative', False),
+        (False, 'spm + derivative + dispersion', True)
     ]
 )
 def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
                          confounds_file, brainmask_file, use_nibabel,
-                         hrf_model):
+                         hrf_model, return_tstat):
     """Test lss interface with nibabel nifti images
     """
     if use_nibabel:
@@ -48,6 +48,7 @@ def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
                                 selected_confounds=selected_confounds,
                                 signal_scaling=0,
                                 hrf_model=hrf_model,
+                                return_tstat=return_tstat,
                                 smoothing_kernel=None,
                                 high_pass=0.008)
     res = beta_series.run()
@@ -102,6 +103,7 @@ def test_fs_beta_series(sub_metadata, preproc_file, sub_events,
                                 signal_scaling=0,
                                 hrf_model=hrf_model,
                                 fir_delays=fir_delays,
+                                return_tstat=False,
                                 smoothing_kernel=None,
                                 high_pass=0.008)
     res = beta_series.run()
@@ -131,9 +133,17 @@ def test_fs_beta_series(sub_metadata, preproc_file, sub_events,
     os.remove(res.outputs.residual)
 
 
-@pytest.mark.parametrize("use_nibabel", [(True), (False)])
+@pytest.mark.parametrize(
+    "use_nibabel,hrf_model,return_tstat",
+    [
+        (True, 'spm', True),
+        (False, 'spm + derivative', False),
+        (False, 'spm + derivative + dispersion', True)
+    ]
+)
 def test_lsa_beta_series(sub_metadata, preproc_file, sub_events,
-                         confounds_file, brainmask_file, use_nibabel):
+                         confounds_file, brainmask_file, use_nibabel,
+                         hrf_model, return_tstat):
     if use_nibabel:
         bold_file = nib.load(str(preproc_file))
         mask_file = nib.load(str(brainmask_file))
@@ -142,7 +152,7 @@ def test_lsa_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)
 
@@ -154,6 +164,7 @@ def test_lsa_beta_series(sub_metadata, preproc_file, sub_events,
                                 signal_scaling=0,
                                 selected_confounds=selected_confounds,
                                 hrf_model=hrf_model,
+                                return_tstat=return_tstat,
                                 smoothing_kernel=None,
                                 high_pass=0.008)
     res = beta_series.run()
@@ -289,15 +300,16 @@ def test_select_confounds(confounds_file, selected_confounds, nan_confounds,
 
 
 @pytest.mark.parametrize(
-    "hrf_model",
+    "hrf_model,return_tstat",
     [
-        ("glover"),
-        ("glover + derivative"),
-        ("glover + derivative + dispersion"),
+        ("glover", True),
+        ("glover + derivative", False),
+        ("glover + derivative + dispersion", True),
     ],
 )
 def test_calc_beta_map(sub_metadata, preproc_file, sub_events,
-                       confounds_file, brainmask_file, hrf_model):
+                       confounds_file, brainmask_file, hrf_model,
+                       return_tstat):
 
     model = first_level_model.FirstLevelModel(
             t_r=2,
@@ -319,6 +331,6 @@ def test_calc_beta_map(sub_metadata, preproc_file, sub_events,
     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)
+    beta_map = _calc_beta_map(model, t_name, hrf_model, return_tstat)
 
     assert beta_map.shape == nib.load(str(brainmask_file)).shape