差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:survival:code [2018/12/09 22:00] – watalur:survival:code [2020/10/20 16:16] (現在) – jin
行 1: 行 1:
 <code> <code>
-draw.reliability.data.diagram <- function(x,sort=TRUE,...) {+plot.reliability.data.diagram <- function(x,sort=TRUE,...) {
   if(is.numeric(x)==TRUE) {   if(is.numeric(x)==TRUE) {
     n <- length(x)     n <- length(x)
行 39: 行 39:
        type="n",        type="n",
        ylab="#",        ylab="#",
-       xlab="Lifetime data",yaxt="n")+       xlab="Lifetime data", 
 +       yaxt="n") 
 +  axis(2,1:n,labels=c(n:1))
   points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)   points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)
   for( i in c(1:n) ) {   for( i in c(1:n) ) {
行 46: 行 48:
 } }
  
-compare.histgram.and.densities <- function(X,hist=TRUE) {+compare.histgram.and.densities <- function(X,hist=TRUE,plot=TRUE) {
   require(MASS)   require(MASS)
   est.weibull <- fitdistr(X$time,densfun="weibull")   est.weibull <- fitdistr(X$time,densfun="weibull")
   est.lnorm <- fitdistr(X$time,densfun="lognormal")   est.lnorm <- fitdistr(X$time,densfun="lognormal")
   est.gamma <- fitdistr(X$time,densfun="gamma")   est.gamma <- fitdistr(X$time,densfun="gamma")
-  if(hist==TRUE) { +  if(plot==TRUE) { 
-    hist(X$time,freq=FALSE,ylim=c(0,0.3), +    if(hist==TRUE) { 
-         xlab="Lifetime") +      hist(X$time,freq=FALSE,ylim=c(0,0.3), 
-    plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], +           xlab="Lifetime") 
-                                   sdlog=est.lnorm$estimate[2]))},+      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], 
 +                                     sdlog=est.lnorm$estimate[2]))}, 
 +           xlim=c(0,12), 
 +           col="blue", 
 +           add=TRUE) 
 +    } else { 
 +      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], 
 +                                     sdlog=est.lnorm$estimate[2]))}, 
 +           xlim=c(0,12), 
 +           col="blue",xlab="Lifetime",ylab="Density") 
 +    } 
 +    plot(function(x){return(dweibull(x,shape=est.weibull$estimate[1], 
 +                                     scale=est.weibull$estimate[2]))},
          xlim=c(0,12),          xlim=c(0,12),
-         col="blue",+         col="red",
          add=TRUE)          add=TRUE)
-  } else { +    plot(function(x){return(dgamma(x,shape=est.gamma$estimate[1], 
-    plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], +                                   rate=est.gamma$estimate[2]))},
-                                   sdlog=est.lnorm$estimate[2]))},+
          xlim=c(0,12),          xlim=c(0,12),
-         col="blue",xlab="Lifetime",ylab="Density")+         col="green", 
 +         add=TRUE) 
 +    legend("topright", 
 +           col=c("blue","red","green"), 
 +           lty=c(1,1,1), 
 +           legend=c("lognormal","weibull","gamma"))
   }   }
-  plot(function(x){return(dweibull(x,shape=est.weibull$estimate[1], 
-                                   scale=est.weibull$estimate[2]))}, 
-       xlim=c(0,12), 
-       col="red", 
-       add=TRUE) 
-  plot(function(x){return(dgamma(x,shape=est.gamma$estimate[1], 
-                                 rate=est.gamma$estimate[2]))}, 
-       xlim=c(0,12), 
-       col="green", 
-       add=TRUE) 
-  legend("topright", 
-         col=c("blue","red","green"), 
-         lty=c(1,1,1), 
-         legend=c("lognormal","weibull","gamma")) 
   return(list(weibull=est.weibull,lnorm=est.lnorm,gamma=est.gamma))   return(list(weibull=est.weibull,lnorm=est.lnorm,gamma=est.gamma))
 } }
-</code> 
  
 +compare.KM.and.distributions <- function(fits,data) {
 +  require(survival)
 +  km <- survfit(Surv(time,status)~1,data=data)
 +  plot(km,
 +       xlab="Time",
 +       ylab="Survival Function")
 +  plot(function(x) pweibull(x,shape=fits$weibull$estimate[1],scale=fits$weibull$estimate[2],
 +                            lower.tail=FALSE),
 +       xlim=c(0,max(km$time)),add=TRUE,col="green")
 +  plot(function(x) plnorm(x,meanlog=fits$lnorm$estimate[1],sdlog=fits$lnorm$estimate[2],
 +                          lower.tail=FALSE),
 +       xlim=c(0,max(km$time)),add=TRUE,col="blue")
 +  plot(function(x) pgamma(x,shape=fits$gamma$estimate[1],rate=fits$gamma$estimate[2],
 +                          lower.tail=FALSE),
 +       xlim=c(0,max(km$time)),add=TRUE,col="orange")
 +  legend("topright",lty=c(1,2,1,1,1),col=c("black","black","green","blue","orange"),
 +         legend=c("Kaplan-Meier Estimate","95% Confidence Interval",
 +                  paste("Weibull",round(fits$weibull$loglik,5)),
 +                  paste("Lognormal",round(fits$lnorm$loglik,5)),
 +                  paste("Gamma",round(fits$gamma$loglik,5))),
 +         cex=0.5)
 +}
 +
 +compare.KM.and.distributions.2 <- function(data) {
 +  require(survival)
 +  km <- survfit(Surv(time,status)~1,data=data)
 +  fits <- compare.histgram.and.densities(data,hist=FALSE,plot=FALSE)
 +  survreg.weibull <- survreg(Surv(time, status)~1, data=data, dist="weibull")
 +  fits$weibull$estimate[1] <- 1/survreg.weibull$scale
 +  fits$weibull$estimate[2] <- exp(survreg.weibull$coefficients[1])
 +  survreg.lognormal <- survreg(Surv(time, status)~1, data=data, dist="lognormal")
 +  fits$lnorm$estimate[1] <- survreg.lognormal$coefficients[1]
 +  fits$lnorm$estimate[2] <- survreg.lognormal$scale
 +  plot(km,
 +       xlab="Time",
 +       ylab="Survival Function")
 +  plot(function(x) pweibull(x,shape=fits$weibull$estimate[1],scale=fits$weibull$estimate[2],
 +                            lower.tail=FALSE),
 +       xlim=c(0,max(km$time)),add=TRUE,col="green")
 +  plot(function(x) plnorm(x,meanlog=fits$lnorm$estimate[1],sdlog=fits$lnorm$estimate[2],
 +                          lower.tail=FALSE),
 +       xlim=c(0,max(km$time)),add=TRUE,col="blue")
 +  legend("topright",lty=c(1,2,1,1,1),col=c("black","black","green","blue","orange"),
 +         legend=c("Kaplan-Meier Estimate","95% Confidence Interval",
 +                  paste("Weibull",round(fits$weibull$loglik,5)),
 +                  paste("Lognormal",round(fits$lnorm$loglik,5))),
 +         cex=0.5)
 +}
 +</code>