survFit: support models fitted with case weights

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)
 

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}.

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.}

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)

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)

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))

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))