[ENH] pad t and F 'short' contrasts with zeros

doc/changes/latest.rst

@@ -20,6 +20,7 @@ Enhancements
 ------------
 
 - Allow setting ``vmin`` in :func:`~nilearn.plotting.plot_glass_brain` and :func:`~nilearn.plotting.plot_stat_map` (:gh:`3993` by `Michelle Wang`_).
+- :bdg-success:`API` Support passing t and F contrasts to :func:`~nilearn.glm.contrasts.compute_contrast` that that have fewer columns than the number of estimated parameters. Remaining columns are padded with zero (:gh:`4067` by `Remi Gau`_).
 
 Changes
 -------

nilearn/glm/_utils.py

@@ -343,3 +343,63 @@ def positive_reciprocal(X):
     """
     X = np.asarray(X)
     return np.where(X <= 0, 0, 1.0 / X)
+
+
+def pad_contrast(con_val, theta, contrast_type):
+    """Pad contrast with zeros if necessary.
+
+    If the contrast is shorter than the number of parameters,
+    it is padded with zeros.
+
+    If the contrast is longer than the number of parameters,
+    a ValueError is raised.
+
+    Parameters
+    ----------
+    con_val : numpy.ndarray of shape (p) or (n, p)
+        Where p = number of regressors
+        with a value explicitly passed by the user.
+        p must be <= P,
+        where P is the total number of regressors in the design matrix.
+
+    theta : numpy.ndarray with shape (P,m)
+        theta of RegressionResults instances
+        where P is the total number of regressors in the design matrix.
+
+    contrast_type : {'t', 'F'}, optional
+        Type of the :term:`contrast`.
+    """
+    nb_cols = con_val.shape[0] if con_val.ndim == 1 else con_val.shape[1]
+    if nb_cols > theta.shape[0]:
+        if contrast_type == "t":
+            raise ValueError(
+                f"t contrasts should be of length P={theta.shape[0]}, "
+                f"but it has length {nb_cols}."
+            )
+        if contrast_type == "F":
+            raise ValueError(
+                f"F contrasts should have {theta.shape[0]} columns, "
+                f"but it has {nb_cols}."
+            )
+
+    pad = False
+    if nb_cols < theta.shape[0]:
+        pad = True
+        if contrast_type == "t":
+            warn(
+                f"t contrasts should be of length P={theta.shape[0]}, "
+                f"but it has length {nb_cols}. "
+                "The rest of the contrast was padded with zeros."
+            )
+        if contrast_type == "F":
+            warn(
+                f"F contrasts should have {theta.shape[0]} colmuns, "
+                f"but it has only {nb_cols}. "
+                "The rest of the contrast was padded with zeros."
+            )
+
+    if pad:
+        padding = np.zeros((nb_cols, theta.shape[0] - nb_cols))
+        con_val = np.hstack((con_val, padding))
+
+    return con_val

nilearn/glm/contrasts.py

@@ -10,7 +10,7 @@
 import pandas as pd
 import scipy.stats as sps
 
-from nilearn.glm._utils import z_score
+from nilearn.glm._utils import pad_contrast, z_score
 from nilearn.maskers import NiftiMasker
 
 DEF_TINY = 1e-50
@@ -95,17 +95,25 @@ def compute_contrast(labels, regression_result, con_val, contrast_type=None):
         var_ = np.zeros(labels.size)
         for label_ in regression_result:
             label_mask = labels == label_
-            resl = regression_result[label_].Tcontrast(con_val)
-            effect_[:, label_mask] = resl.effect.T
-            var_[label_mask] = (resl.sd**2).T
+            reg = regression_result[label_].Tcontrast(con_val)
+            effect_[:, label_mask] = reg.effect.T
+            var_[label_mask] = (reg.sd**2).T
+
     elif contrast_type == "F":
         from scipy.linalg import sqrtm
 
         effect_ = np.zeros((dim, labels.size))
         var_ = np.zeros(labels.size)
+        # TODO
+        # explain why we cannot simply do
+        # reg = regression_result[label_].Tcontrast(con_val)
+        # like above or refactor the code so it can be done
         for label_ in regression_result:
             label_mask = labels == label_
             reg = regression_result[label_]
+            con_val = pad_contrast(
+                con_val=con_val, theta=reg.theta, contrast_type=contrast_type
+            )
             cbeta = np.atleast_2d(np.dot(con_val, reg.theta))
             invcov = np.linalg.inv(
                 np.atleast_2d(reg.vcov(matrix=con_val, dispersion=1.0))

nilearn/glm/model.py

@@ -7,7 +7,7 @@
 from scipy.linalg import inv
 from scipy.stats import t as t_distribution
 
-from nilearn.glm._utils import positive_reciprocal
+from nilearn.glm._utils import pad_contrast, positive_reciprocal
 
 # Inverse t cumulative distribution
 inv_t_cdf = t_distribution.ppf
@@ -196,11 +196,9 @@ def Tcontrast(self, matrix, store=("t", "effect", "sd"), dispersion=None):
             matrix = matrix[None]
         if matrix.shape[0] != 1:
             raise ValueError("t contrasts should have only one row")
-        if matrix.shape[1] != self.theta.shape[0]:
-            raise ValueError(
-                f"t contrasts should be length P={self.theta.shape[0]}, "
-                f"but this is length {matrix.shape[1]}"
-            )
+        matrix = pad_contrast(
+            con_val=matrix, theta=self.theta, contrast_type="t"
+        )
         store = set(store)
         if not store.issubset(("t", "effect", "sd")):
             raise ValueError(f"Unexpected store request in {store}")
@@ -268,6 +266,9 @@ def Fcontrast(self, matrix, dispersion=None, invcov=None):
                 f"F contrasts should have shape[1] P={self.theta.shape[0]}, "
                 f"but this has shape[1] {matrix.shape[1]}"
             )
+        matrix = pad_contrast(
+            con_val=matrix, theta=self.theta, contrast_type="F"
+        )
         ctheta = np.dot(matrix, self.theta)
         if matrix.ndim == 1:
             matrix = matrix.reshape((1, matrix.shape[0]))

nilearn/glm/tests/test_contrasts.py

@@ -229,3 +229,19 @@ def test_invalid_contarst_type():
     variance = np.ones(1)
     with pytest.raises(ValueError, match="is not a valid contrast_type."):
         Contrast(effect, variance, contrast_type="foo")
+
+
+def test_contrast_padding(rng):
+    n, p, q = 100, 80, 10
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
+    labels, results = run_glm(Y, X, "ar1")
+
+    con_val = [1]
+
+    with pytest.warns(
+        UserWarning, match="The rest of the contrast was padded with zeros."
+    ):
+        compute_contrast(labels, results, con_val).z_score()
+
+    con_val = np.eye(q)[:3, :3]
+    compute_contrast(labels, results, con_val, contrast_type="F").z_score()

nilearn/glm/tests/test_model.py

@@ -66,8 +66,6 @@ def test_t_contrast():
     assert_array_almost_equal(RESULTS.Tcontrast([1, 0]).t, 3.25)
     assert_array_almost_equal(RESULTS.Tcontrast([0, 1]).t, 7.181, 3)
     # Input matrix checked for size
-    with pytest.raises(ValueError):
-        RESULTS.Tcontrast([1])
     with pytest.raises(ValueError):
         RESULTS.Tcontrast([1, 0, 0])
     # And shape

nilearn/glm/tests/test_second_level.py

@@ -1059,11 +1059,6 @@ def test_second_level_contrast_computation_errors(tmp_path, rng):
         model.compute_contrast(cnull)
 
     # passing wrong parameters
-    with pytest.raises(
-        ValueError,
-        match=("t contrasts should be length P=1, but this is length 0"),
-    ):
-        model.compute_contrast(second_level_contrast=[])
     with pytest.raises(ValueError, match="Allowed types are .*'t', 'F'"):
         model.compute_contrast(
             second_level_contrast=c1, second_level_stat_type=""
@@ -1083,11 +1078,40 @@ def test_second_level_contrast_computation_errors(tmp_path, rng):
         ValueError, match="No second-level contrast is specified"
     ):
         model.compute_contrast(None)
+
+
+def test_second_level_t_contrast_length_errors(tmp_path):
+    func_img, mask = fake_fmri_data(file_path=tmp_path)
+
+    model = SecondLevelModel(mask_img=mask)
+
+    func_img, mask = fake_fmri_data(file_path=tmp_path)
+    Y = [func_img] * 4
+    X = pd.DataFrame([[1]] * 4, columns=["intercept"])
+    model = model.fit(Y, design_matrix=X)
+
+    with pytest.raises(
+        ValueError,
+        match=("t contrasts should be of length P=1, but it has length 2."),
+    ):
+        model.compute_contrast(second_level_contrast=[1, 2])
+
+
+def test_second_level_F_contrast_length_errors(tmp_path):
+    func_img, mask = fake_fmri_data(file_path=tmp_path)
+
+    model = SecondLevelModel(mask_img=mask)
+
+    func_img, mask = fake_fmri_data(file_path=tmp_path)
+    Y = [func_img] * 4
+    X = pd.DataFrame([[1]] * 4, columns=["intercept"])
+    model = model.fit(Y, design_matrix=X)
+
     with pytest.raises(
         ValueError,
-        match=("t contrasts should be length P=2, but this is length 1"),
+        match=("F contrasts should have .* columns, but it has .*"),
     ):
-        model.compute_contrast([1])
+        model.compute_contrast(second_level_contrast=np.eye(2))
 
 
 def test_non_parametric_inference_contrast_computation(tmp_path):