Skip to content

Carry the model's weights and tie handling into ggforest_subgroup() refits - #870

Open
kassambara wants to merge 5 commits into
masterfrom
fix/ggforest-subgroup-weights
Open

Carry the model's weights and tie handling into ggforest_subgroup() refits#870
kassambara wants to merge 5 commits into
masterfrom
fix/ggforest-subgroup-weights

Conversation

@kassambara

@kassambara kassambara commented Jul 29, 2026

Copy link
Copy Markdown
Owner

ggforest_subgroup() gets each subgroup hazard ratio by refitting the Cox model on
that subset, but the refit was a plain coxph(): the case weights, tie handling and
clustering of the supplied model were dropped. A weighted model reported unweighted
subgroup estimates, and a clustered model drew intervals from the naive variance.

Those settings are now carried into every refit. Weights and clustering ride along as
data columns so subsetting the data subsets them too, and both are read from the fitted
object rather than re-evaluated from its call.

The interaction test is affected the same way, and it cannot simply be refit weighted:
coxph() attaches a robust variance to any non-integer weighted fit, and anova()
declines to compare likelihoods for such a fit because their differences are no longer
chi-square. A fit carrying a robust variance now gets a Wald chi-square on the
treatment-by-subgroup coefficients instead; everything else, integer weights included,
keeps the likelihood-ratio test.

A data whose rows cannot be matched to the fit is now refused when weights or
clustering have to be carried across.

library(survminer)
library(survival)

cc <- colon[colon$etype == 2 & colon$rx %in% c("Obs", "Lev+5FU"), ]
cc$rx    <- droplevels(cc$rx)
cc$sex   <- factor(cc$sex, labels = c("Female", "Male"))
cc$node4 <- factor(cc$node4, labels = c("<= 4 nodes", "> 4 nodes"))

set.seed(9)
cc$w <- runif(nrow(cc), 0.2, 3)          # inverse-probability weights

m <- coxph(Surv(time, status) ~ rx + age, data = cc, weights = w)

ggforest_subgroup(m, data = cc, treatment = "rx",
                  subgroups = c(Sex = "sex", Nodes = "node4"))

Before — the refits ignore the weights:

Subgroup forest plot with overall HR 0.69 and Female 0.87 (0.63-1.19)

After — the overall row reproduces exp(coef(m))["rxLev+5FU"] = 0.70, and the female
subgroup moves from an apparent 13% reduction to exactly no effect:

The same plot with overall HR 0.70 and Female 1.00 (0.70-1.44)

Models specifying none of weights, ties, robust or clustering are unchanged.

…efits

Every subgroup hazard ratio in ggforest_subgroup() comes from refitting the
Cox model on that subset, but the refit was a plain coxph() on the subset:
case weights and the tie handling of the supplied model were both dropped.
A weighted model therefore reported unweighted subgroup estimates, with
nothing on the plot to say so -- on a weighted colon fit the overall hazard
ratio read 0.689 where the model's own value is 0.703, and the female
subgroup read 0.863 against 0.991.

The weights now ride along as a data column, so subsetting the data subsets
them too, and ties/robust are passed through to each refit.

The interaction test is affected the same way and was reporting p = 0.042
from unweighted fits where the weighted likelihood-ratio test gives 8e-05.
Its two comparison models are now fit the same way as the rest, minus the
robust variance coxph() attaches to any weighted fit -- the likelihood-ratio
test compares log-likelihoods, which the variance estimator does not change,
and anova() declines a fit that carries one.

Unweighted models are unaffected.
Following up the weights work: forcing the two comparison models to fit
without a robust variance made anova() run, but routed around a guard that
exists for a reason. ?anova.coxph declines a robust fit because differences
in the log-partial-likelihood no longer have a chi-square distribution, and
that is exactly the case here -- under inverse-probability weights the
model-based information treats the sample as if it held the sum of the
weights. Simulated under a true null with IPT weights, the likelihood-ratio
p-value fired at 24% of replicates at the 5% level; a Wald chi-square on the
interaction coefficients taken from the robust variance fired at 4.9%.

So the interaction is now tested that way whenever the fit carries a robust
variance, and by the likelihood-ratio test otherwise. Integer (frequency)
weights leave the variance model-based and keep the likelihood-ratio test.

