242 lines
7.6 KiB
R
242 lines
7.6 KiB
R
|
#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)
|