Cluster-aware bootstrap resampling for repeated measures designs

CHANGELOG.md

@@ -2,6 +2,17 @@
 
 <!-- do not remove -->
 
+## Unreleased
+
+### New Features
+1. **Cluster-aware bootstrap and permutation tests**: `dabest.load()` accepts a new `cluster_col` argument naming the column that identifies the independent sampling unit (cluster) each observation belongs to, such as a participant who contributes several observations or several pairs of paired observations. When it is set, the bootstrap resamples whole clusters with replacement (a cluster bootstrap, stratified by the pattern of groups each cluster appears in) and the permutation test reshuffles labels at the cluster level, so that confidence intervals and permutation p-values account for the correlation between observations from the same cluster. This works for unpaired data, for paired data (`paired` with `id_col`, where `id_col` identifies the pairs and `cluster_col` the units the pairs are nested in), for shared-control and multi-group `idx`, and for delta-delta and mini-meta analyses. The results table gains an `n_clusters` column and `TwoGroupsEffectSize`/`PermutationTest` accept `control_clusters`/`test_clusters` directly. The parametric and rank-based tests in `statistical_tests` are unchanged and still ignore clustering.
+2. **Cluster count on plots**: when `cluster_col` is set, each group's axis label reports the number of distinct clusters on a second line below the number of observations, `(N=<observations>,` / ` n=<clusters>)` instead of the usual `(N=<observations>)`, since the cluster count is what the bootstrap actually resamples. To make room, the gap between the raw-data and contrast axes of vertical Cumming plots is then sized to the height of the labels (including the taller labels of two-column Sankey plots), so they clear the contrast axes. Unclustered plots are unchanged.
+3. **Small-sample expansion of cluster-bootstrap intervals**: bootstrap intervals are too narrow when there are few independent units, so with `cluster_col` set, confidence intervals are now expanded for the number of clusters by default, using the expanded percentile method of Hesterberg (2015). Each interval is read further into the tails of the same bootstrap distribution, at a level derived from a t distribution with degrees of freedom based on the number of clusters (for nested or mixed designs, whose groups of clusters are resampled independently, combined with the Welch-Satterthwaite approximation). This applies to the percentile and bias-corrected and accelerated intervals, and to those of the baseline error curve, delta-delta and mini-meta analyses. The reported confidence level (`ci`) is unchanged, and the level actually read is reported as `ci_expanded` in the results and in the printed summary, with the unexpanded limits alongside (`bca_low_unexpanded` and so on). Because the plotted bootstrap distribution is unchanged, estimation plots and forest plots draw an expanded interval in two parts: the unexpanded interval as the usual thick bar, and the expansion beyond it as a thinner line, which can be styled with `contrast_expanded_errorbar_kwargs` in `.plot()` (or `expanded_errorbar_kwargs` in `forest_plot()`). In simulations, expanded intervals came within about 3 percentage points of nominal coverage for every effect size once there were at least 6 clusters, and at least 4 in each independently resampled group of clusters; below that, a warning is given. Pass `cluster_ci_expansion=False` to `dabest.load()` to turn the expansion off. `cluster_col` can also no longer be the same column as `x`.
+
+### Documentation
+1. **Baseline error curve, explained**: the [Plot Aesthetics tutorial](nbs/tutorials/08-plot_aesthetics.ipynb) and the `show_baseline_ec` docstring now spell out what the baseline error curve (`show_baseline_ec=True`) actually computes, and call out that it is always an *unpaired* self-comparison of the control group, regardless of `paired`. This matters with `cluster_col`: paired real comparisons largely cancel between-cluster variation, but the always-unpaired baseline curve does not, so it can become much wider than the real contrasts once clustering is on. A worked example with and without `cluster_col` is included.
+2. **New tutorial: [Cluster-Robust Bootstrap for Repeated Measures](nbs/tutorials/11-cluster_robust_bootstrap.ipynb)**: a self-contained, simulation-based worked example of `cluster_col` for the common case of a participant contributing several sets of paired observations. Its central example is tuned so that, for the identical data and point estimate, the naive dummy-ID bootstrap's 95% interval lies entirely above zero while the cluster-aware interval spans it, making the practical stakes of pseudoreplication concrete rather than abstract. A 200-dataset coverage simulation then shows this is systematic, not a fluke of one dataset: against a known true effect, the naive interval covers the truth only about 82% of the time, the unexpanded cluster-aware interval about 90%, and the default expanded cluster-aware interval about 95%. A section on the small-sample expansion explains when and why it applies.
+
 ## v2025.10.20
 
 ### New Features

dabest/_api.py

