Kaplan-Meier平滑曲线起始值大于1的问题及修复难题

问题背景 #

我参考Stack Overflow上的示例实现了R中平滑Kaplan-Meier曲线的绘制，代码如下：

library(tidyverse)
library(survival)
library(purrr)
library(ggsurvfit)

## Data
df <- survfit(Surv(time, status) ~ surg, data = ggsurvfit::df_colon) |> ggsurvfit::tidy_survfit(type = "survival")

df_split <- split(df,df$strata)

df_smoothed <- purrr::reduce(lapply(c("estimate","conf.low", "conf.high"), function(j) {
  do.call(rbind,
          lapply(seq_along(df_split), function(i) {
            nms <- names(df_split)[i]
            y <-
              predict(mgcv::gam(as.formula(paste0(
                j[[1]], " ~ s(time, bs = 'cs')"
              )), data = df_split[[i]]))
            df <- data.frame(df_split[[i]]$time, y, nms)
            names(df) <- c("time", paste0(j[[1]], ".smooth"), "strata")
            df
          }))
}),dplyr::full_join) |> full_join(df)
#> Joining with `by = join_by(time, strata)`
#> Joining with `by = join_by(time, strata)`
#> Joining with `by = join_by(time, strata)`

## Plotting
ggplot(data=df_smoothed) + 
  geom_line(aes(x=time, y=estimate.smooth, color = strata))+
  geom_ribbon(aes(x=time, ymin = conf.low.smooth, ymax = conf.high.smooth, fill = strata), alpha = 0.50)

这段代码可生成4个分层的平滑Kaplan-Meier曲线，但因不能使用阶梯函数（会泄露微观数据），只能采用平滑曲线。不过各分层曲线的起始值略大于1（初始时刻样本存活率应为100%），推测是平滑逻辑导致的问题。

尝试解决遇到的问题 #

我尝试引入cobs包约束曲线起始于1并保持水平，但修复后仅显示一个分层，且不知如何在cobs中融入survfit函数中的权重参数，尝试代码如下：

library(tidyverse)
library(survival)
library(purrr)
library(ggsurvfit)

## Data
df <- survfit(Surv(time, outcome) ~ exposure+strata(sex), data = my_data) |> ggsurvfit::tidy_survfit(type = "survival")

df_split <- split(df,df$strata)

df_smoothed <- purrr::reduce(lapply(c("estimate","conf.low", "conf.high"), function(j) {
  do.call(rbind,
          lapply(seq_along(df_split), function(i) {
            nms <- names(df_split)[i]
            y <-
              predict(mgcv::gam(as.formula(paste0(
                j[[1]], " ~ s(time, bs = 'cs')"
              )), data = df_split[[i]]))
            df <- data.frame(df_split[[i]]$time, y, nms)
            names(df) <- c("time", paste0(j[[1]], ".smooth"), "strata")
            df
          }))
}),dplyr::full_join) |> full_join(df)
#> Joining with `by = join_by(time, strata)`
#> Joining with `by = join_by(time, strata)`
#> Joining with `by = join_by(time, strata)`


library(cobs) ##NEW CODING STARTS FROM HERE

pw <- rbind(c( 1,min(df_smoothed1$time),1), 
          c(-1,max(df_smoothed$time),0)) 
x <- df_smoothed$time
y <- df_smoothed$estimate.smooth
ft <- cobs(x,y, constraint="decrease", nknots=4,pointwise= con2,
        degree = 2)
fit <- predict(ft, x) [, 'fit']

df_smoothed$x <- x
df_smoothed$y <- y
df_smoothed$fit <- fit

## Plotting
ggplot(data=df_smoothed, aes(x,y, color = strata)) + 
  geom_line(aes(y=fit))+
  geom_ribbon(aes(x=time, ymin = conf.low.smooth, ymax = conf.high.smooth, fill = strata), alpha = 0.50)

现寻求解决曲线起始值异常及cobs包适配问题的方案。

内容的提问来源于stack exchange，提问作者Devi Sita