rfx_center and rfx_sigma input parameters to better control RFX

Author: EtienneCmb

frites/stats/stats_param.py

@@ -8,8 +8,25 @@
 logger = logging.getLogger("frites")
 
 
-RECENTER = dict(mean=np.mean, median=np.median,
-                trimmed=lambda x, axis=0: trim_mean(x, .2, axis=axis))
+def _trimmed(x, prop=.2, axis=0, keepdims=False):
+    trm = trim_mean(x, prop, axis=axis)
+    if keepdims:
+        ax = [slice(None)] * x.ndim
+        ax[axis] = np.newaxis
+        trm = trm[tuple(ax)]
+    return trm
+
+
+RECENTER = {
+    'mean': np.mean, 'median': np.median, 'zscore': np.mean,
+    'trimmed': _trimmed
+}
+
+
+def _recenter(x, fcn_mean, zscore, axis=-1):
+    """Recentering function."""
+    _std = x.std(axis=axis, keepdims=True) if zscore else 1.
+    return (x - fcn_mean(x, axis=axis, keepdims=True)) / _std
 
 
 def ttest_1samp(x, pop_mean, axis=0, implementation='mne', sigma=0.001, **kw):
@@ -50,7 +67,7 @@ def fcn(x, pop_mean, axis):  # noqa
     return fcn(x, pop_mean, axis)
 
 
-def rfx_ttest(mi, mi_p, center=False, zscore=False, ttested=False):
+def rfx_ttest(mi, mi_p, center=False, sigma=0.001, ttested=False):
     """Perform the t-test across subjects.
 
     Parameters
@@ -61,12 +78,12 @@ def rfx_ttest(mi, mi_p, center=False, zscore=False, ttested=False):
     mi_p : array_like
         A list of array of permuted mutual information of shape
         (n_perm, n_suj, n_times). If `ttested` is True, n_suj shoud be 1.
-    center : {'mean', 'median', 'trimmed'} | False
-        If True, substract the mean of the surrogates to the true and permuted
-        mi. The median or the 20% trimmed mean can also be removed
-        :cite:`wilcox2018guide`
-    zscore : bool | False
-        Apply z-score normalization to the true and permuted mutual information
+    sigma : float | 0.001
+        Hat adjustment method, a value of 1e-3 may be a reasonable choice
+    center : {False, 'mean', 'median', 'trimmed', 'zscore'}
+        Re-center the time-series of effect arround 0 before computing the
+        t-test. This parameters can be useful in case of a different number
+        of data per brain region.
     ttested : bool | False
         Specify if the inputs have already been t-tested
 
@@ -93,14 +110,17 @@ def rfx_ttest(mi, mi_p, center=False, zscore=False, ttested=False):
     n_roi = len(mi_p)
 
     # remove the mean / median / trimmed
+    zscore = center == 'zscore'
     if center in RECENTER.keys():
-        logger.info(f"    RFX recenter distributions (center={center}, "
-                    f"z-score={zscore})")
-        for k in range(len(mi)):
-            _med = RECENTER[center](mi_p[k], axis=0)
-            _std = mi_p[k].std(axis=0) if zscore else 1.
-            mi[k] = (mi[k] - _med) / _std
-            mi_p[k] = (mi_p[k] - _med) / _std
+        # get the centering function
+        fcn_mean = RECENTER[center]
+
+        # here, we need to make a copy of the effect sizes to avoid changing
+        # the ouputs
+        mi, mi_p = mi.copy(), mi_p.copy()
+        for k in range(n_roi):
+            mi[k] = _recenter(mi[k], fcn_mean, zscore, axis=-1).copy()
+            mi_p[k] = _recenter(mi_p[k], fcn_mean, zscore, axis=-1).copy()
 
     # get the mean of surrogates (low ram method)
     n_element = np.sum([np.prod(k.shape) for k in mi_p])
@@ -112,12 +132,14 @@ def rfx_ttest(mi, mi_p, center=False, zscore=False, ttested=False):
     that the MNE t-test is going to evaluate one sigma per roi. To fix that,
     we estimate this sigma using the variance of all of the data
     """
-    from frites.config import CONFIG
-    s_hat = CONFIG['TTEST_MNE_SIGMA']
+    s_hat = sigma
 
     # sigma on true mi and permuted mi
