Repository: HBClab/NiBetaSeries

Issue:
Output normalized betaseries
**Is your feature request related to a problem? Please describe.**
Estimating betaseries can result in poor fits and bizarre parameter estimates. These outlier estimates can then influence subsequent analyses. 

**Describe the solution you'd like**
One approach to dealing with this issue is to normalize the betaseries by the variance of the estimate. Adding this should be a fairly simple extension. In interfaces/nistats.py, could additionally calculate the effect variance and divide the beta_map by the variance:

`beta_var = model.compute_contrast(trial_type, output_type='effect_variance')`
`beta_mapn = beta_map/sqrt(beta_var)`

I'm actually not sure if the beta_map object can be directly used in the calculations like this, but that would be the idea. And these normalized betas could then be saved alongside the betaseries.


Discussion:
jdkent: Hi @michaelmack,

here is the pull request, to normalize the betas (e.g., return t statistics as opposed to raw estimates). If would like you may review it before I merge it into the master branch. 

I will cut a release after one more pull request (`v0.5.0`)
codecov[bot]: # [Codecov](https://codecov.io/gh/HBClab/NiBetaSeries/pull/299?src=pr&el=h1) Report
> Merging [#299](https://codecov.io/gh/HBClab/NiBetaSeries/pull/299?src=pr&el=desc) into [master](https://codecov.io/gh/HBClab/NiBetaSeries/commit/245575d89682ed810e6979018486a7e191e9d0d6?src=pr&el=desc) will **increase** coverage by `0.15%`.
> The diff coverage is `100%`.


```diff
@@            Coverage Diff             @@
##           master     #299      +/-   ##
==========================================
+ Coverage   89.85%   90.01%   +0.15%     
==========================================
  Files          10       10              
  Lines         759      771      +12     
  Branches      115      116       +1     
==========================================
+ Hits          682      694      +12     
  Misses         45       45              
  Partials       32       32
```


Patch:
diff --git a/docs/workflows.rst b/docs/workflows.rst
index da8a1c2..37e7da5 100644
--- a/docs/workflows.rst
+++ b/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=[''])
diff --git a/src/nibetaseries/cli/run.py b/src/nibetaseries/cli/run.py
index 5a0a81b..46edce5 100644
--- a/src/nibetaseries/cli/run.py
+++ b/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,
diff --git a/src/nibetaseries/interfaces/nistats.py b/src/nibetaseries/interfaces/nistats.py
index 225efff..69203c6 100644
--- a/src/nibetaseries/interfaces/nistats.py
+++ b/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 @@ class LSSBetaSeries(NistatsBaseInterface, SimpleInterface):
                     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 @@ class LSABetaSeries(NistatsBaseInterface, SimpleInterface):
             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
diff --git a/src/nibetaseries/workflows/base.py b/src/nibetaseries/workflows/base.py
index b3da68d..2c13d32 100644
--- a/src/nibetaseries/workflows/base.py
+++ b/src/nibetaseries/workflows/base.py
@@ -32,7 +32,7 @@ from ..interfaces.nilearn import CensorVolumes
 def init_nibetaseries_participant_wf(
     estimator, atlas_img, atlas_lut, bids_dir,
     database_path, derivatives_pipeline_dir, exclude_description_label,
-    fir_delays, hrf_model, high_pass, output_dir, run_label,
+    fir_delays, hrf_model, high_pass, norm_betas, output_dir, run_label,
     selected_confounds, session_label, signal_scaling, smoothing_kernel,
     space_label, subject_list, task_label, description_label, work_dir,
         ):
