[wip] replace Brain_Data.regress with nilearn and use new ResultsContainer class

Author: ejolly

nltools/data/__init__.py

@@ -5,11 +5,13 @@
 from .brain_data import Brain_Data, Groupby
 from .adjacency import Adjacency
 from .design_matrix import Design_Matrix, Design_Matrix_Series
+from .results import ResultsContainer
 
 __all__ = [
     "Brain_Data",
     "Adjacency",
     "Groupby",
     "Design_Matrix",
     "Design_Matrix_Series",
+    "ResultsContainer",
 ]

nltools/data/brain_data.py

@@ -10,6 +10,7 @@
 # Need to figure out how to speed up loading and resampling of data
 
 from nilearn.signal import clean
+from nilearn.glm.first_level import FirstLevelModel
 from scipy.stats import ttest_1samp, pearsonr, spearmanr
 from scipy.spatial.distance import cdist
 from scipy.stats import t as t_dist
@@ -63,7 +64,6 @@
     find_spikes,
     regress_permutation,
 )
-from nltools.stats import regress as regression
 from .adjacency import Adjacency
 from nltools.prefs import MNI_Template
 from nilearn.decoding import SearchLight
@@ -72,6 +72,9 @@
 from contextlib import redirect_stdout
 
 
+warnings.filterwarnings("ignore", category=UserWarning, module="nilearn")
+warnings.filterwarnings("ignore", category=RuntimeWarning, module="nilearn")
+
 # Optional dependencies
 nx = attempt_to_import("networkx", "nx")
 tables = attempt_to_import("tables")
@@ -879,26 +882,26 @@ def iplot(self, threshold=0, surface=False, anatomical=None, **kwargs):
             self, threshold=threshold, surface=surface, anatomical=anatomical, **kwargs
         )
 
-    def regress(self, mode="ols", **kwargs):
-        """Run a mass-univariate regression across voxels. Three types of regressions can be run:
-        1) Standard OLS (default)
-        2) Robust OLS (heteroscedasticty and/or auto-correlation robust errors), i.e. OLS with "sandwich estimators"
-        3) ARMA (auto-regressive and moving-average lags = 1 by default; experimental)
+    def regress(self, noise_model="ols", **kwargs):
+        """Runs a mass-univariate GLM analyses using the `Design_Matrix` supplied to `.X`
 
