Discover, then confirm

The problem this workflow solves

A conditional inference tree chooses its splits by searching the data for the variables and cut points that separate the outcome best. Testing those same splits on the same data is circular: the partition was selected to maximize exactly the separation being tested. With multiply imputed data there is a second problem, specific to this package. The node-level test of ctree_stacked() is not calibrated to its nominal level when the imputation model conditions on the outcome, which is the recommended practice. Under a true null it rejects more often than alpha implies.

Both problems have one remedy. Treat the tree as a hypothesis about subgroup structure, and test that hypothesis on data the tree has never seen. discover_confirm() does the whole sequence: split the data, impute each half separately, grow the tree on one half, and test every split on the other.

A worked example

Six hundred observations, two real predictors and one noise predictor, 15% of x1 and 10% of x2 missing at random.

library(ctreeMI)
n <- 600
d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
d$y <- rnorm(n) + 1.2 * (d$x1 > 0) + 0.9 * (d$x1 > 0) * (d$x2 > 0.5)
d$x1[sample(n, 90)] <- NA
d$x2[sample(n, 60)] <- NA

One call runs the workflow. m = 5 keeps this vignette quick; use 20 to 30 in practice.

res <- discover_confirm(y ~ x1 + x2 + x3, data = d, m = 5, seed = 11)
res
#> Discover-then-confirm Stack/M analysis
#>   discovery n = 300, confirmation n = 300
#> 
#> Discovery tree:
#> -- ctreeMI: Stacked-Imputation Conditional Inference Tree --
#>   Imputations (M)      : 5
#>   Rows per imputation  : 300
#>   Stacked rows         : 1500
#>   Alpha                : 0.0500
#>   Correction           : statistic/M
#>   Outcome dimension    : 1
#>   Degrees of freedom   : derived per node and candidate variable
#>   Splits kept          : 3 of 11
#>   Formula              : y ~ x1 + x2 + x3
#> ------------------------------------------------------------
#> 
#> 
#> Model formula:
#> y ~ x1 + x2 + x3
#> 
#> Fitted party:
#> [1] root
#> |   [2] x1 <= -0.025: 0.047 (n = 706, err = 676.5)
#> |   [3] x1 > -0.025
#> |   |   [4] x2 <= 0.48539
#> |   |   |   [5] x2 <= -0.45439: 0.848 (n = 278, err = 245.1)
#> |   |   |   [6] x2 > -0.45439: 1.420 (n = 301, err = 277.2)
#> |   |   [7] x2 > 0.48539: 2.102 (n = 215, err = 291.7)
#> 
#> Number of inner nodes:    3
#> Number of terminal nodes: 4
#> 
#> Confirmation of a Stack/M partition on independent data
#>   300 confirmation observations, 5 imputations, pooled by Rubin's rules
#> 
#> Omnibus test that the outcome differs across terminal nodes (Li-Raghunathan-Rubin D1):
#>  outcome      F df1   df2      p   riv
#>        y 54.112   3 192.9 <2e-16 0.112
#> 
#> Per-split tests (children of each internal node compared within it):
#>  node depth           split outcome             difference    p_adj   n
#>     1     1 x1 <= -0.024998       y   1.638 [1.378, 1.899]  < 2e-16 300
#>     3     2   x2 <= 0.48539       y   0.764 [0.350, 1.178] 0.000682 131
#>     4     3  x2 <= -0.45439       y -0.046 [-0.519, 0.426] 0.845297  85
#>   difference: pooled contrast between children with 95% interval; p_adj: holm-adjusted across splits.
#>   A split beneath a non-confirmed split is not interpretable.
#>   scale: difference in means (continuous).
#> 
#> Pooled node estimates:
#>  node_id outcome estimate    se n_confirm tested
#>        2       y   -0.076 0.087     168.6   TRUE
#>        5       y    1.316 0.172      44.2   TRUE
#>        6       y    1.270 0.179      40.8   TRUE
#>        7       y    2.057 0.171      46.4   TRUE

Four things happened:

split_holdout() divided the rows at random into a discovery half and a confirmation half, before any imputation.
mice::mice() imputed each half separately, so that neither half’s outcomes informed the other’s imputed predictor values. That separation is what makes the confirmation independent.
ctree_stacked() stacked the discovery imputations and grew one tree with the Stack/M correction.
confirm_ctreeMI() tested every split of that tree on the confirmation half.

What the confirmation tests

Each split in the discovery tree is a hypothesis. For every internal node, the confirmation observations that reach that node are divided by the node’s own rule, the rule found in discovery, and the outcome is compared between the two children. Because the parent and the rule were fixed before the confirmation data were seen, this is a pre-specified comparison. A split whose children do not differ in independent data is one the tree should not have made.

The comparison is a Wald test of the group difference: the model is fitted in each imputed confirmation dataset, the estimates and their covariance are pooled by Rubin’s rules, and the pooled estimate is tested once, by the multivariate Wald procedure of Li, Raghunathan and Rubin (1991) with small-sample degrees of freedom. The p-values are then adjusted across the splits of the tree, by Holm’s method by default.

Every split is tested on its own, on the participants who reach it, whether or not its parent passed. A split beneath an unconfirmed split has no clear interpretation, which is why the depth column is reported.

