diff --git a/mne/stats/_adjacency.py b/mne/stats/_adjacency.py
index c81344a8aef..8b0abf5bb67 100644
--- a/mne/stats/_adjacency.py
+++ b/mne/stats/_adjacency.py
@@ -28,6 +28,25 @@ def combine_adjacency(*structure):
     -------
     adjacency : scipy.sparse.coo_matrix, shape (n_features, n_features)
         The adjacency matrix.
+
+    Notes
+    -----
+    For 4-dimensional data with shape ``(n_obs, n_times, n_freqs, n_chans)``,
+    you can specify **no** connections among elements in a particular
+    dimension by passing a matrix of zeros. For example:
+
+    >>> import numpy as np
+    >>> from scipy.sparse import diags
+    >>> from mne.stats import combine_adjacency
+    >>> n_times, n_freqs, n_chans = (50, 7, 16)
+    >>> chan_adj = diags([1., 1.], offsets=(-1, 1), shape=(n_chans, n_chans))
+    >>> combine_adjacency(
+    ...     n_times,  # regular lattice adjacency for times
+    ...     np.zeros((n_freqs, n_freqs)),  # no adjacency between freq. bins
+    ...     chan_adj,  # custom matrix, or use mne.channels.find_ch_adjacency
+    ...     )  # doctest: +NORMALIZE_WHITESPACE
+    <5600x5600 sparse matrix of type '<class 'numpy.float64'>'
+            with 27076 stored elements in COOrdinate format>
     """
     from scipy import sparse
     structure = list(structure)
@@ -35,15 +54,11 @@ def combine_adjacency(*structure):
         name = f'structure[{di}]'
         _validate_type(dim, ('int-like', np.ndarray, sparse.spmatrix), name)
         if isinstance(dim, int_like):
-            dim = int(dim)
-            # Don't add the diagonal, because we explicitly remove it later:
-            # dim = sparse.eye(dim, format='coo')
-            # dim += sparse.eye(dim.shape[0], k=1, format='coo')
-            # dim += sparse.eye(dim.shape[0], k=-1, format='coo')
-            ii, jj = np.arange(0, dim - 1), np.arange(1, dim)
-            edges = np.vstack([np.hstack([ii, jj]), np.hstack([jj, ii])])
-            dim = sparse.coo_matrix(
-                (np.ones(edges.shape[1]), edges), (dim, dim), float)
+            # Don't add the diagonal, because we explicitly remove it later
+            dim = sparse.diags([1, 1],
+                               offsets=(-1, 1),
+                               shape=(int(dim), int(dim)),
+                               dtype=float).tocoo()
         else:
             _check_option(f'{name}.ndim', dim.ndim, [2])
             if dim.shape[0] != dim.shape[1]:
diff --git a/mne/utils/docs.py b/mne/utils/docs.py
index 435fcc3149e..9a4160bb532 100644
--- a/mne/utils/docs.py
+++ b/mne/utils/docs.py
@@ -1782,7 +1782,7 @@
 stat_fun : callable | None
     Function called to calculate the test statistic. Must accept 1D-array as
     input and return a 1D array. If ``None`` (the default), uses
-    :func:`mne.stats.{}`.
+    `mne.stats.{}`.
 """
 docdict['clust_stat_f'] = docdict['clust_stat'].format('f_oneway')
 docdict['clust_stat_t'] = docdict['clust_stat'].format('ttest_1samp_no_p')
@@ -1803,10 +1803,11 @@
        'a square matrix with dimension ``{x}.shape[-1]`` (n_vertices) to save '
        'memory and computation, and to use ``max_step`` to define the extent '
        'of temporal adjacency to consider when clustering.')
+comb = ' The function `mne.stats.combine_adjacency` may be useful for 4D data.'
 st = dict(sp='spatial', lastdim='', parone='(n_vertices)',
           partwo='(n_times * n_vertices)', memory=mem)
 tf = dict(sp='', lastdim=' (or the last two dimensions if ``{x}`` is 2D)',
-          parone='', partwo='', memory='')
+          parone='(for 3D data)', partwo='(for 4D data)', memory=comb)
 nogroups = dict(eachgrp='', x='X')
 groups = dict(eachgrp='each group ', x='X[k]')
 docdict['clust_adj_st1'] = docdict['clust_adj'].format(**st).format(**nogroups)