Utilizza oggetto Surv in ggplot o lattice

Question

Qualcuno sa come trarre vantaggio da ggplot o reticolo nel fare analisi di sopravvivenza? Sarebbe bello fare un traliccio o grafici di sopravvivenza simili a faccette.

---

Quindi alla fine ho giocato e ho trovato una soluzione per una trama di Kaplan-Meier. Mi scuso per il codice disordinato nel prendere gli elementi della lista in un dataframe, ma non riuscivo a capire un altro modo.

Nota: Funziona solo con due livelli di strati. Se qualcuno sa come posso usare x<-length(stratum) per farlo, fammi sapere (in Stata potrei aggiungere una macro-incertezza su come funziona in 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

Sono stato con il seguente codice nel lattice. La prima funzione disegna KM-curve per un gruppo e viene in genere utilizzato come funzione panel.group, mentre il secondo aggiunge il p-value log-rank test per l'intero pannello:

indicatore

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

la censura deve essere 0-1 per far funzionare questo codice. L'utilizzo sarebbe lungo le seguenti linee:

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) 

Se si basta usare 'pannello = panel.superpose', allora non sarà possibile ottenere il p-value.

Answer 2

Ho iniziato seguendo quasi esattamente l'approccio utilizzato nella risposta aggiornata. Ma la cosa che è irritante per il sopravvissuto è che segna solo i cambiamenti, non ogni tick - per esempio, ti darà 0 - 100%, 3 - 88% invece di 0 - 100%, 1 - 100%, 2 - 100 %, 3 - 88%. Se lo inserisci in ggplot, le tue linee saranno inclinate da 0 a 3, invece di rimanere piatte e cadere dritte verso il basso a 3. Ciò potrebbe andare bene a seconda dell'applicazione e delle ipotesi, ma non è la classica trama KM. Ecco come ho gestito la numero variabile di strati:

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 

Per quel che vale, è così che ho finito per farlo - ma questo non è davvero un complotto KM, o, perché io non sono il calcolo la stima KM di per sé (anche se non ho censure, quindi questo è equivalente ... credo).

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

}