Few clusters: the within-cluster design

When splitting whole clusters leaves too few to confirm on

The cluster argument of discover_confirm() assigns whole clusters to one half of the split or the other. That keeps the two halves independent at the cluster level, but the confirmation is left with half the clusters, and its reference degrees of freedom are that number minus one. With eight cohorts, the confirmation would rest on four, three degrees of freedom. The package warns below twenty clusters in the confirmation half, which means about forty in total.

The within-cluster design is the alternative for data nested in a small number of clusters:

  1. Split within clusters. Every cluster is cut in half at random (strata = cluster), so every cluster contributes participants to both halves and the confirmation keeps them all.
  2. Discover on the centered outcome. Before the tree is grown, the outcome is centered on its cluster mean within each imputed dataset (center_by = cluster), so that cluster-level variables cannot stand in for differences between clusters.
  3. Confirm three ways. The standard cluster-robust confirmation on the raw outcome; the same confirmation with cluster fixed effects; and a separate confirmation within each cluster.

discover_confirm(cluster = , within = TRUE) runs the sequence.

Why the extra steps are needed

Because both halves contain the same clusters, a split that merely sorts clusters with different outcome levels can replicate for the wrong reason. Suppose cohorts differ in their average outcome, and some predictor tracks cohort membership, a site characteristic, say. A tree grown on the raw outcome will tend to split on that predictor, because it predicts the outcome through the cohort differences; and the split will then pass a cluster-robust confirmation, because the same cohorts, with the same differences, sit in the confirmation half.

Centering the outcome at discovery removes the cluster differences before the tree sees them. The fixed-effects confirmation removes them again at confirmation, so each split’s contrast is estimated from variation within clusters only. The per-cluster confirmations show whether a split points the same way in most clusters or rests on one.

A worked example

Eight cohorts with their own intercepts, a real within-cohort effect of x1, and a cohort proxy x4 that tracks the cohort intercepts but has no within-cohort relation to the outcome.

library(ctreeMI)
G <- 8; n_per <- 100; n <- G * n_per
cohort <- factor(rep(sprintf("c%02d", 1:G), each = n_per))
u  <- rnorm(G, 0, 0.6)[as.integer(cohort)]
x1 <- rnorm(n); x2 <- rnorm(n); x4 <- u + rnorm(n, 0, 0.3)
y  <- u + 1.2 * (x1 > 0) + rnorm(n)
d  <- data.frame(y, x1, x2, x4, cohort)
d$x1[sample(n, 120)] <- NA
res <- discover_confirm(y ~ x1 + x2 + x4, data = d, m = 5,
                        cluster = "cohort", within = TRUE, seed = 7)
#> Warning: Only 8 clusters in the confirmation data. Cluster-robust variance is
#> unreliable with few clusters, and the reference distributions, whose degrees of
#> freedom are bounded by the number of clusters, will have little power. Treat
#> p-values as approximate and rely on the contrasts and their intervals.
#> Warning: Only 8 clusters in the confirmation data. Cluster-robust variance is
#> unreliable with few clusters, and the reference distributions, whose degrees of
#> freedom are bounded by the number of clusters, will have little power. Treat
#> p-values as approximate and rely on the contrasts and their intervals.
res$within
#> Confirmation within clusters, three ways
#>   8 clusters (cohort), all present in the confirmation data
#>   robust: cluster-robust variance on the raw outcome (report these estimates)
#>   within: the same test with cluster fixed effects
#>   by_cluster: each cluster confirmed on its own
#> 
#>  node split outcome contrast    p_adj within p_adj_within same_sign
#>     1    x1       y    0.914 0.000446  0.925     0.000199       8/8
#>                    reading
#>  confirmed within clusters
#> 
#> 'confirmed within clusters': adjusted p < 0.05 in both the robust and the within test.
#> 'between-cluster difference': robust test passed, within-cluster test did not.
#> same_sign: clusters whose own contrast shares the robust contrast's sign / clusters tested.

The tree split on x1, not on the proxy. The summary lines the three confirmations up by split:

The full confirmation objects are there for inspection: res$within$robust, res$within$within, and res$within$by_cluster, a list with one confirmation per cohort. res$confirmation is the robust one, so prune_unconfirmed() and report_confirm() work as they do for any other analysis.

report_confirm(res)
#> The sample of 800 observations was divided at random within each cohort
#> into a discovery set of 400 and a confirmation set of 400 observations, so
#> that all 8 cohorts contributed to both. Each was multiply imputed
#> separately, with M = 5 imputations, so that neither set's outcomes informed
#> the other's imputed predictor values. A conditional inference tree was
#> fitted to the stacked discovery imputations, with the outcome centered on
#> its cohort mean so that cohort-level differences could not drive the
#> splits, with the Stack/M correction at alpha = 0.05, yielding 2 terminal
#> nodes defined by 1 split. Because the node-level test of that procedure is
#> not calibrated to its nominal level under outcome-conditioned imputation,
#> and because testing a partition on the data used to select it is invalid in
#> any case, each split was then tested on the confirmation set. The
#> confirmation observations falling within each internal node were divided by
#> that node's own rule and the resulting children compared, pooling across
#> imputations by the method of Li, Raghunathan and Rubin (1991) with variance
#> estimated robust to clustering by cohort and adjusting across splits by
#> Holm's method. 1 of 1 split met the 0.05 level on the confirmation data.
#> Because the two sets shared cohorts, each split was also tested with cohort
#> fixed effects, so that its contrast reflected variation within cohorts
#> only, and within each of the 8 cohorts separately. 1 of 1 split met the
#> 0.05 level with fixed effects as well.

What the claim is

With a whole-cluster split, a confirmed split is evidence about clusters the tree never saw. With the within-cluster design, the claim is replication in held-out participants from these clusters. How far that extends to other clusters depends on how consistent the effect is across them, which is what the per-cluster confirmations show: an effect that points the same way in most clusters is unlikely to be an artifact of one or two.

What the design has been checked against

The design was examined in a simulation with 5 to 20 cohorts, equal and unequal cohort sizes, and null, moderate and weak within-cohort effects, against the same split with the tree grown on the raw outcome and against splitting whole cohorts. In short: growing the tree on the centered outcome kept a cohort-level proxy out of the tree, nothing was confirmed within clusters under the null, the fixed-effects test cost almost no power, and unequal cohort sizes changed none of this. The script and results are in the package’s GitHub repository, under simulations/. Effects that differ in size between cohorts were not examined.

Doing it by hand

The pieces are separate functions, for when the imputation model needs control.

parts <- split_holdout(d, strata = "cohort", seed = 42)
imp_d <- mice::mice(parts$discover, m = 30, seed = 1)   # cohort stays in
imp_c <- mice::mice(parts$confirm,  m = 30, seed = 2)   # the imputation model
tree  <- ctree_stacked(y ~ x1 + x2 + x4, data = imp_d, center_by = "cohort")
cw    <- confirm_within(tree, imp_c, cluster = "cohort")

center_outcome() is the centering step on its own; it scores an ordinal outcome and expands a nominal one into centered indicators, and ctree_stacked() grows a multivariate tree on the expansion while keeping the original formula for the confirmation.