diff --git a/mne/stats/regression.py b/mne/stats/regression.py
index bb696f6ef05..d6343cfccba 100644
--- a/mne/stats/regression.py
+++ b/mne/stats/regression.py
@@ -297,9 +297,10 @@ def solver(X, y):
     coefs = solver(X, data.T)
     if coefs.shape[0] != data.shape[0]:
         raise ValueError(
-            "solver output has unexcepted shape. Supply a "
+            f"solver output has unexcepted shape {coefs.shape}. Supply a "
             "function that returns coefficients in the form "
-            "(n_targets, n_features), where targets == channels."
+            "(n_targets, n_features), where "
+            f"n_targets == n_channels == {data.shape[0]}."
         )
 
     # construct Evoked objects to be returned from output
diff --git a/mne/stats/tests/test_regression.py b/mne/stats/tests/test_regression.py
index 388d9f462e1..84c23a33487 100644
--- a/mne/stats/tests/test_regression.py
+++ b/mne/stats/tests/test_regression.py
@@ -139,7 +139,9 @@ def test_continuous_regression_with_overlap():
     from sklearn.linear_model import ridge_regression
 
     def solver(X, y):
-        return ridge_regression(X, y, alpha=0.0, solver="cholesky")
+        # Newer scikit-learn returns 1D array for ridge_regression, so ensure
+        # 2D output
+        return np.atleast_2d(ridge_regression(X, y, alpha=0.0, solver="cholesky"))
 
     assert_allclose(
         effect,