#survival # survfit experiments # library(ggsurvfit) ## This is what I wanted (!) (p <- survfit2(Surv(time, status) ~ surg, data = df_colon) |> ggsurvfit(linewidth = 1,) + add_confidence_interval() + add_risktable() + add_quantile(y_value = 0.6, color = "gray50", linewidth = 0.75) + # limit plot to show 8 years and less coord_cartesian(xlim = c(0, 8)) + # update figure labels/titles labs( y = "Percentage Survival", title = "Recurrence by Time From Surgery to Randomization", ) + # reduce padding on edges of figure and format axes scale_y_continuous(label = scales::percent, breaks = seq(0, 1, by = 0.2), expand = c(0.015, 0)) + scale_x_continuous(breaks = 0:10, expand = c(0.02, 0))) head(df_colon) library("survival") library("survminer") variables <- c("sex", "age", "adhere", "extent", "surg") cox_model <- coxph(as.formula(paste("Surv(time, status) ~",paste(variables,collapse="+"))), data = df_colon) summary(cox_model) ggsurvplot(survfit(cox_model), data=df_colon, palette = "#2E9FDF", ggtheme = theme_minimal()) # Checks (ph_check <- survival::cox.zph(cox_model)) survminer::ggcoxzph(ph_check, var=c("sex", "age", "adhere", "surg"), font.main = 10, font.x = 10, font.y = 10) ## Below is the nest best solution ## Try at smoothing the survival curve ## Working code to get satisfying object df <- survfit2(Surv(time, status) ~ surg, data = df_colon) |> 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) 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) + # geom_smooth(aes(x=time, y=estimate, color = strata), method = "gam", formula = y ~ s(x, bs = "cs")) + # reduce padding on edges of figure and format axes scale_y_continuous(label = scales::percent, breaks = seq(0, 1, by = 0.2), expand = c(0.015, 0), limits = c(0,1)) + scale_x_continuous(breaks = 0:10, expand = c(0.02, 0))+ labs( y = "Percentage Survival", title = "Recurrence by Time From Surgery to Randomization", ) + # limit plot to show 8 years and less coord_cartesian(xlim = c(0, 8)) ## Dendrogram df_colon[do.call(c, lapply(seq_len(ncol(df_colon)), function(i) { is.double(df_colon[[i]]) }))][-1] |> scale() |> dist() |> hclust(method="average") |> ggdendro::ggdendrogram() ## Better example?? ## library(tidyverse) library(survival) library(purrr) library(ggsurvfit) library(cobs) ## Data plot.type <- "survival" x <- survfit2(Surv(time, status) ~ surg, data = df_colon) df <- tidy_survfit(x, type = plot.type) %>% dplyr::mutate(survfit = c(list(x), rep_len(list(), dplyr::n() - 1L))) method <- "gam" 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] x = df_split[[i]]$time if (method=="loess"){ y <- predict(loess(as.formula(paste0( j[[1]], " ~ time" )), data = df_split[[i]])) } else if (method=="gam"){ y <- predict(mgcv::gam(as.formula(paste0( j[[1]], " ~ s(time, bs = 'cs')" )), data = df_split[[i]])) } else if (method=="cobs") { if (plot.type=="survival"){ ## This will make the plot start in (0,1) con <- rbind(c( 0,min(x),1)) ## This ensures a monotonic decreasing slope ## for the estimate, not the CIs if (j[[1]]=="estimate"){ direction="decrease" } else {direction="none"} } else if (plot.type=="risk"){ con <- rbind(c( 0,min(x),0)) if (j[[1]]=="estimate"){ direction="increase" } else {direction="none"} } m <- cobs(x,df_split[[i]][[j]], constraint=direction, nknots = 4, pointwise= con, degree = 2,) y <- predict(m, x)[, 'fit'] } df <- data.frame(x, y, nms) names(df) <- c("time", paste0(j[[1]], ".smooth"), "strata") df })) }),dplyr::full_join) |> full_join(df) ## 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) ## Weighted GAM approach ## https://stackoverflow.com/a/66705556/21019325 ## It does not work. Gonna stop here due to lack of time. ## Apparantly dat_orig <- df_split[[1]][,c("time","estimate")] x1=0 y1=1 # set.seed(123) # N = 100 # x <- sort(runif(N) * 4 - 1) # f <- exp(4*x)/(1+exp(4*x)) # y <- f + rnorm(N) * 0.1 # x = c(-1, x) # y = c(-0.1, y) # dat = data.frame(x = x, y= y) x <- do.call(c,c(x1,dat_orig[1])) y <- do.call(c,c(y1,dat_orig[,2])) dat <- data.frame(x=x,y=y) k <- 13 library(mgcv) fit0 <- gam(y ~ s(x, k = k, bs = "cr"),data=dat) # predict from unconstrained GAM fit newdata <- data.frame(x = x) newdata$y_pred_fit0 <- predict(fit0, newdata = newdata) # Show regular spline fit (and save fitted object) # f.ug <- gam(y~s(x,k=k,bs="cr")) # explicitly construct smooth term's design matrix sm <- smoothCon(s(x,k=k,bs="cr"),dat,knots=NULL)[[1]] # find linear constraints sufficient for monotonicity of a cubic regression spline # it assumes "cr" is the basis and its knots are provided as input f.mono <- mono.con(sm$xp,up = FALSE) G <- list( X=sm$X, C=matrix(0,0,0), # [0 x 0] matrix (no equality constraints) sp=fit0$sp, # smoothing parameter estimates (taken from unconstrained model) p=sm$xp, # array of feasible initial parameter estimates y=dat[,2], w= c(1e8, rep(1,nrow(dat_orig))), # weights for data Ain=f.mono$A, # matrix for the inequality constraints bin=f.mono$b, # vector for the inequality constraints S=sm$S, # list of penalty matrices; The first parameter it penalizes is given by off[i]+1 off=0 # Offset values locating the elements of M$S in the correct location within each penalty coefficient matrix. (Zero offset implies starting in first location) ) p <- pcls(G) # fit spline (using smoothing parameter estimates from unconstrained fit) # predict newdata$y_pred_fit2 <- Predict.matrix(sm, data.frame(x = newdata$x)) %*% p # plot ggplot(data=newdata) + geom_line(aes(x=x, y=y_pred_fit2)) # geom_ribbon(aes(x=time, ymin = conf.low.smooth, ymax = conf.high.smooth, fill = strata), alpha = 0.50) # plot(y ~ x, data = dat) lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2) lines(y_pred_fit2 ~ x, data = newdata, col = 4, lwd = 2) abline(v = -1) abline(h = -0.1)