robin_surv(): factor label collision silently merges distinct joint strata, corrupting estimates

Summary

robin_surv() groups joint strata by pasted factor labels rather than by level identity. When two distinct strata combinations paste to the same string, they silently collapse into one stratum. This affects two independent sites in R/survival_score.R, one of which corrupts the point estimate and standard error.

Because interaction() and split() default to sep = ".", ordinary labels such as "a.b" are enough to trigger this — no unusual characters required. Note that interaction() output itself contains dots, so a pre-computed joint-stratum column is a realistic trigger.

This is pre-existing on main (43ba7c1), not introduced by #114. The equivalent bug in robin_mh() was fixed in #114 via a new internal helper that derives stratum identity from level codes instead of labels; the same approach applies here.

Site 1 — wrong estimate and SE (R/survival_score.R:182)

df_split <- split(df, f = as.formula(strata_formula), drop = TRUE)

split() groups by pasted labels, so distinct strata merge. Same data, only the factor labels differ:

library(RobinCar2)
library(survival)

set.seed(9); n <- 120
d <- data.frame(
  s1  = factor(rep(c("a.b", "a"), each = 2, length.out = n)),
  s2  = factor(rep(c("c", "b.c"), each = n / 2)),
  trt = factor(rep(c("A", "B"), n / 2))
)
rate <- ifelse(d$s1 == "a.b", 0.05, 0.4)   # strong stratum effect
d$time <- rexp(n, rate); d$status <- rbinom(n, 1, 0.8)
f <- Surv(time, status) ~ strata(s1) + strata(s2)

r1 <- robin_surv(f, data = d, treatment = trt ~ pb(s1, s2))

d2 <- d                                     # identical data, renamed levels
levels(d2$s1) <- c("P", "Q"); levels(d2$s2) <- c("R", "S")
r2 <- robin_surv(f, data = d2, treatment = trt ~ pb(s1, s2))

colliding labels : est = 0.2234474041  se = 0.2193840365
renamed levels   : est = 0.2374141150  se = 0.2216864755

The estimate differs by 1.4e-02 and the SE by 2.3e-03 purely because of how the levels are spelled. Root cause:

d <- data.frame(s1 = factor(c("a.b", "a.b", "a", "a")),
                s2 = factor(c("c", "c", "b.c", "b.c")), v = 1:4)
length(split(d, f = ~ s1 + s2, drop = TRUE))   # 1  -- should be 2

Site 2 — randomization-strata warning suppressed (R/survival_score.R:173)

!h_first_fct_nested_in_second(interaction(df[strata]), interaction(df[randomization_strata]))

The collision merges randomization levels, which makes the nesting test spuriously succeed, so a warning that should fire is silently dropped:

set.seed(4); n <- 100
d <- data.frame(
  g   = factor(rep("g1", n)),                    # analysis stratum, coarser
  s1  = factor(rep(c("a.b", "a"), each = n / 2)),
  s2  = factor(rep(c("c", "b.c"), each = n / 2)),
  trt = factor(rep(c("A", "B"), n / 2))
)
d$time <- rexp(n, 0.1); d$status <- rbinom(n, 1, 0.7)

robin_surv(Surv(time, status) ~ strata(g), data = d, treatment = trt ~ pb(s1, s2))

colliding      -> NO WARNING      (wrong: s1, s2 are not covered by g)
non-colliding  -> WARNED (correct)

Minimal check:

g <- factor(rep("g1", 4))
h_first_fct_nested_in_second(g, interaction(data.frame(
  s1 = factor(c("a.b", "a.b", "a", "a")),
  s2 = factor(c("c", "c", "b.c", "b.c"))), drop = TRUE))
#> TRUE   -- should be FALSE

Suggested fix

Derive stratum identity from level codes, not labels — labels are cosmetic. #114 added h_mh_joint_strata() in R/robin_mh.R for exactly this: it builds a mixed-radix index from the per-column integer codes, reproducing interaction()'s level order and labels (made unique via make.unique()) while keeping colliding combinations distinct. It is verified identical() to interaction(..., sep = ":") in both labels and codes across 200 randomized no-collision cases.

To reuse it here it should be renamed (it is not MH-specific) and moved to R/utils.R, then applied at both sites — with split() keyed on the resulting factor rather than on a formula.

Changing the grouping will require reviewing tests/testthat/_snaps/survival.md and re-validating the coxph / RobinCar equivalence tests, so this is deliberately left out of #114.

Note that robin_glm() / robin_lm() are unaffected — their h_interaction() calls are about model formula interactions, unrelated to joint strata.