Version 0.3.3 (#262)

* [FIX] Remove pyembree requirement * [FIX] Remove pyembree requirement (#255) * [FIX] Fix data inclusion * [DOC] Update tutorial 2. * styling * [DOC] Update docs for pyembree fix (#257) * [FIX] Remove pyembree requirement * [FIX] Fix data inclusion * [DOC] Update tutorial 2. * styling * [FIX] Fix data inclusion * [FIX] Fix data inclusion * [FIX] Remove pyembree requirement * [FIX] Fix data inclusion * [DOC] Update tutorial 2. * styling * [FIX] Fix data inclusion * [FIX] Fix readthedocs generation * [ENH] Add CIVET to histology module (#259) * [ENH] Add CIVET to histology module * styling * styling * readthedocs fix * [ENH] Enable hotelling T2 for more than 3 variates (#261) * [FIX] Fix automatic naming of mixed effects (#266) * [FIX] Fix error in automatic naming for mixed effects * fix mypy * [FIX] Remove a mustbefolder (#267) * [ENH] Optimisation of finding peaks and clusters (#268) * Optimisation of finding peaks and clusters Significantly faster way of finding cluster ids - especially for large number of vertex above thresh * [FIX] FIX clus_id and black lint Co-authored-by: Reinder Vos de Wael <[email protected]> * [ENH, DOC] Added MICs, updated tutorial * update MATLAB tutorials & functions * add rater check to abide data * styling * docsting updates * add MICS dataset * add verbose flag to _download_file * styling * add mask loading for all templates * matlab tutorial 1 new * flatten slm.P.pval.C * tutorial update * slm now allows pandas input, updatepython tutorial 2 * fix slm.contrast * save for checkout * tutorial updates * version bump Co-authored-by: Lasse Madsen <[email protected]>
MICA-MNI · Dec 17, 2021 · dff5577 · dff5577
1 parent ac3c356
commit dff5577
Show file tree

Hide file tree

Showing 76 changed files with 2,889 additions and 1,409 deletions.
diff --git a/brainstat/__init__.py b/brainstat/__init__.py
@@ -1,2 +1,2 @@
 """Neuroimaging statistics toolbox."""
-__version__ = "0.3.2"
+__version__ = "0.3.3"
diff --git a/brainstat/_utils.py b/brainstat/_utils.py
@@ -22,6 +22,7 @@
     "ABIDE_DATA_DIR": BRAINSTAT_DATA_DIR / "abide_data",
     "BIGBRAIN_DATA_DIR": BRAINSTAT_DATA_DIR / "bigbrain_data",
     "GRADIENT_DATA_DIR": BRAINSTAT_DATA_DIR / "gradient_data",
+    "MICS_DATA_DIR": BRAINSTAT_DATA_DIR / "mics_data",
     "NEUROSYNTH_DATA_DIR": BRAINSTAT_DATA_DIR / "neurosynth_data",
     "PARCELLATION_DATA_DIR": BRAINSTAT_DATA_DIR / "parcellation_data",
     "SURFACE_DATA_DIR": BRAINSTAT_DATA_DIR / "surface_data",
@@ -81,6 +82,15 @@ def generate_data_fetcher_json() -> None:
             "civet164k": {
                 "url": "https://box.bic.mni.mcgill.ca/s/rei5HtTDvexlEPA/download"
             },
+            "fsaverage5": {
+                "url": "https://box.bic.mni.mcgill.ca/s/jsSDyiDcm9KEQpf/download"
+            },
+            "fsaverage": {
+                "url": "https://box.bic.mni.mcgill.ca/s/XiZx9HfHaufb4kD/download"
+            },
+            "fslr32k": {
+                "url": "https://box.bic.mni.mcgill.ca/s/cFXCrSkfiJFjUJ0/download"
+            },
         },
         "spheres": {
             "civet41k": {
@@ -92,6 +102,14 @@ def generate_data_fetcher_json() -> None:
                 "url": "https://s3.amazonaws.com/fcp-indi/data/Projects/ABIDE_Initiative/Phenotypic_V1_0b_preprocessed1.csv"
             },
         },
+        "mics_tutorial": {
+            "demographics": {
+                "url": "https://box.bic.mni.mcgill.ca/s/7bW9JIpvQvSJeuf/download"
+            },
+            "thickness": {
+                "url": "https://box.bic.mni.mcgill.ca/s/kMi6lU91piwdaCf/download"
+            },
+        },
     }
     with open(json_file, "w") as f:
         json.dump(data, f, indent=4, sort_keys=True)
@@ -129,7 +147,9 @@ def deprecated_func(*args, **kwargs):
     return deprecated_decorator
 
 