@@ -63,6 +63,8 @@ def init_nibetaseries_participant_wf(
         high_pass : float
             High pass filter to apply to bold (in Hertz).
             Reminder - frequencies _higher_ than this number are kept.
+        norm_betas : Bool
+            If True, beta estimates are divided by the square root of their variance
         output_dir : str
             Directory where derivatives are saved
         run_label : str or None
@@ -172,6 +174,7 @@ It is released under the [CC0]\
             hrf_model=hrf_model,
             high_pass=high_pass,
             name='single_subject' + subject_label + '_wf',
+            norm_betas=norm_betas,
             output_dir=output_dir,
             preproc_img_list=preproc_img_list,
             selected_confounds=selected_confounds,
@@ -194,7 +197,7 @@ It is released under the [CC0]\
 def init_single_subject_wf(
     estimator, atlas_img, atlas_lut, bold_metadata_list, brainmask_list,
     confound_tsv_list, events_tsv_list, fir_delays, hrf_model, high_pass,
-    name, output_dir, preproc_img_list, selected_confounds,
+    name, norm_betas, output_dir, preproc_img_list, selected_confounds,
     signal_scaling, smoothing_kernel,
         ):
     """
@@ -217,6 +220,7 @@ def init_single_subject_wf(
             hrf_model='',
             high_pass='',
             name='subtest',
+            norm_betas=False,
             output_dir='.',
             preproc_img_list=[''],
             selected_confounds=[''],
@@ -245,6 +249,8 @@ def init_single_subject_wf(
         high_pass : float
             High pass filter to apply to bold (in Hertz).
             Reminder - frequencies _higher_ than this number are kept.
+        norm_betas : Bool
+            If True, beta estimates are divided by the square root of their variance
         name : str
             name of the workflow (e.g. ``subject-01_wf``)
         output_dir : str
@@ -313,6 +319,7 @@ def init_single_subject_wf(
                                        fir_delays=fir_delays,
                                        hrf_model=hrf_model,
                                        high_pass=high_pass,
+                                       norm_betas=norm_betas,
                                        selected_confounds=selected_confounds,
                                        signal_scaling=signal_scaling,
                                        smoothing_kernel=smoothing_kernel)
diff --git a/src/nibetaseries/workflows/model.py b/src/nibetaseries/workflows/model.py
index fb13b15..5cf2b42 100644
--- a/src/nibetaseries/workflows/model.py
+++ b/src/nibetaseries/workflows/model.py
@@ -23,6 +23,7 @@ def init_betaseries_wf(name="betaseries_wf",
                        fir_delays=None,
                        hrf_model='glover',
                        high_pass=0.0078125,
+                       norm_betas=False,
                        signal_scaling=0,
                        selected_confounds=None,
                        smoothing_kernel=None,
@@ -55,6 +56,8 @@ def init_betaseries_wf(name="betaseries_wf",
     high_pass : float
         high pass filter to apply to bold (in Hertz).
         Reminder - frequencies _lower_ than this number are kept.
+    norm_betas : Bool
+        If True, beta estimates are divided by the square root of their variance
     selected_confounds : list or None
         the list of confounds to be included in regression.
     signal_scaling : False or 0
@@ -91,6 +94,7 @@ def init_betaseries_wf(name="betaseries_wf",
         fwhm=smoothing_kernel,
         hrf=hrf_model,
         hpf=high_pass,
+        norm_betas=norm_betas,
         selected_confounds=selected_confounds,
         signal_scaling=signal_scaling,
         estimator=estimator,
@@ -111,6 +115,7 @@ def init_betaseries_wf(name="betaseries_wf",
                 selected_confounds=selected_confounds,
                 signal_scaling=signal_scaling,
                 hrf_model=hrf_model,
+                return_tstat=norm_betas,
                 smoothing_kernel=smoothing_kernel,
                 high_pass=high_pass),
             name='betaseries_node')
@@ -119,6 +124,7 @@ def init_betaseries_wf(name="betaseries_wf",
                 selected_confounds=selected_confounds,
                 signal_scaling=signal_scaling,
                 hrf_model=hrf_model,
+                return_tstat=norm_betas,
                 smoothing_kernel=smoothing_kernel,
                 high_pass=high_pass),
             name='betaseries_node')
@@ -141,7 +147,7 @@ def init_betaseries_wf(name="betaseries_wf",
 
 def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
                        selected_confounds, signal_scaling,
-                       estimator, fir_delays=None):
+                       estimator, norm_betas, fir_delays=None):
     from textwrap import dedent
 
     smooth_str = ('smoothed with a Gaussian kernel with a FWHM of {fwhm} mm,'
@@ -153,7 +159,7 @@ def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
                    .format(smooth_str=smooth_str, signal_scale_str=signal_scale_str))
 
     beta_series_tmp = dedent("""
-        After fitting {n_models} model, the parameter estimate (i.e., beta) map
+        After fitting {n_models} model, the{normed} parameter estimate (i.e., beta) map
         associated with the target trial's regressor was retained and concatenated
         into a 4D image with all other trials from the same condition, resulting
         in a set of N 4D images where N refers to the number of conditions in the task.
@@ -182,14 +188,15 @@ def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
             """.format(fir_delays=[str(d) for d in fir_delays]))
 
         beta_series_str = dedent("""\
-            After fitting each model, the parameter estimate (i.e., beta) map associated
+            After fitting each model, the{normed} parameter estimate (i.e., beta) map associated
             with each of the target trial's {n_delays} delay-specific FIR regressors
             was retained and concatenated into delay-specific 4D images with all other
             trials from the same condition, resulting in a set of N * {n_delays} 4D
             images where N refers to the number of conditions in the task.
             The number of volumes in each 4D image represents the number of trials in that
             condition.\
-            """.format(n_delays=len(fir_delays)))
+            """.format(n_delays=len(fir_delays),
+                       normed=' normalized' if norm_betas else ''))
 
     elif estimator == "lss" and hrf != "fir":
         estimator_intro = dedent("""\
@@ -204,7 +211,8 @@ def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
             modeled in their own regressors.\
             """)
 
