差分

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

この比較画面へのリンク

次のリビジョン		前のリビジョン
r:survival:code [2018/12/09 21:58] – created 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) ) {
行 45: 行 47:
   }   }
 } }
-</code> 
  
 +compare.histgram.and.densities <- function(X,hist=TRUE,plot=TRUE) {
 +  require(MASS)
 +  est.weibull <- fitdistr(X$time,densfun="weibull")
 +  est.lnorm <- fitdistr(X$time,densfun="lognormal")
 +  est.gamma <- fitdistr(X$time,densfun="gamma")
 +  if(plot==TRUE) {
 +    if(hist==TRUE) {
 +      hist(X$time,freq=FALSE,ylim=c(0,0.3),
 +           xlab="Lifetime")
 +      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),
 +         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))
 +}
 +
 +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>