@@ -25,6 +25,8 @@ def load(
     x1_level=None,
     mini_meta=False,
     ps_adjust=False,
+    cluster_col=None,
+    cluster_ci_expansion=True,
 ):
     """
     Loads data in preparation for estimation statistics.
@@ -88,6 +90,57 @@ def load(
     ps_adjust : boolean, default False
         Indicator of whether to adjust calculated p-value according to Phipson & Smyth (2010)
         # https://doi.org/10.2202/1544-6115.1585
+    cluster_col : string, default None
+        Name of the column identifying the independent sampling unit (cluster)
+        that each observation belongs to, for example a participant who
+        contributes several observations, or several pairs of paired
+        observations. When supplied, the bootstrap resamples whole clusters
+        with replacement (a cluster bootstrap) and the permutation test
+        reshuffles labels at the cluster level (swapping the control and test
+        observations of whole clusters, or reassigning whole clusters between
+        groups when clusters are nested within groups), so that the confidence
+        intervals and permutation p-values account for the correlation between
+        observations from the same cluster. This works with both unpaired data
+        and paired data (`paired` with `id_col`): for paired data, `id_col`
+        identifies the pairs and `cluster_col` the units the pairs are nested
+        in, and every pair must belong to a single cluster. The parametric and
+        rank-based tests reported in `statistical_tests` do not account for
+        clustering. Note that the baseline error curve shown by
+        `.plot(show_baseline_ec=True)` is always an unpaired comparison (see
+        that argument's docstring), so with `cluster_col` set it can become
+        much wider than a paired analysis's real effect-size curves.
+    cluster_ci_expansion : boolean, default True
+        Only used when `cluster_col` is set. Bootstrap confidence intervals are
+        too narrow when there are few independent units: with around 15 to 30
+        clusters, a nominal 95% cluster-bootstrap interval typically covers the
+        true effect only 90 to 94% of the time. When True, the percentile and
+        bias-corrected and accelerated intervals (including those of the
+        baseline error curve, delta-delta and mini-meta analyses) are therefore
+        expanded for the number of clusters with the expanded percentile
+        method of Hesterberg (2015, The American Statistician, 69(4),
+        371-386): each interval is read further into the tails of the same
+        bootstrap distribution, at a level derived from a t distribution with
+        degrees of freedom based on the number of clusters. When clusters are
+        nested within groups, or some clusters appear in only some groups, the
+        groups of clusters are resampled independently, and their degrees of
+        freedom are combined with the Welch-Satterthwaite approximation. The
+        reported confidence level (`ci`) is unchanged; the level actually read
+        is reported as `ci_expanded` in the results, alongside the unexpanded
+        limits (`bca_low_unexpanded` and so on). Plots draw the unexpanded
+        interval (the nominal `ci`% interval of the plotted bootstrap
+        distribution) as the usual thick bar, and the expansion beyond it as a
+        thinner line (see `contrast_expanded_errorbar_kwargs` in `plot()`). The
+        correction fades as the number of clusters grows. In simulations the
+        expanded intervals came within about 3 percentage points of nominal
+        coverage for every effect size once there were at least 6 clusters in
+        all and at least 4 in each independently resampled group of clusters
+        (for example 6 participants who each take part in every condition, or
+        8 participants split between two conditions); with fewer, no bootstrap
+        interval is reliable, and a warning is given. Because expanded
+        intervals are read further into the
+        tails of the bootstrap distribution, consider increasing `resamples`
+        (to 20000, say) when there are few clusters. Set to False to report
+        unexpanded cluster-bootstrap intervals.
 
     Returns
     -------
@@ -112,6 +165,8 @@ def load(
         x1_level,
         mini_meta,
         ps_adjust,
+        cluster_col=cluster_col,
+        cluster_ci_expansion=cluster_ci_expansion,
     )
 
 # %% ../nbs/API/load.ipynb #570ff65a

dabest/_dabest_object.py

@@ -41,6 +41,8 @@ def __init__(
         x1_level,
         mini_meta,
         ps_adjust,
+        cluster_col=None,
+        cluster_ci_expansion=True,
     ):
         """
         Parses and stores pandas DataFrames in preparation for estimation
