Используйте объект Surv в ggplot или решетке
Кто-либо знает, как использовать в своих интересах ggplot или решетку в выполнении анализа выживания? Было бы хорошо сделать решетку или подобные фасету графики выживания.
Таким образом в конце я играл вокруг и вид найденных решение для графика Kaplan-Meier. Я приношу извинения за грязный код во взятии элементов списка в кадр данных, но я не мог выяснить иначе.
Примечание: Это только работает с двумя уровнями слоев. Если кто-либо знает, как я могу использовать x<-length(stratum) чтобы сделать это, сообщите мне (в Stata, который я мог добавить к макронеуверенному, как это работает в R).
ggkm<-function(time,event,stratum) {
m2s<-Surv(time,as.numeric(event))
fit <- survfit(m2s ~ stratum)
f$time <- fit$time
f$surv <- fit$surv
f$strata <- c(rep(names(fit$strata[1]),fit$strata[1]),
rep(names(fit$strata[2]),fit$strata[2]))
f$upper <- fit$upper
f$lower <- fit$lower
r <- ggplot (f, aes(x=time, y=surv, fill=strata, group=strata))
+geom_line()+geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.3)
return(r)
}
2 ответа
Я использовал следующий код в решетке . Первая функция рисует KM-кривые для одной группы и обычно используется как функция panel.group , тогда как вторая добавляет p-значение теста логарифмического ранга для всей панели:
km.panel <- function(x,y,type,mark.time=T,...){
na.part <- is.na(x)|is.na(y)
x <- x[!na.part]
y <- y[!na.part]
if (length(x)==0) return()
fit <- survfit(Surv(x,y)~1)
if (mark.time){
cens <- which(fit$time %in% x[y==0])
panel.xyplot(fit$time[cens], fit$surv[cens], type="p",...)
}
panel.xyplot(c(0,fit$time), c(1,fit$surv),type="s",...)
}
logrank.panel <- function(x,y,subscripts,groups,...){
lr <- survdiff(Surv(x,y)~groups[subscripts])
otmp <- lr$obs
etmp <- lr$exp
df <- (sum(1 * (etmp > 0))) - 1
p <- 1 - pchisq(lr$chisq, df)
p.text <- paste("p=", signif(p, 2))
grid.text(p.text, 0.95, 0.05, just=c("right","bottom"))
panel.superpose(x=x,y=y,subscripts=subscripts,groups=groups,...)
}
Цензура индикатор должен быть 0-1, чтобы этот код работал. Использование будет следующим:
library(survival)
library(lattice)
library(grid)
data(colon) #built-in example data set
xyplot(status~time, data=colon, groups=rx, panel.groups=km.panel, panel=logrank.panel)
Если вы просто используете 'panel = panel.superpose', вы не получите p-значение.
Я начал почти в точности следовать подходу, который вы используете в своем обновленном ответе. Но в survfit раздражает то, что он отмечает только изменения, а не каждый тик - например, он даст вам 0 - 100%, 3 - 88% вместо 0 - 100%, 1 - 100%, 2 - 100%, 3 - 88%. Если вы передадите это в ggplot, ваши линии будут иметь наклон от 0 до 3, а не оставаться плоскими и падать прямо вниз в точке 3. Это может быть хорошо в зависимости от вашего приложения и предположений, но это не классический график KM. Вот как я справился с разным количеством страт:
groupvec <- c()
for(i in seq_along(x$strata)){
groupvec <- append(groupvec, rep(x = names(x$strata[i]), times = x$strata[i]))
}
f$strata <- groupvec
Если уж на то пошло, то вот как я это сделал - но это тоже не совсем KM-график, потому что я не вычисляю оценку KM как таковую (хотя у меня нет цензуры, так что это эквивалентно... я полагаю).
survcurv <- function(surv.time, group = NA) {
#Must be able to coerce surv.time and group to vectors
if(!is.vector(as.vector(surv.time)) | !is.vector(as.vector(group))) {stop("surv.time and group must be coercible to vectors.")}
#Make sure that the surv.time is numeric
if(!is.numeric(surv.time)) {stop("Survival times must be numeric.")}
#Group can be just about anything, but must be the same length as surv.time
if(length(surv.time) != length(group)) {stop("The vectors passed to the surv.time and group arguments must be of equal length.")}
#What is the maximum number of ticks recorded?
max.time <- max(surv.time)
#What is the number of groups in the data?
n.groups <- length(unique(group))
#Use the number of ticks (plus one for t = 0) times the number of groups to
#create an empty skeleton of the results.
curves <- data.frame(tick = rep(0:max.time, n.groups), group = NA, surv.prop = NA)
#Add the group names - R will reuse the vector so that equal numbers of rows
#are labeled with each group.
curves$group <- unique(group)
#For each row, calculate the number of survivors in group[i] at tick[i]
for(i in seq_len(nrow(curves))){
curves$surv.prop[i] <- sum(surv.time[group %in% curves$group[i]] > curves$tick[i]) /
length(surv.time[group %in% curves$group[i]])
}
#Return the results, ordered by group and tick - easier for humans to read.
return(curves[order(curves$group, curves$tick), ])
}