-    sigma = s_hat * max([np.var(k, axis=0, ddof=1).max() for k in mi])
-    sigma_p = s_hat * max([np.var(k, axis=1, ddof=1).max() for k in mi_p])
+    if s_hat > 0:
+        sigma = s_hat * max([np.var(k, axis=0, ddof=1).max() for k in mi])
+        sigma_p = s_hat * max([np.var(k, axis=1, ddof=1).max() for k in mi_p])
+    else:
+        sigma = sigma_p = 0.
     logger.debug(f"sigma_true={sigma}; sigma_permuted={sigma_p}")
     kw = dict(implementation='mne', method='absolute', sigma=sigma)
     kw_p = dict(implementation='mne', method='absolute', sigma=sigma_p)

frites/stats/tests/test_stats_param.py

@@ -24,18 +24,10 @@ def test_rfx_ttest(self):
         assert tv.shape == (n_roi, n_times)
         assert tv_p.shape == (n_perm, n_roi, n_times)
         # center
-        for center in [False, 'mean', 'median', 'trimmed']:
+        for center in [False, True, 'mean', 'median', 'trimmed', 'zscore']:
             tv, tv_p, _ = rfx_ttest(x, x_p, center=center)
             assert tv.shape == (n_roi, n_times)
             assert tv_p.shape == (n_perm, n_roi, n_times)
-        # zscore
-        tv, tv_p, _ = rfx_ttest(x, x_p, zscore=True)
-        assert tv.shape == (n_roi, n_times)
-        assert tv_p.shape == (n_perm, n_roi, n_times)
-        # center + zscore
-        tv, tv_p, _ = rfx_ttest(x, x_p, zscore=True, center=True)
-        assert tv.shape == (n_roi, n_times)
-        assert tv_p.shape == (n_perm, n_roi, n_times)
         # t-tested
         tv, tv_p, _ = rfx_ttest(x, x_p, ttested=True)
         assert tv.shape == (n_roi * n_suj, n_times)

frites/workflow/wf_mi.py

@@ -451,7 +451,7 @@ def confidence_interval(self, dataset, ci=95, n_boots=200, rfx_es='mi',
             pbar = ProgressBar(range(n_boots), mesg='Estimating CI')
 
             # get t-test related variables
-            s_hat = CONFIG['TTEST_MNE_SIGMA']
+            s_hat = self._wf_stats.attrs['ttest_sigma']
 
             tt = []
             for n_p in range(n_boots):

frites/workflow/wf_stats.py

@@ -28,7 +28,8 @@ def __init__(self, verbose=None):  # noqa
         logger.info("Definition of a non-parametric statistical workflow")
 
     def fit(self, effect, perms, inference='rfx', mcp='cluster', tail=1,
-            cluster_th=None, cluster_alpha=0.05, ttested=False):
+            cluster_th=None, cluster_alpha=0.05, ttested=False,
+            rfx_sigma=0.001, rfx_center=False):
         """Fit the workflow on true data.
 
         Parameters
@@ -71,6 +72,12 @@ def fit(self, effect, perms, inference='rfx', mcp='cluster', tail=1,
             the 95th percentile of the permutations is used.
         ttested : bool | False
             Specify if the inputs have already been t-tested
+        rfx_sigma : float | 0.001
+            Hat adjustment method, a value of 1e-3 may be a reasonable choice
+        rfx_center : {False, 'mean', 'median', 'trimmed', 'zscore'}
+            Re-center the time-series of effect arround 0 before computing the
+            t-test. This parameters can be useful in case of a different number
+            of data per brain region.
 
         Returns
         -------
@@ -123,11 +130,15 @@ def fit(self, effect, perms, inference='rfx', mcp='cluster', tail=1,
                 rfx_suj = np.min(nb_suj_roi) > 1
                 assert rfx_suj, "For RFX, `n_subjects` should be > 1"
                 # modelise how subjects are distributed
-                es, es_p, pop_mean = rfx_ttest(effect, perms)
-                from frites.config import CONFIG
-                sigma = CONFIG['TTEST_MNE_SIGMA']
-                self.attrs.update(dict(ttest_pop_mean=pop_mean,
-                                       ttest_sigma=sigma))
+                rfx_center = 'mean' if isinstance(
+                    rfx_center, bool) and rfx_center else rfx_center
+                es, es_p, pop_mean = rfx_ttest(
+                    effect, perms, sigma=rfx_sigma, center=rfx_center
+                )
+                self.attrs.update(dict(
+                    ttest_pop_mean=pop_mean, ttest_sigma=rfx_sigma,
+                    ttest_center=rfx_center
+                ))
             tvalues = es
 
         # ---------------------------------------------------------------------