@@ -60,6 +62,8 @@ def __init__(
         self.__is_proportional = proportional
         self.__is_mini_meta = mini_meta
         self.__ps_adjust = ps_adjust
+        self.__cluster_col = cluster_col
+        self.__cluster_ci_expansion = cluster_ci_expansion
 
         # after this call the attributes self.__experiment_label and self.__x1_level are updated
         self._check_errors(x, y, idx, experiment, experiment_label, x1_level)
@@ -128,6 +132,16 @@ def __repr__(self):
         resamples_line2 = "will be used to generate the effect size bootstraps."
         out.append(resamples_line1 + resamples_line2)
 
+        if self.__cluster_col is not None:
+            cluster_line1 = "Whole clusters, as defined by `{}`, ".format(self.__cluster_col)
+            cluster_line2 = "will be resampled by the bootstrap and reshuffled by the permutation test."
+            out.append(cluster_line1 + cluster_line2)
+            if self.__cluster_ci_expansion:
+                out.append(
+                    "Confidence intervals will be expanded for the number of clusters "
+                    "(set `cluster_ci_expansion=False` to turn this off)."
+                )
+
         return "\n".join(out)
 
 
@@ -341,6 +355,23 @@ def id_col(self):
         """
         return self.__id_col
 
+    @property
+    def cluster_col(self):
+        """
+        Returns the cluster column declared to `dabest.load()`, if any.
+        When set, the bootstrap resamples whole clusters of observations and
+        the permutation test reshuffles labels at the cluster level.
+        """
+        return self.__cluster_col
+
+    @property
+    def cluster_ci_expansion(self):
+        """
+        Returns whether confidence intervals of clustered data are expanded for
+        a small number of clusters, as declared to `dabest.load()`.
+        """
+        return self.__cluster_ci_expansion
+
     @property
     def ci(self):
         """
@@ -582,6 +613,33 @@ def _check_errors(self, x, y, idx, experiment, experiment_label, x1_level):
             if self.__id_col not in self.__output_data.columns:
                 err = "`id_col` was given as '{}'; however, '{}' is not a column in `data`.".format(self.__id_col, self.__id_col)
                 raise IndexError(err)
+
+        # Check if `cluster_ci_expansion` is valid
+        if not isinstance(self.__cluster_ci_expansion, (bool, np.bool_)):
+            raise TypeError("`cluster_ci_expansion` must be True or False.")
+
+        # Check if `cluster_col` is valid
+        if self.__cluster_col is not None:
+            if self.__cluster_col not in self.__output_data.columns:
+                err = "`cluster_col` was given as '{}'; however, '{}' is not a column in `data`.".format(self.__cluster_col, self.__cluster_col)
+                raise IndexError(err)
+
+            if y is not None and self.__cluster_col == y:
+                err = "`cluster_col` cannot be the same column as `y`."
+                raise ValueError(err)
+
+            x_columns = x if isinstance(x, (list, tuple)) else [x]
+            if self.__cluster_col in x_columns:
+                err1 = "`cluster_col` cannot be the same column as `x`: every group would "
+                err2 = "then be a single cluster, leaving nothing to resample."
+                raise ValueError(err1 + err2)
+
+            if x is None and idx is not None:
+                # Wide format: the cluster column cannot also be one of the groups.
+                groups = [g for item in idx for g in (item if isinstance(item, (tuple, list)) else (item,))]
+                if self.__cluster_col in groups:
+                    err = "`cluster_col` ('{}') cannot also be one of the groups in `idx`.".format(self.__cluster_col)
+                    raise ValueError(err)
 
         # Check if x and y are supplied (relevant to long format data)
         if x is None and y is not None:
@@ -677,8 +735,32 @@ def _get_plot_data(self, x, y, all_plot_groups):
                 plot_data[self.__xvar], categories=all_plot_groups, ordered=True
             )
 
+        if self.__cluster_col is not None:
+            self._check_clusters(plot_data)
+
         return plot_data
 
+    def _check_clusters(self, plot_data):
+        """
+        Check that the cluster labels are complete and, for paired data,
+        consistent within each `id_col` value.
+        """
+        clusters = plot_data[self.__cluster_col]
+        if clusters.isnull().any():
+            err1 = "`cluster_col` ('{}') contains missing values.".format(self.__cluster_col)
+            err2 = " Every observation must belong to a cluster."
+            raise ValueError(err1 + err2)
+
+        if self.__is_paired:
+            clusters_per_id = plot_data.groupby(self.__id_col, observed=True)[self.__cluster_col].nunique()
+            inconsistent = clusters_per_id.index[clusters_per_id > 1].tolist()
+            if inconsistent:
+                err1 = "Each value of `id_col` must belong to a single cluster in `cluster_col`,"
+                err2 = " but the following values of '{}' have more than one cluster label: {}.".format(
+                    self.__id_col, inconsistent[:10]
+                )
+                raise ValueError(err1 + err2)
+
     def _compute_effectsize_dfs(self):
         '''
         Function to compute all attributes based on EffectSizeDataFrame.
@@ -698,6 +780,7 @@ def _compute_effectsize_dfs(self):
             x2=self.__x2,
             mini_meta=self.__is_mini_meta,
             ps_adjust=self.__ps_adjust,
+            cluster_ci_expansion=self.__cluster_ci_expansion,
         )
 
         self.__mean_diff = EffectSizeDataFrame(