diff --git a/docs/changes/newsfragments/287.enh b/docs/changes/newsfragments/287.enh
new file mode 100644
index 000000000..e389e2a99
--- /dev/null
+++ b/docs/changes/newsfragments/287.enh
@@ -0,0 +1 @@
+Add ``mode`` as an aggregation function option in :func:`.get_aggfunc_by_name` by `Synchon Mandal`_
diff --git a/junifer/api/tests/test_api_utils.py b/junifer/api/tests/test_api_utils.py
index f466cdda3..de5aad532 100644
--- a/junifer/api/tests/test_api_utils.py
+++ b/junifer/api/tests/test_api_utils.py
@@ -35,6 +35,7 @@ def test_get_dependency_information_short() -> None:
     assert list(dependency_information.keys()) == [
         "click",
         "numpy",
+        "scipy",
         "datalad",
         "pandas",
         "nibabel",
@@ -51,6 +52,7 @@ def test_get_dependency_information_long() -> None:
     for key in [
         "click",
         "numpy",
+        "scipy",
         "datalad",
         "pandas",
         "nibabel",
diff --git a/junifer/stats.py b/junifer/stats.py
index bf4b133e2..aa411f832 100644
--- a/junifer/stats.py
+++ b/junifer/stats.py
@@ -7,7 +7,7 @@
 from typing import Any, Callable, Dict, List, Optional
 
 import numpy as np
-from scipy.stats import trim_mean
+from scipy.stats import mode, trim_mean
 from scipy.stats.mstats import winsorize
 
 from .utils import logger, raise_error
@@ -24,10 +24,11 @@ def get_aggfunc_by_name(
         Name to identify the function. Currently supported names and
         corresponding functions are:
 
-        * ``winsorized_mean`` -> :func:`scipy.stats.mstats.winsorize`
         * ``mean`` -> :func:`numpy.mean`
-        * ``std`` -> :func:`numpy.std`
+        * ``winsorized_mean`` -> :func:`scipy.stats.mstats.winsorize`
         * ``trim_mean`` -> :func:`scipy.stats.trim_mean`
+        * ``mode`` -> :func:`scipy.stats.mode`
+        * ``std`` -> :func:`numpy.std`
         * ``count`` -> :func:`.count`
         * ``select`` -> :func:`.select`
 
@@ -40,6 +41,7 @@ def get_aggfunc_by_name(
     -------
     function
         Respective function with ``func_params`` parameter set.
+
     """
     from functools import partial  # local import to avoid sphinx error
 
@@ -51,6 +53,7 @@ def get_aggfunc_by_name(
         "trim_mean",
         "count",
         "select",
+        "mode",
     }
     if func_params is None:
         func_params = {}
@@ -93,6 +96,8 @@ def get_aggfunc_by_name(
         elif pick is not None and drop is not None:
             raise_error("Either pick or drop must be specified, not both.")
         func = partial(select, **func_params)
+    elif name == "mode":
+        func = partial(mode, **func_params)
     else:
         raise_error(
             f"Function {name} unknown. Please provide any of "
@@ -115,6 +120,7 @@ def count(data: np.ndarray, axis: int = 0) -> np.ndarray:
     -------
     numpy.ndarray
         Number of elements along the given axis.
+
     """
     ax_size = data.shape[axis]
     if axis < 0:
@@ -137,7 +143,7 @@ def winsorized_mean(
         The axis to calculate winsorized mean on (default None).
     **win_params : dict
         Dictionary containing the keyword arguments for the winsorize function.
-        E.g. ``{'limits': [0.1, 0.1]}``.
+        E.g., ``{'limits': [0.1, 0.1]}``.
 
     Returns
     -------
@@ -149,6 +155,7 @@ def winsorized_mean(
     --------
     scipy.stats.mstats.winsorize :
         The winsorize function used in this function.
+
     """
     win_dat = winsorize(data, axis=axis, **win_params)
     win_mean = win_dat.mean(axis=axis)
@@ -180,6 +187,13 @@ def select(
     numpy.ndarray
         Subset of the inputted data with the select settings
         applied as specified in ``select_params``.
+
+    Raises
+    ------
+    ValueError
+        If both ``pick`` and ``drop`` are None or
+        if both ``pick`` and ``drop`` are not None.
+
     """
 
     if pick is None and drop is None:
diff --git a/junifer/tests/test_stats.py b/junifer/tests/test_stats.py
index d8aed8390..900c765e6 100644
--- a/junifer/tests/test_stats.py
+++ b/junifer/tests/test_stats.py
@@ -21,6 +21,8 @@ from junifer.stats import count, get_aggfunc_by_name, select, winsorized_mean
         ("count", None),
         ("trim_mean", None),
         ("trim_mean", {"proportiontocut": 0.1}),
+        ("mode", None),
+        ("mode", {"keepdims": True}),
     ],
 )
 def test_get_aggfunc_by_name(name: str, params: Optional[Dict]) -> None:
diff --git a/pyproject.toml b/pyproject.toml
index 13e09a6df..213bc37f7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -37,6 +37,7 @@ classifiers = [
 dependencies = [
     "click>=8.1.3,<8.2",
     "numpy>=1.24,<1.27",
+    "scipy>=1.9.0,<=1.11.4",
     "datalad>=0.15.4,<0.20",
     "pandas>=1.4.0,<2.2",
     "nibabel>=3.2.0,<5.11",
@@ -177,6 +178,7 @@ known-first-party = ["junifer"]
 known-third-party =[
     "click",
     "numpy",
+    "scipy",
     "datalad",
     "pandas",
     "nibabel",