Repository: naplab/naplib-python (language: python)

Issue:
Add "average" parameter to stats.responsive_ttest
<!--
If you'd like to report a bug in naplib-python, fill out the template below and provide any extra information that may be useful / related to your problem.
-->

Would be great to have a parameter in naplab.stats.responsive_ttest to optionally control whether each segment (e.g. 1 second clip before and after onset) is averaged before being concatenated. Averaging may produce more consistent results, though not as good if you only have a few trials.



Discussion / hints:

Fix patch:
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 0778ffb..448ce89 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -28,6 +28,11 @@ Change tags (adopted from `sklearn <https://scikit-learn.org/stable/whats_new/v0
 
 - |API| : you will need to change your code to have the same effect in the future; or a feature will be removed in the future.
 
+Version 1.3.0
+-------------
+- |Feature| : Added the function ``peak_rate`` to ``naplib.features`` which can be used to extract peak rate events from an acoustic stimulus. See documentation for details.
+- |Enhancement| : Added the parameter ``average`` to ``naplib.stats.responsive_ttest`` to enable more robust t-test statistics when there are enough trials. See documentation for details.
+
 Version 1.2.0
 -------------
 - |Enhancement| : Added the parameter ``makedirs`` to ``naplib.io.save`` to enable the automatic creation of directories in the path for the filename provided.
diff --git a/naplib/__init__.py b/naplib/__init__.py
index 179d103..468d309 100644
--- a/naplib/__init__.py
+++ b/naplib/__init__.py
@@ -44,5 +44,5 @@ def set_logging(level: Union[int, str]):
 from .data import Data, join_fields, concat
 import naplib.naplab
 
-__version__ = "1.2.0"
+__version__ = "1.3.0"
 


Test patch:
diff --git a/naplib/stats/responsive_ttest.py b/naplib/stats/responsive_ttest.py
index adb5ed4..bd40b76 100644
--- a/naplib/stats/responsive_ttest.py
+++ b/naplib/stats/responsive_ttest.py
@@ -7,7 +7,7 @@
 from naplib.data import Data
 
 
-def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre_post=[1, 1], alpha=0.05, fdr_method='indep', alternative='two-sided', equal_var=True):
+def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre_post=[1, 1], alpha=0.05, average=False, fdr_method='indep', alternative='two-sided', equal_var=True):
     '''
     Identify responsive electrodes by performing a t-test between response
     values during silence (before stimulus) compared to during speech/sound
@@ -43,7 +43,7 @@ def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre
     sfreq : str | int, default='dataf'
         The sampling frequency of the responses. If a string, specifies field of
         the Data containing the sampling frequency.
-    pre_post : np.ndarray | list
+    pre_post : np.ndarray | list, default=[1, 1]
         List or array of length 2 or 4, giving the time windows (in seconds) to compare
         before and after stimulus onset. If length 2 (such as [x,y]), both numbers must be positive floats, and
         the first number specifies the window [-x, 0) relative to sound onset while the second number
@@ -55,10 +55,13 @@ def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre
         0.8 seconds immediately preceeding stimulus onset are compared to the responses in
         the 1 second immediately following stimulus onset. Sometimes, sound does not onset
         until a little later after stimulus onset, and there is also neural delay for many electrodes,
-        so one `pre_post` option could be [-1, 0, 0.2, 1.2] to add a buffer of 200ms for the
+        so one good `pre_post` option could be [-1, 0, 0.2, 1.2] to add a buffer of 200ms for the
         electrodes to begin responding more strongly.
     alpha : float, default=0.05
         Error rate.
+    average : bool, default=False
+        Whether to average each segment before t-test or not. Should only be used if have a sufficiently
+        high number of trials, because this will create one comparison per trial.
     fdr_method : str, {'indep', 'negcorr', None}, default='indep'
         Method of correction for multiple comparisons. If 'indep' it implements 
         Benjamini/Hochberg for independent or if 'negcorr' it corresponds
@@ -89,6 +92,19 @@ def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre
         - 'significant': True if null hypothesis was rejected (response is significantly different), False if not, shape (num_channels,)
         - 'alpha': error rate of the test
 
+    Examples
+    --------
+    >>> import naplib as nl
+    >>> data = nl.io.load_speech_task_data()
+    >>> # Get responsiveness of the 10 electrodes, assuming that they must show an increase
+    >>> # response to the stimulus, and averaging segment windows
+    >>> new_data, stats = responsive_ttest(data, average=True, alternative='less')
+    >>> # All 10 electrodes have significant responses
+    >>> stats['significant']
+    array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
+        True])
+
+
     References
     ----------
     .. [1] Mesgarani, N., & Chang, E. F. (2012). Selective cortical representation
@@ -144,8 +160,12 @@ def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre
             raise ValueError(f'befaft period is too short, there must be at least 5 samples of response before stimulus onset.')
         before_tmp = resp[t][pre_start:pre_end] - resp[t].mean(0, keepdims=True)
         after_tmp = resp[t][post_start:post_end] - resp[t].mean(0, keepdims=True)
-        before_samples.append(before_tmp)
-        after_samples.append(after_tmp)
+        if average:
+            before_samples.append(before_tmp.mean(0, keepdims=True))
+            after_samples.append(after_tmp.mean(0, keepdims=True))
+        else:
+            before_samples.append(before_tmp)
+            after_samples.append(after_tmp)
 
     before_samples_cat = np.concatenate(before_samples, axis=0)
     after_samples_cat = np.concatenate(after_samples, axis=0)
@@ -163,4 +183,4 @@ def responsive_ttest(data=None, resp='resp', befaft='befaft', sfreq='dataf', pre
 
     stats = {'pval': pval_corrected, 'stat': statistics, 'significant': reject, 'alpha': alpha}
 
-    return resp_corrected, stats
\ No newline at end of file
+    return resp_corrected, stats
diff --git a/tests/stats/test_responsive_elecs.py b/tests/stats/test_responsive_elecs.py
index 25d8fd3..1609254 100644
--- a/tests/stats/test_responsive_elecs.py
+++ b/tests/stats/test_responsive_elecs.py
@@ -25,6 +25,13 @@ def test_responsive_ttest_picks_correct_electrode_from_outstruct(outstruct):
     assert np.allclose(stats['stat'], np.array([-1.67113052, -63.20943622,   0.23296204,  -1.78460185]), atol=1e-7)
     assert np.allclose(stats['pval'], np.array([1.26958793e-001, 2.51319873e-266, 8.15870594e-001, 1.26958793e-001]), atol=1e-7)
 
+def test_responsive_ttest_picks_correct_electrode_from_outstruct_average_segment_windows(outstruct):
+    new_out, stats = responsive_ttest(data=outstruct, resp='resp', sfreq='dataf', befaft='befaft', average=True)
+    assert new_out[0].shape[1] == 1
+    assert np.array_equal(stats['significant'], np.array([0,1,0,0]).astype('bool'))
+    assert np.allclose(stats['stat'], np.array([-1.59007915, -5.28643157,  0.28904705, -1.57361706]), atol=1e-7)
+    assert np.allclose(stats['pval'], np.array([0.25424966, 0.02457324, 0.78690711, 0.25424966]), atol=1e-7)
+
 def test_responsive_ttest_picks_correct_electrode_from_outstruct_alternative_less(outstruct):
     new_out, stats = responsive_ttest(data=outstruct, resp='resp', sfreq='dataf', befaft='befaft', alternative='less')
     assert new_out[0].shape[1] == 1