The output leads with estimates. Each split comes with the pooled difference between its children and a 95% interval; the adjusted p-value says only whether the split survived.

res$confirmation$splits[, c("node_id", "depth", "rule", "contrast", "lower", "upper", "p_adj")]
#>   node_id depth            rule   contrast      lower     upper        p_adj
#> 1       1     1 x1 <= -0.024998  1.6384314  1.3781224 1.8987404 2.984814e-28
#> 2       3     2   x2 <= 0.48539  0.7642759  0.3502183 1.1783334 6.824860e-04
#> 3       4     3  x2 <= -0.45439 -0.0460490 -0.5185912 0.4264932 8.452966e-01

Acting on the result

prune_unconfirmed() collapses every split that did not confirm, working from the bottom of the tree up, and returns the partition the confirmation data support.

pruned <- prune_unconfirmed(res$tree, res$confirmation)
pruned
#> -- ctreeMI: Stacked-Imputation Conditional Inference Tree --
#>   Imputations (M)      : 5
#>   Rows per imputation  : 300
#>   Stacked rows         : 1500
#>   Alpha                : 0.0500
#>   Correction           : statistic/M
#>   Outcome dimension    : 1
#>   Degrees of freedom   : derived per node and candidate variable
#>   Splits kept          : 3 of 11
#>   Formula              : y ~ x1 + x2 + x3
#> ------------------------------------------------------------
#> 
#> 
#> Model formula:
#> y ~ x1 + x2 + x3
#> 
#> Fitted party:
#> [1] root
#> |   [2] x1 <= -0.025: 0.047 (n = 706, err = 676.5)
#> |   [3] x1 > -0.025
#> |   |   [4] x2 <= 0.48539: 1.146 (n = 579, err = 569.5)
#> |   |   [5] x2 > 0.48539: 2.102 (n = 215, err = 291.7)
#> 
#> Number of inner nodes:    2
#> Number of terminal nodes: 3

report_confirm() writes the methods paragraph: how the sample was divided, how each half was imputed, what was discovered, and what survived.

report_confirm(res)
#> The sample of 600 observations was divided at random into a discovery set
#> of 300 and a confirmation set of 300. 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 Stack/M correction at alpha = 0.05,
#> yielding 4 terminal nodes defined by 3 splits. 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) and adjusting across splits by Holm's method.
#> 2 of 3 splits met the 0.05 level on the confirmation data; the remainder
#> was not supported and the corresponding node was collapsed.

Outcome types

The confirmation handles every outcome a tree can be grown on.

Outcome	Model in the confirmation	The contrast reported for a split
numeric	linear regression	difference in means
two-level factor or logical	logistic regression	log odds ratio
ordered factor, 3+ levels	proportional-odds regression (`MASS::polr`)	log cumulative odds ratio
unordered factor, 3+ levels	multinomial regression (`nnet::multinom`)	log odds ratio against the first level, for the level that differs most; the test covers all levels

Several outcomes can be confirmed together, as in y1 + y2 ~ x1 + x2; each split is tested for each outcome. For an ordinal or nominal outcome, node estimates are the proportions in each level, one row per node and level.

outcome_type = "auto" resolves the type from the column. It can be set explicitly, one entry per outcome, when a numeric 0/1 column should be treated as binary, for example.

d$y_ord <- cut(d$y, c(-Inf, 0, 1.5, Inf), labels = c("low", "mid", "high"),
               ordered_result = TRUE)
res_ord <- discover_confirm(y_ord ~ x1 + x2 + x3, data = d, m = 5, seed = 12)
res_ord$confirmation$splits[, c("node_id", "rule", "contrast", "p_adj")]
#>   node_id            rule  contrast        p_adj
#> 1       1 x1 <= -0.022559 2.1016531 9.699531e-05
#> 2       3    x2 <= 1.0489 0.7891126 1.103380e-01
head(res_ord$confirmation$nodes)
#>   node_id outcome level  estimate         se n_confirm n_clusters tested
#> 1       2   y_ord   low 0.5200486 0.03759153     161.6         NA   TRUE
#> 2       2   y_ord   mid 0.3776375 0.04054167     161.6         NA   TRUE
#> 3       2   y_ord  high 0.1023138 0.03935261     161.6         NA   TRUE
#> 4       4   y_ord   low 0.1156758 0.04277880     115.8         NA   TRUE
#> 5       4   y_ord   mid 0.4771262 0.05115256     115.8         NA   TRUE
#> 6       4   y_ord  high 0.4071980 0.04434025     115.8         NA   TRUE

Clustered data

When observations are nested in cohorts, sites or classrooms, two things must change: the split should keep each cluster on one side, and the confirmation should use a variance that respects the nesting. Both are the cluster argument.

res <- discover_confirm(y ~ x1 + x2 + x3, data = d, cluster = "cohort",
                        m = 30, seed = 1)

With cluster, whole clusters go to one half or the other, so a confirmed split is evidence about clusters the tree never saw. The cost is that only the clusters in the confirmation half are available for testing: the reference degrees of freedom are that number minus one, and the package warns below twenty. When there are too few clusters to split whole, see the “Few clusters” vignette for the within-cluster design.