Sử dụng đối tượng Surv trong ggplot hoặc mạng

Question

Bất cứ ai cũng biết cách tận dụng lợi thế của ggplot hoặc mạng trong phân tích sự sống còn? Nó sẽ là tốt đẹp để làm một trellis hoặc facet giống như đồ thị tồn tại.

---

Vì vậy, cuối cùng tôi chơi xung quanh và tìm thấy một giải pháp cho âm mưu Kaplan-Meier. Tôi xin lỗi vì mã lộn xộn trong việc đưa các phần tử danh sách vào một khung dữ liệu, nhưng tôi không thể tìm ra cách khác.

Lưu ý: Nó chỉ hoạt động với hai cấp tầng. Nếu ai biết làm thế nào tôi có thể sử dụng x<-length(stratum) để làm điều này xin vui lòng cho tôi biết (trong Stata tôi có thể thêm vào một vĩ mô không chắc chắn như thế nào này hoạt động trong 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) 
} 

Answer 1

Tôi đã sử dụng mã sau đây trong lattice. Chức năng đầu tiên rút KM-đường cong cho một nhóm và thường sẽ được sử dụng làm panel.group chức năng, trong khi thứ hai cho biết thêm các test thứ hạng log p-giá trị cho toàn bộ bảng điều khiển:

chỉ

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

sự kiểm duyệt phải được 0-1 để mã này hoạt động. Mức sử dụng sẽ nằm dọc theo các dòng sau:

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) 

Nếu bạn chỉ sử dụng 'panel = panel.superpose' thì bạn sẽ không nhận được giá trị p.

Answer 2

Tôi bắt đầu theo dõi gần như chính xác cách tiếp cận bạn sử dụng trong câu trả lời cập nhật của mình. Nhưng điều gây khó chịu cho người sống sót là nó chỉ đánh dấu những thay đổi, không phải mỗi lần đánh dấu - ví dụ, nó sẽ cho bạn 0 - 100%, 3 - 88% thay vì 0 - 100%, 1 - 100%, 2 - 100 %, 3 - 88%. Nếu bạn cho nó vào ggplot, đường của bạn sẽ dốc từ 0 đến 3, chứ không phải bằng phẳng và thả thẳng xuống 3. Điều đó có thể tốt tùy thuộc vào ứng dụng và giả định của bạn, nhưng nó không phải là cốt truyện KM cổ điển. Đây là cách tôi xử lý các con số khác nhau của các tầng lớp nhân dân:

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 

Đối với những gì nó có giá trị, đây là cách tôi đã kết thúc làm việc đó - nhưng điều này là không thực sự là một âm mưu KM, một trong hai, bởi vì tôi không tính toán ra ước tính KM mỗi lần (mặc dù tôi không kiểm duyệt, vì vậy điều này là tương đương ... Tôi tin).

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

}