-        For more information see the help for nltools.stats.regress
+        This is a wrapper around [`nilearn.glm.first_level.FirstLevelModel`](https://nilearn.github.io/stable/modules/generated/nilearn.glm.first_level.FirstLevelModel.html#nilearn.glm.first_level.FirstLevelModel) which you can reference for additional information about what `**kwargs` are supported.
 
-        ARMA notes: This experimental mode is similar to AFNI's 3dREMLFit but without spatial smoothing of voxel auto-correlation estimates. It can be **very computationally intensive** so parallelization is used by default to try to speed things up. Speed is limited because a unique ARMA model is fit to *each voxel* (like AFNI/FSL), but unlike SPM, which assumes the same AR parameters (~0.2) at each voxel. While coefficient results are typically very similar to OLS, std-errors and so t-stats, dfs and and p-vals can differ greatly depending on how much auto-correlation is explaining the response in a voxel
-        relative to other regressors in the design matrix.
+        However, we override some defaults:
+        - no smoothing (use `.smooth()`)
+        - no scaling (use `.scale()`
+        - no drift model (should already be in the `Design_Matrix` set to `.X`)
+        - OLS noise model (use `noise_model = 'ar1'` to switch but takes more time)
 
         Args:
-            mode (str): kind of model to fit; must be one of 'ols' (default), 'robust', or 'arma'
-            kwargs (dict): keyword arguments to nltools.stats.regress
+            noise_model (str, optional): temporal variance model. Defaults to "ols"
 
-        Returns:
-            out: dictionary of regression statistics in Brain_Data instances
-                {'beta','t','p','df','residual'}
 
+        Returns:
+            ResultsContainer: with keys for each convolved column of `.X` and values as `Brain_Data` objects of the GLM statistics
         """
+        # Avoid circular import
+        from .results import ResultsContainer
 
         if not isinstance(self.X, pd.DataFrame):
             raise ValueError("Make sure self.X is a pandas DataFrame.")
@@ -909,43 +912,34 @@ def regress(self, mode="ols", **kwargs):
         if self.data.shape[0] != self.X.shape[0]:
             raise ValueError("self.X does not match the correct size of self.data")
 
-        b, se, t, p, df, res = regression(self.X, self.data, mode=mode, **kwargs)
-
-        # Prevent copy of all data in self multiple times; instead start with an empty instance and copy only needed attributes from self, and use this as a template for other outputs
-        b_out = self.__class__()
-        b_out.mask = deepcopy(self.mask)
-        b_out.nifti_masker = deepcopy(self.nifti_masker)
-
-        # Use this as template for other outputs before setting data
-        se_out = b_out.copy()
-        t_out = b_out.copy()
-        p_out = b_out.copy()
-        df_out = b_out.copy()
-        res_out = b_out.copy()
-        (
-            b_out.data,
-            se_out.data,
-            t_out.data,
-            p_out.data,
-            df_out.data,
-            res_out.data,
-        ) = (
-            b,
-            se,
-            t,
-            p,
-            df,
-            res,
+        # Nilearn FirstLevelModel default overrides
+        smoothing_fwhm = kwargs.get("smoothing_fwhm", None)
+        drift_model = kwargs.get("drift_model", None)
+        signal_scaling = kwargs.get("signal_scaling", False)
+        stat_type = kwargs.get("stat_type", "t")
+        output_type = kwargs.get("output_type", "all")
+
+        # Run GLM
+        glm = FirstLevelModel(
+            t_r=1 / self.X.sampling_freq,
+            mask_img=self.mask,
+            smoothing_fwhm=smoothing_fwhm,
+            drift_model=drift_model,
+            signal_scaling=signal_scaling,
+            noise_model=noise_model,
         )
-
-        return {
-            "beta": b_out,
-            "t": t_out,
-            "p": p_out,
-            "df": df_out,
-            "sigma": se_out,
-            "residual": res_out,
+        glm.fit(self.to_nifti(), design_matrices=self.X)
+        self.glm = glm
+
+        # Assemble results
+        regressors_of_interest = self.X.convolved
+        results = {
+            r: ResultsContainer(
+                glm.compute_contrast(r, stat_type=stat_type, output_type=output_type)
+            )
+            for r in regressors_of_interest
         }
+        return results
 
     def randomise(
         self, n_permute=5000, threshold_dict=None, return_mask=False, **kwargs

nltools/data/results.py

@@ -0,0 +1,63 @@
+from typing import Dict, Any
+import nibabel as nib
+from nltools.data import Brain_Data
+
+
+class ResultsContainer(object):
+    """A genericcontainer that dynamically creates attributes from a dictionary of string: Nifti1Image entries, where each attribute is a Brain_Data instance initialized from the corresponding Nifti image.
+
+    Args:
+        images_dict (Dict[str, nib.Nifti1Image]): Dictionary mapping attribute names to Nifti images.
+
+    Example:
+        >>> rc = ResultsContainer({'foo': img1, 'bar': img2})
+        >>> rc.foo  # Brain_Data instance
+        >>> rc.bar  # Brain_Data instance
+    """
+
+    def __init__(self, images=None):
+        self._is_single = True
+        self.data = []
+        if isinstance(images, dict):
+            for key, img in images.items():
+                if not isinstance(img, nib.Nifti1Image):
+                    raise TypeError(
+                        f"Value for key '{key}' is not a Nifti1Image it's a {type(img)}."
+                    )
+                if key == "stat":
+                    new_key = "t"
+                elif key == "p_value":
+                    new_key = "p"
+                elif key == "effect_size":
+                    new_key = "beta"
+                elif key == "effect_variance":
+                    new_key = "se"
+                else:
+                    new_key = key
+                setattr(self, new_key, Brain_Data(img))
+        elif isinstance(images, list):
+            self._is_single = False
+            for img in images:
+                self.data.append(ResultsContainer(img))
+
+    def __repr__(self):
+        attr_names = [k for k in self.__dict__.keys()]
+        return f"ResultsContainer(attributes={attr_names})"
+
+    def append(self, image_dict):
+        self._is_single = False
+        # image_dict = (
+        #     ResultsContainer(image_dict)
+        #     if not isinstance(image_dict, ResultsContainer)
+        #     else image_dict
+        # )
+        self.data.append(image_dict)
+
+    def __getitem__(self, index):
+        return self.data[index]
+
+    def __len__(self):
+        return len(self.data)
+
+    def __iter__(self):
+        return iter(self.data)

nltools/file_reader.py

@@ -86,7 +86,9 @@ def onsets_to_dm(
             polys = [c for c in dm.columns if "drift" in c or "constant" in c]
         else:
             convolved, polys = [], []
-        dm = Design_Matrix(dm, convolved=convolved, sampling_freq=1 / TR, polys=polys)
+        dm = Design_Matrix(
+            dm, convolved=convolved, sampling_freq=1 / TR, polys=polys
+        ).reset_index(drop=True)
         out.append(dm)
 
     return out if len(out) > 1 else out[0]