-def _download_file(url: str, output_file: Path, overwrite: bool = False) -> None:
+def _download_file(
+    url: str, output_file: Path, overwrite: bool = False, verbose=True
+) -> None:
     """Downloads a file.
 
     Parameters
@@ -140,12 +160,15 @@ def _download_file(url: str, output_file: Path, overwrite: bool = False) -> None
         Path object of the output file.
     overwrite : bool
         If true, overwrite existing files, defaults to False.
+    verbose : bool
+        If true, print a download message, defaults to True.
     """
 
     if output_file.exists() and not overwrite:
         return
 
-    logger.info("Downloading " + str(output_file) + " from " + url + ".")
+    if verbose:
+        logger.info("Downloading " + str(output_file) + " from " + url + ".")
     with urllib.request.urlopen(url) as response, open(output_file, "wb") as out_file:
         shutil.copyfileobj(response, out_file)
 

diff --git a/brainstat/context/histology.py b/brainstat/context/histology.py
@@ -137,7 +137,7 @@ def read_histology_profile(
     data_dir = Path(data_dir) if data_dir else data_directories["BIGBRAIN_DATA_DIR"]
 
     if template[:5] == "civet":
-        logger.warn(
+        logger.info(
             "CIVET histology profiles were not included with BigBrainWarp. Interpolating from fsaverage using nearest neighbor interpolation."
         )
         civet_template = template

diff --git a/brainstat/data/data_urls.json b/brainstat/data/data_urls.json
@@ -26,6 +26,23 @@
         },
         "civet41k": {
             "url": "https://box.bic.mni.mcgill.ca/s/9kzBetBCZkkqN6w/download"
+        },
+        "fsaverage": {
+            "url": "https://box.bic.mni.mcgill.ca/s/XiZx9HfHaufb4kD/download"
+        },
+        "fsaverage5": {
+            "url": "https://box.bic.mni.mcgill.ca/s/jsSDyiDcm9KEQpf/download"
+        },
+        "fslr32k": {
+            "url": "https://box.bic.mni.mcgill.ca/s/cFXCrSkfiJFjUJ0/download"
+        }
+    },
+    "mics_tutorial": {
+        "demographics": {
+            "url": "https://box.bic.mni.mcgill.ca/s/7bW9JIpvQvSJeuf/download"
+        },
+        "thickness": {
+            "url": "https://box.bic.mni.mcgill.ca/s/kMi6lU91piwdaCf/download"
         }
     },
     "neurosynth_precomputed": {

diff --git a/brainstat/datasets/base.py b/brainstat/datasets/base.py
@@ -65,7 +65,7 @@ def fetch_parcellation(
     data_dir.mkdir(parents=True, exist_ok=True)
 
     if template == "civet41k" or template == "civet164k":
-        logger.warn(
+        logger.info(
             "CIVET parcellations were not included with the toolbox. Interpolating parcellation from the fsaverage surface with a nearest neighbor interpolation."
         )
         civet_template = template
@@ -169,7 +169,8 @@ def fetch_mask(
     Parameters
     ----------
     template : str
-        Name of the surface template. Valid templates are "civet41k" and "civet164k".
+        Name of the surface template. Valid templates are: "fsaverage5",
+        "fsaverage", "fslr32k", "civet41k", and "civet164k".
     join : bool, optional
         If true, returns a numpy array containing the mask. If false, returns a
         tuple containing the left and right hemispheric masks, respectively, by
@@ -183,6 +184,7 @@ def fetch_mask(
         Midline mask, either as a single array or a tuple of a left and right hemispheric
         array.
     """
+
     data_dir = Path(data_dir) if data_dir else data_directories["SURFACE_DATA_DIR"]
     data_dir.mkdir(parents=True, exist_ok=True)
 
@@ -235,7 +237,7 @@ def fetch_gradients(
 
     hf = h5py.File(gradients_file, "r")
     if template == "civet41k" or template == "civet164k":
-        logger.warn(
+        logger.info(
             "CIVET gradients were not included with the toolbox. Interpolating gradients from the fsaverage surface with a nearest interpolation."
         )
         fsaverage_left, fsaverage_right = fetch_template_surface(

diff --git a/brainstat/stats/SLM.py b/brainstat/stats/SLM.py
@@ -31,7 +31,7 @@ class SLM:
     def __init__(
         self,
         model: Union[FixedEffect, MixedEffect],
-        contrast: Union[ArrayLike, FixedEffect],
+        contrast: ArrayLike,
         surf: Optional[Union[str, dict, BSPolyData]] = None,
         mask: Optional[ArrayLike] = None,
         *,
@@ -49,7 +49,7 @@ def __init__(
         ----------
         model : brainstat.stats.terms.FixedEffect, brainstat.stats.terms.MixedEffect
             The linear model to be fitted of dimensions (observations, predictors).
-        contrast : array-like, brainstat.stats.terms.FixedEffect
+        contrast : array-like
             Vector of contrasts in the observations.
         surf : str, dict, BSPolyData, optional
             A surface provided as either a dictionary with keys 'tri' for its
@@ -85,7 +85,7 @@ def __init__(
         """
         # Input arguments.
         self.model = model
-        self.contrast = contrast
+        self.contrast = np.array(contrast)
 
         if isinstance(surf, str):
             self.surf_name = surf
@@ -133,6 +133,10 @@ def fit(self, Y: np.ndarray) -> None:
                 raise NotImplementedError(
                     "One-tailed tests are not implemented for multivariate data."
                 )
+            if Y.shape[2] > 3 and "rft" in self.correction:
+                raise NotImplementedError(
+                    "Random Field Theory corrections are not implemented for more than three variates."
+                )
             student_t_test = Y.shape[2] == 1
         else:
             student_t_test = True

diff --git a/brainstat/stats/_multiple_comparisons.py b/brainstat/stats/_multiple_comparisons.py
@@ -229,6 +229,8 @@ def _random_field_theory(self) -> valid_rft_output:
     )[0]
     tlim = tlim[1]
     pval["P"] = pval["P"] * (self.t[0, :] > tlim) + (self.t[0, :] <= tlim)
+    if pval["C"] is not None:
+        pval["C"] = pval["C"].flatten()
 
     return pval, peak, clus, clusid
 
@@ -916,20 +918,37 @@ def peak_clus(
     edg = voxid[edg[np.all(excurset[edg], 1), :]]
     nf = np.arange(1, n + 1)
 
-    # Find cluster id's in nf (from Numerical Recipes in C, page 346):
-    for el in range(1, edg.shape[0] + 1):
-        j = edg[el - 1, 0]
-        k = edg[el - 1, 1]
-        while nf[j - 1] != j:
-            j = nf[j - 1]
-        while nf[k - 1] != k:
-            k = nf[k - 1]
-        if j != k:
-            nf[j - 1] = k
-
-    for j in range(1, n + 1):
-        while nf[j - 1] != nf[nf[j - 1] - 1]:
-            nf[j - 1] = nf[nf[j - 1] - 1]
+    # Find cluster id's in nf
+    not_used_indexes = set(np.arange(1, n + 1))
+
+    cluster_val = 1
+    for i in range(1, n + 1):
+        if i not in not_used_indexes:
+            continue
+        in_cluster = set([i])
+        not_used_indexes.remove(i)
+
+        neighbours = (
+            set(edg[np.isin(edg, list(in_cluster)).any(axis=1)].ravel())
+            & not_used_indexes
+        )
+        in_cluster = in_cluster | neighbours
+
+        while True:
+            neighbours = (
+                set(edg[np.isin(edg, list(neighbours)).any(axis=1)].ravel())
+                & not_used_indexes
+            )
+            not_used_indexes = not_used_indexes - neighbours
+            if len(neighbours) == 0:
+                break
+            else:
+                in_cluster = in_cluster | neighbours
+
+        in_cluster_idx = [x - 1 for x in list(in_cluster)]
+        nf[in_cluster_idx] = cluster_val
+
+        cluster_val = cluster_val + 1
 
     vox = np.argwhere(excurset) + 1
     ivox = np.argwhere(np.in1d(vox, lmvox)) + 1

diff --git a/brainstat/stats/_t_test.py b/brainstat/stats/_t_test.py
@@ -147,98 +147,32 @@ def _t_test(self) -> None:
         vf = np.divide(np.sum(np.square(np.dot(c.T, pinvX)), axis=1), self.df)
         self.sd = np.sqrt(vf * self.SSE[jj, :])
 
-        if k == 2:
-            det = np.multiply(self.SSE[0, :], self.SSE[2, :]) - np.square(
-                self.SSE[1, :]
-            )
-
-            self.t = (
-                np.multiply(np.square(self.ef[0, :]), self.SSE[2, :])
-                + np.multiply(np.square(self.ef[1, :]), self.SSE[0, :])
-                - np.multiply(
-                    np.multiply(2 * self.ef[0, :], self.ef[1, :]), self.SSE[1, :]
-                )
-            )
-
-        if k == 3:
-            det = (
-                np.multiply(
-                    self.SSE[0, :],
-                    (
-                        np.multiply(self.SSE[2, :], self.SSE[5, :])
-                        - np.square(self.SSE[4, :])
-                    ),
-                )
-                - np.multiply(self.SSE[5, :], np.square(self.SSE[1, :]))
-                + np.multiply(
-                    self.SSE[3, :],
-                    (
-                        np.multiply(self.SSE[1, :], self.SSE[4, :]) * 2
-                        - np.multiply(self.SSE[2, :], self.SSE[3, :])
-                    ),
-                )
-            )
-
-            self.t = np.multiply(
-                np.square(self.ef[0, :]),
-                (
-                    np.multiply(self.SSE[2, :], self.SSE[5, :])
-                    - np.square(self.SSE[4, :])
-                ),
-            )
-
-            self.t = self.t + np.multiply(
-                np.square(self.ef[1, :]),
-                (
-                    np.multiply(self.SSE[0, :], self.SSE[5, :])
-                    - np.square(self.SSE[3, :])
-                ),
-            )
-
-            self.t = self.t + np.multiply(
-                np.square(self.ef[2, :]),
-                (
-                    np.multiply(self.SSE[0, :], self.SSE[2, :])
-                    - np.square(self.SSE[1, :])
-                ),
-            )
-
-            self.t = self.t + np.multiply(
-                2 * self.ef[0, :],
-                np.multiply(
-                    self.ef[1, :],
-                    (
-                        np.multiply(self.SSE[3, :], self.SSE[4, :])
-                        - np.multiply(self.SSE[1, :], self.SSE[5, :])
-                    ),
-                ),
-            )
-
-            self.t = self.t + np.multiply(
-                2 * self.ef[0, :],
-                np.multiply(
-                    self.ef[2, :],
-                    (
-                        np.multiply(self.SSE[1, :], self.SSE[4, :])
-                        - np.multiply(self.SSE[2, :], self.SSE[3, :])
-                    ),
-                ),
-            )
-
-            self.t = self.t + np.multiply(
-                2 * self.ef[1, :],
-                np.multiply(
-                    self.ef[2, :],
-                    (
-                        np.multiply(self.SSE[1, :], self.SSE[3, :])
-                        - np.multiply(self.SSE[0, :], self.SSE[4, :])
-                    ),
-                ),
-            )
-
-        if k > 3:
-            sys.exit("Hotelling" "s T for k>3 not programmed yet")
-
-        self.t = np.multiply(np.divide(self.t, (det + (det <= 0))), (det > 0)) / vf
-        self.t = np.multiply(np.sqrt(self.t + (self.t <= 0)), (self.t > 0))
+        # Initialize some variables
+        self.t = np.zeros(v)
+        indices_lower = np.tril_indices(self.k)
+        sse_indices = np.zeros((self.k, self.k), dtype=int)
+        sse_indices[indices_lower] = np.arange(0, self.k * (self.k + 1) / 2)
+        sse_indices += np.tril(sse_indices, -1).T
+
+        M = np.zeros((self.k + 1, self.k + 1))
+        M_ef = np.zeros((self.k + 1, self.k + 1), dtype=bool)
+        M_ef[1:, 0] = True
+        M_ef[0, 1:] = True
+        M_sse = np.zeros((self.k + 1, self.k + 1), dtype=bool)
+        M_sse[1:, 1:] = True
+
+        ef_duplicate = np.concatenate((self.ef, self.ef), axis=0)
+        for i in range(v):
+            sse_vertex = self.SSE[:, i]
+            sse_matrix = sse_vertex[sse_indices]
+
+            det_sse = np.linalg.det(sse_matrix)
+            if det_sse <= 0:
+                self.t[i] = 0
+            else:
+                M[M_ef] = ef_duplicate[:, i]
+                M[M_sse] = sse_matrix.flatten()
+
+                self.t[i] = np.sqrt(-np.linalg.det(M) / det_sse / vf)
+
     self.t = np.atleast_2d(self.t)
diff --git a/brainstat/stats/terms.py b/brainstat/stats/terms.py
@@ -1,4 +1,5 @@
 """Classes for fixed, mixed, and random effects."""
+import difflib
 import re
 from typing import Any, List, Optional, Sequence, Union
 
@@ -484,12 +485,19 @@ def __init__(
                     if v != 1:
                         name_ran += f"{v}**2"
                 else:
-                    name = check_names(ran)
-                    if name is not None and name_ran is None:
-                        name_ran = name
+                    if name_ran is None:
+                        name = check_names(ran)
+                        if name:
+                            name_ran = name[0]
+                            for i in range(1, len(name)):
+                                sm = difflib.SequenceMatcher(None, name_ran, name[i])
+                                match = sm.find_longest_match(
+                                    0, len(name_ran), 0, len(name[i])
+                                )
+                                name_ran = name_ran[match.a : match.a + match.size]
+
                 ran = ran @ ran.T
                 ran = ran.values.ravel()
-
             self.variance = FixedEffect(ran, names=name_ran, add_intercept=False)
         self.mean = FixedEffect(fix, names=name_fix, add_intercept=add_intercept)