2010-06-01 9 views
5

कोई भी जानता है कि जीवित रहने के विश्लेषण में ggplot या जाली का लाभ कैसे लें? ट्रेली या पहलू-जैसे अस्तित्व ग्राफों को करना अच्छा लगेगा।ggplot या जाली में सर्व ऑब्जेक्ट का उपयोग


तो अंत में मैंने चारों ओर खेला और एक कपलान-मेयर प्लॉट के लिए एक समाधान मिला। सूची तत्वों को डेटाफ्रेम में लेने में गन्दा कोड के लिए मैं क्षमा चाहता हूं, लेकिन मैं एक और तरीका नहीं समझ सका।

नोट: यह केवल स्तर के दो स्तरों के साथ काम करता है। यदि कोई जानता है कि मैं ऐसा करने के लिए x<-length(stratum) का उपयोग कैसे कर सकता हूं तो कृपया मुझे बताएं (स्टेटा में मैं एक मैक्रो-अनिश्चितता में शामिल हो सकता हूं कि यह आर में कैसे काम करता है)।

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) 
} 
+3

रेमन Saccilotto लिखा एक ggplot2 ट्यूटोरियल कि ggplot2 में के.एम. भूखंडों के लिए कार्य शामिल हैं: http://www.ceb-institute.org/bbs/wp-content/uploads/2011/09/handout_ggplot2.pdf – MattBagg

उत्तर

4

मैं निम्नलिखित कोड का उपयोग lattice में कर रहा हूं। पहले समारोह एक समूह के लिए ड्रॉ KM-घटता है और आम तौर panel.group समारोह के रूप में इस्तेमाल किया जाएगा, जबकि दूसरा पूरे पैनल के लिए लॉग-रैंक परीक्षण पी-मूल्य जोड़ता है:

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.superpose' तो आप पी-मूल्य प्राप्त नहीं होगा।

1

मैंने आपके अपडेट किए गए उत्तर में उपयोग किए जाने वाले लगभग बिल्कुल सही दृष्टिकोण का पालन करना शुरू कर दिया। लेकिन जो चीज जीवित रहने के बारे में चिंतित है वह यह है कि यह केवल परिवर्तनों को चिह्नित करता है, न कि प्रत्येक टिक - उदाहरण के लिए, यह आपको 0 - 100%, 3 - 88% 0 - 100%, 1 - 100%, 2 - 100 के बजाय देगा %, 3 - 88%। यदि आप इसे जीजीप्लॉट में खिलाते हैं, तो आपकी लाइनें फ्लैट से बनी रहेंगी और सीधे 3 पर गिरने के बजाय 0 से 3 तक ढल जाएंगी। यह आपके आवेदन और धारणाओं के आधार पर ठीक हो सकता है, लेकिन यह क्लासिक केएम प्लॉट नहीं है। यह मैं कैसे तबके के अलग-अलग संख्याओं को संभाला है:

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 

क्या इसके लायक है के लिए, यह कैसे मैं यह कर समाप्त हो गया है - लेकिन यह नहीं सच में एक किलोमीटर की साजिश, या तो, क्योंकि मैं की गणना नहीं कर रहा हूँ है केएम अनुमान प्रति से बाहर (हालांकि मेरे पास कोई सेंसरिंग नहीं है, इसलिए यह बराबर है ... मुझे विश्वास है)।

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

}