-        beta_series_str = beta_series_tmp.format(n_models='each')
+        beta_series_str = beta_series_tmp.format(n_models='each',
+                                                 normed=' normalized' if norm_betas else '')
 
     elif estimator == "lsa":
         estimator_intro = dedent("""\
@@ -217,7 +225,8 @@ def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
             was modeled in its own regressor.\
             """)
 
-        beta_series_str = beta_series_tmp.format(n_models='the')
+        beta_series_str = beta_series_tmp.format(n_models='the',
+                                                 normed=' normalized' if norm_betas else '')
 
     else:
         raise ValueError("{est} not a supported estimator".format(est=estimator))


Test patch:
diff --git a/src/nibetaseries/cli/tests/test_run.py b/src/nibetaseries/cli/tests/test_run.py
index c205561..5c8fd75 100644
--- a/src/nibetaseries/cli/tests/test_run.py
+++ b/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:
diff --git a/src/nibetaseries/interfaces/tests/test_nistats.py b/src/nibetaseries/interfaces/tests/test_nistats.py
index b1f7a94..9765894 100644
--- a/src/nibetaseries/interfaces/tests/test_nistats.py
+++ b/src/nibetaseries/interfaces/tests/test_nistats.py
@@ -17,16 +17,16 @@ from ..nistats import (LSSBetaSeries, LSABetaSeries,
 
 
 @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
diff --git a/src/nibetaseries/workflows/tests/test_base.py b/src/nibetaseries/workflows/tests/test_base.py
index c5f2c6f..d6f668f 100644
--- a/src/nibetaseries/workflows/tests/test_base.py
+++ b/src/nibetaseries/workflows/tests/test_base.py
@@ -8,14 +8,18 @@ import pytest
 from ..base import init_nibetaseries_participant_wf, _check_bs_len
 
 
-@pytest.mark.parametrize("estimator,fir_delays,hrf_model,signal_scaling",
-                         [('lsa', None, 'spm', 0),
-                          ('lss', None, 'spm', False),
-                          ('lss', [0, 1, 2, 3, 4], 'fir', False)])
+@pytest.mark.parametrize(
+    "estimator,fir_delays,hrf_model,signal_scaling,norm_betas",
+    [
+        ('lsa', None, 'spm', 0, True),
+        ('lss', None, 'spm', False, False),
+        ('lss', [0, 1, 2, 3, 4], 'fir', False, True),
+    ]
+)
 def test_valid_init_nibetaseries_participant_wf(
         bids_dir, deriv_dir, sub_fmriprep, sub_top_metadata, bold_file, preproc_file,
         sub_events, confounds_file, brainmask_file, atlas_file, atlas_lut, bids_db_file,
-        estimator, fir_delays, hrf_model, signal_scaling):
+        estimator, fir_delays, hrf_model, signal_scaling, norm_betas):
 
     output_dir = op.join(str(bids_dir), 'derivatives', 'atlasCorr')
     work_dir = op.join(str(bids_dir), 'derivatives', 'work')
@@ -39,6 +43,7 @@ def test_valid_init_nibetaseries_participant_wf(
         exclude_description_label=None,
         hrf_model=hrf_model,
         high_pass=0.008,
+        norm_betas=norm_betas,
         output_dir=output_dir,
         run_label=None,
         selected_confounds=['white_matter', 'csf'],
@@ -84,6 +89,7 @@ def test_filters_init_nibetaseries_participant_wf(
             exclude_description_label=None,
             hrf_model='spm',
             high_pass=0.008,
+            norm_betas=False,
             output_dir=output_dir,
             run_label=run_label,
             selected_confounds=['white_matter', 'csf'],