diff --git a/R/survival.R b/R/survival.R index 071eed4..5a50829 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 e1268e9..f2bc351 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 c85c7ee..dc9d80b 100644 --- a/man/survFit.Rd +++ b/man/survFit.Rd @@ -23,13 +23,20 @@ \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 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 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 8336e10..8feea66 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 821dcf1..73b1760 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 d8348d1..a3ef2b7 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 3c9e4cf..dc1cc01 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))