Also in this pass:

- Carry the clustering of the model into the refits. coxph() moves a
  cluster() term out of the formula and into the call, so it was being
  dropped along with the weights; a clustered model was drawing intervals
  from the naive variance, which can be several times too narrow.

- Read the weights from the fitted object instead of re-evaluating the
  weights expression, and evaluate ties/robust in the model's own
  environment. Re-evaluating resolved names wherever this function happened
  to be standing: a model fitted inside a wrapper function errored outright,
  and an unrelated object of the same name on the search path was silently
  used in place of the model's real setting.

- Refuse a `data` whose rows cannot be matched to the fit when weights or
  clustering have to be carried across, rather than attaching a weight to
  the wrong subject or falling back to an unweighted plot.

Models that specify none of weights, ties, robust or clustering are
unchanged.
A frequency weight of w means w copies of that row, so replicating each row
and fitting unweighted must reproduce the weighted fit exactly. Adds that as
a test of the subgroup estimates, which needs no second implementation to be
believed. Also merges a comment that had been left duplicated.
Four fixes to the subgroup refits.

`coxph()` records a `cluster()` term and a `cluster` argument the same way, but
`id = ` is stored separately, so an id-clustered model had its clustering
dropped: intervals came from the naive variance, and where `robust = TRUE` was
also given the refits failed outright ("one of cluster or id is needed") and
the plot came back with labels and no hazard ratios at all. `id` is now carried
alongside `cluster`.

The row-matching guard compared the stored response with a rebuilt one at zero
tolerance, but coxph's timefix snaps near-tied times, so the two are not
bit-identical even for the very frame the model was fitted on. On continuous
times that rejected roughly 7% of datasets with an error telling the user to
supply data they had already supplied. It now compares at the ordinary
tolerance, which is what it needs: it is looking for rows in a different order,
not for floating-point noise.

A model fitted with `y = FALSE` stores no response, so nothing can be checked
against. That was treated as "fine" and let mismatched weights through
silently; it is now refused like any other unverifiable case.

An empty factor level, or one present in a single treatment arm, leaves an
aliased interaction coefficient. The Wald test returned NA for the whole term,
which dropped the interaction column from the figure entirely. It now tests the
estimable coefficients, giving the same p-value as dropping the unused level
by hand.

The documentation now lists exactly what the refits carry and says plainly what
they do not -- `subset`, `offset()`, `control` and `tt()` -- notes that the
integrality rule coxph uses to attach a robust variance is about the numbers
rather than what the weights mean, and warns that a Wald interaction test with
few clusters is approximate.
`ties` has an alias, `method`, and coxph records whichever name the user
typed. Reading the call therefore missed the `method =` spelling completely,
so a model fitted `method = "breslow"` was refit with Efron and the plot
disagreed with the model it came from -- 0.688 against the model's 0.702 on
tied data. The resolved value is stored on the fit, so it is now taken from
there and both spellings behave the same.

Whether coxph attaches a robust variance depends on the weights being whole
numbers, which does not distinguish a frequency count from a sampling weight
that happens to be integral. A likelihood-ratio test is only valid for counts:
on integral sampling weights it rejected a true null 17% of the time at the 5%
level. Every weighted fit is now tested with a Wald chi-square on a robust
variance, which is valid either way -- 5.7% on the same simulation.

The row-matching check compared times at a relative tolerance, but timefix
snaps on an absolute budget proportional to the time range, so on data whose
times start near zero and spread widely it still refused the frame the model
was fitted on -- 30% of datasets at n = 10000. Times are now compared on their
own scale and the status exactly.

Two things that were silent are now reported: a fit with fewer than fifty
clusters, whose cluster-robust intervals are too narrow and whose interaction
test rejects too often, and a whole-data refit that fails to reproduce the
model's own treatment coefficient -- which catches a `subset`, an `offset()` or
`control` settings, none of which can be carried into a subset fit.

The documentation now says which parts of a fit travel with the formula and
which are lost, that a weighted or clustered model has to keep its response,
and that rows the model dropped for missing covariates stay out of the subset
fits.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant