From 95b952ded3561c1cf91947fd190461f2b2430ec9 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 11:20:27 +0000 Subject: [PATCH 1/4] survFit: support models fitted with case weights Use the weighted Breslow estimator of the cumulative baseline hazard instead of stopping with an error for models fitted with non-unit weights: events and risk sets are weighted by the case weights, and the baseline survivor function is computed at the weighted mean of the linear predictor (identical to the previous behaviour for unit weights). Results agree with survival::survfit for weighted coxph fits, and fits with integer case weights give results identical to fitting on the correspondingly expanded data. Closes boost-R/mboost#121 Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HqB9XX7JjdQJ1vah4eyWv9 --- R/survival.R | 21 ++++++++++------ inst/NEWS.Rd | 7 ++++++ man/survFit.Rd | 10 +++++++- tests/regtest-blackboost.R | 9 +++++++ tests/regtest-blackboost.Rout.save | 9 +++++++ tests/regtest-glmboost.R | 30 +++++++++++++++++++++++ tests/regtest-glmboost.Rout.save | 39 ++++++++++++++++++++++++++++++ 7 files changed, 116 insertions(+), 9 deletions(-) diff --git a/R/survival.R b/R/survival.R index 071eed4c..5a50829e 100644 --- a/R/survival.R +++ b/R/survival.R @@ -5,20 +5,25 @@ survFit <- function(object, ...) survFit.mboost <- function(object, newdata = NULL, ...) { - n <- length(w <- model.weights(object)) - if (!isTRUE(all.equal(w,rep(1,n)))) - stop("survFit cannot (yet) deal with weights") + w <- model.weights(object) y <- object$response stopifnot(inherits(y, "Surv")) + pr <- predict(object) + ### center linear predictor around its weighted mean such that the + ### baseline curve corresponds to the (weighted) "average" observation + lpc <- weighted.mean(pr, w) + ord <- order(y[,1]) time <- y[ord,1] event <- y[ord,2] - n.event <- aggregate(event, by = list(time), sum)[,2] + w <- w[ord] + ### weighted number of events at each time point + n.event <- aggregate(w * event, by = list(time), sum)[,2] - pr <- predict(object) - linpred <- (pr - mean(pr))[ord] - R <- rev(cumsum(rev(exp(linpred)))) + linpred <- (pr - lpc)[ord] + ### weighted Breslow estimator of the cumulative baseline hazard + R <- rev(cumsum(rev(w * exp(linpred)))) R[R<1e-6] <- Inf d <- c(1,diff(time)) != 0 @@ -29,7 +34,7 @@ survFit.mboost <- function(object, newdata = NULL, ...) devent <- n.event > 0 if (!is.null(newdata)){ S <- exp( tcrossprod( -H, exp(as.numeric(predict(object, - newdata=newdata)- mean(predict(object)))) ))[devent,] + newdata=newdata) - lpc)) ))[devent, , drop = FALSE] colnames(S) <- rownames(newdata) } else S <- matrix(exp(-H)[devent], ncol = 1) diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index e1268e9f..f2bc3517 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -2,6 +2,13 @@ \title{News for Package 'mboost'} \section{Changes in mboost version 2.9-13 (2026-07-16)}{ + \subsection{User-visible changes}{ + \itemize{ + \item \code{survFit} can now deal with models fitted with case + weights by using the weighted Breslow estimator + (\url{https://github.com/boost-R/mboost/issues/121}). + } + } \subsection{Other}{ \itemize{ \item Handle M1/ARM precision issues in \file{tests}. diff --git a/man/survFit.Rd b/man/survFit.Rd index c85c7eef..63778250 100644 --- a/man/survFit.Rd +++ b/man/survFit.Rd @@ -30,6 +30,13 @@ of \code{object}, along with the Breslow estimator, is used for computing the predicted survivor function for each row in \code{newdata}. + If the model was fitted with case weights, the weighted Breslow + estimator is used for computing the baseline survivor function, i.e., + each observation contributes to the number of events and to the risk + sets according to its weight, and the survivor function for + \code{newdata = NULL} is computed at the weighted mean of the linear + predictor. + } \value{ An object of class \code{survFit} containing the following components: @@ -38,7 +45,8 @@ \item{time}{ the time points at which the survivor functions are evaluated. } \item{n.event}{ the number of events observed at each time point given - in \code{time}.} + in \code{time}. For models fitted with case weights, + this is the weighted number of events.} } \seealso{ \code{\link{gamboost}}, \code{\link{glmboost}} and \code{\link{blackboost}} for model fitting.} diff --git a/tests/regtest-blackboost.R b/tests/regtest-blackboost.R index 8336e102..8feea669 100644 --- a/tests/regtest-blackboost.R +++ b/tests/regtest-blackboost.R @@ -84,6 +84,15 @@ newdata <- ovarian[c(1,3,12),] A2 <- survFit(fit2, newdata = newdata) print(A2) +## survFit also works for models fitted with case weights +## see https://github.com/boost-R/mboost/issues/121 +w <- seq(0.5, 1.5, length.out = nrow(ovarian)) +fit2w <- blackboost(Surv(futime,fustat) ~ age + resid.ds + rx + ecog.ps, + data = ovarian, family = CoxPH(), weights = w, + control = boost_control(mstop = 100)) +A2w <- survFit(fit2w, newdata = newdata) +stopifnot(all(A2w$surv >= 0), all(A2w$surv <= 1)) + ### predictions: set.seed(1907) x1 <- rnorm(100) diff --git a/tests/regtest-blackboost.Rout.save b/tests/regtest-blackboost.Rout.save index 821dcf1e..73b17604 100644 --- a/tests/regtest-blackboost.Rout.save +++ b/tests/regtest-blackboost.Rout.save @@ -104,6 +104,15 @@ Loading required package: stabs + A2 <- survFit(fit2, newdata = newdata) + print(A2) + ++ ## survFit also works for models fitted with case weights ++ ## see https://github.com/boost-R/mboost/issues/121 ++ w <- seq(0.5, 1.5, length.out = nrow(ovarian)) ++ fit2w <- blackboost(Surv(futime,fustat) ~ age + resid.ds + rx + ecog.ps, ++ data = ovarian, family = CoxPH(), weights = w, ++ control = boost_control(mstop = 100)) ++ A2w <- survFit(fit2w, newdata = newdata) ++ stopifnot(all(A2w$surv >= 0), all(A2w$surv <= 1)) ++ + ### predictions: + set.seed(1907) + x1 <- rnorm(100) diff --git a/tests/regtest-glmboost.R b/tests/regtest-glmboost.R index d8348d1a..a3ef2b79 100644 --- a/tests/regtest-glmboost.R +++ b/tests/regtest-glmboost.R @@ -197,6 +197,36 @@ A3 <- survFit(fit3, newdata = newdata) max(A1$surv-A2$surv) max(A1$surv-A3$surv) +## survFit for models fitted with case weights +## see https://github.com/boost-R/mboost/issues/121 + +## integer case weights: same results as fitting on the expanded data +w <- rep(1:2, length.out = nrow(ovarian)) +fit4 <- glmboost(fm, data = ovarian, family = CoxPH(), weights = w, + control = boost_control(mstop = 1000), center = FALSE) +fit5 <- glmboost(fm, data = ovarian[rep(1:nrow(ovarian), w),], family = CoxPH(), + control = boost_control(mstop = 1000), center = FALSE) +A4 <- survFit(fit4) +A5 <- survFit(fit5) +print(.all.equal(A4$surv, A5$surv)) +print(.all.equal(A4$time, A5$time)) +print(.all.equal(A4$n.event, A5$n.event)) +A4 <- survFit(fit4, newdata = newdata) +A5 <- survFit(fit5, newdata = newdata) +print(.all.equal(A4$surv, A5$surv)) + +## non-integer case weights: compare with weighted coxph fit +wts <- seq(0.5, 1.5, length.out = nrow(ovarian)) +fitw <- coxph(fmSurv, data = ovarian, weights = wts) +fit6 <- glmboost(fm, data = ovarian, family = CoxPH(), weights = wts, + control = boost_control(mstop = 1000), center = TRUE) +A1 <- survfit(fitw, censor = FALSE) +A6 <- survFit(fit6) +print(max(abs(A1$surv - A6$surv)) < 0.001) +A1 <- survfit(fitw, newdata = newdata, censor = FALSE) +A6 <- survFit(fit6, newdata = newdata) +print(max(abs(A1$surv - A6$surv)) < 0.001) + ### Poisson models df <- data.frame(x1 = runif(100), x2 = runif(100)) diff --git a/tests/regtest-glmboost.Rout.save b/tests/regtest-glmboost.Rout.save index 3c9e4cf6..dc1cc011 100644 --- a/tests/regtest-glmboost.Rout.save +++ b/tests/regtest-glmboost.Rout.save @@ -289,6 +289,45 @@ In glmboost.formula(fm, data = ovarian, family = CoxPH(), control = boost_contro > max(A1$surv-A3$surv) [1] 0.189 > +> ## survFit for models fitted with case weights +> ## see https://github.com/boost-R/mboost/issues/121 +> +> ## integer case weights: same results as fitting on the expanded data +> w <- rep(1:2, length.out = nrow(ovarian)) +> fit4 <- glmboost(fm, data = ovarian, family = CoxPH(), weights = w, ++ control = boost_control(mstop = 1000), center = FALSE) +> fit5 <- glmboost(fm, data = ovarian[rep(1:nrow(ovarian), w),], family = CoxPH(), ++ control = boost_control(mstop = 1000), center = FALSE) +> A4 <- survFit(fit4) +> A5 <- survFit(fit5) +> print(.all.equal(A4$surv, A5$surv)) +[1] TRUE +> print(.all.equal(A4$time, A5$time)) +[1] TRUE +> print(.all.equal(A4$n.event, A5$n.event)) +[1] TRUE +> A4 <- survFit(fit4, newdata = newdata) +> A5 <- survFit(fit5, newdata = newdata) +> print(.all.equal(A4$surv, A5$surv)) +[1] TRUE +> +> ## non-integer case weights: compare with weighted coxph fit +> wts <- seq(0.5, 1.5, length.out = nrow(ovarian)) +> fitw <- coxph(fmSurv, data = ovarian, weights = wts) +> fit6 <- glmboost(fm, data = ovarian, family = CoxPH(), weights = wts, ++ control = boost_control(mstop = 1000), center = TRUE) +Warning message: +In glmboost.formula(fm, data = ovarian, family = CoxPH(), weights = wts, : + model with centered covariates does not contain intercept +> A1 <- survfit(fitw, censor = FALSE) +> A6 <- survFit(fit6) +> print(max(abs(A1$surv - A6$surv)) < 0.001) +[1] TRUE +> A1 <- survfit(fitw, newdata = newdata, censor = FALSE) +> A6 <- survFit(fit6, newdata = newdata) +> print(max(abs(A1$surv - A6$surv)) < 0.001) +[1] TRUE +> > ### Poisson models > > df <- data.frame(x1 = runif(100), x2 = runif(100)) From 11c221e36214ba44f80c51a0bb3cb145d82840c3 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 23 Jul 2026 11:39:19 +0000 Subject: [PATCH 2/4] Initial plan From adc11a92f82e2d2dd12b15edbd5df9614cb628ad Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 23 Jul 2026 11:41:02 +0000 Subject: [PATCH 3/4] docs: clarify survFit newdata=NULL centering behavior --- man/survFit.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/survFit.Rd b/man/survFit.Rd index c85c7eef..f3af1218 100644 --- a/man/survFit.Rd +++ b/man/survFit.Rd @@ -23,7 +23,7 @@ \details{ If \code{newdata = NULL}, the survivor function of the Cox proportional - hazards model is computed for the mean of the covariates used in the + hazards model is computed for a centered linear predictor from the \code{\link{blackboost}}, \code{\link{gamboost}}, or \code{\link{glmboost}} call. The Breslow estimator is used for computing the baseline survivor function. If \code{newdata} is a data frame, the \code{\link{predict}} method From cd37cb4d9bea4f680ed8ec9dbc6bc920598ccf75 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 23 Jul 2026 11:41:38 +0000 Subject: [PATCH 4/4] docs: refine centered linear predictor wording --- man/survFit.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/survFit.Rd b/man/survFit.Rd index f3af1218..a0d67c36 100644 --- a/man/survFit.Rd +++ b/man/survFit.Rd @@ -23,7 +23,7 @@ \details{ If \code{newdata = NULL}, the survivor function of the Cox proportional - hazards model is computed for a centered linear predictor from the + hazards model is computed for a centered linear predictor based on the \code{\link{blackboost}}, \code{\link{gamboost}}, or \code{\link{glmboost}} call. The Breslow estimator is used for computing the baseline survivor function. If \code{newdata} is a data frame, the \code{\link{predict}} method