差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:survival [2018/12/09 21:53] – watalur:survival [2018/12/11 19:42] (現在) – watalu
行 1: 行 1:
 ==== Rによる生存時間解析 ==== ==== Rによる生存時間解析 ====
 +
 +[[r:survival:code]]の方に記したコードを、先にRに貼り付けて読み込ませておくこと。いくつかのグラフを描く関数を用意しておいたものである。参考までに、用いているRの関数を個別に詳細に説明したページも用意している。
 +
 +  * [[r:survival:data.frame]]は寿命データの記し方を説明している
 +  * [[r:mass:fitdistr]]は完全データからの確率分布のパラメータの最尤法に基づいて推定する
 +  * [[r:survival:survfit]]は打ち切りを含む寿命データから、生存関数をノンパラメトリックに推定する
 +
 +ただし、実験の内容を実施するためには、上の３つのページは読まずとも問題なく、下の説明に従って進めれば十分である。
  
 === データの書き方 === === データの書き方 ===
行 94: 行 102:
  
 信頼性データを図示する関数draw.reliability.data.diagramを用意した。 信頼性データを図示する関数draw.reliability.data.diagramを用意した。
- 
-<code> 
-draw.reliability.data.diagram <- function(x,sort=TRUE,...) { 
-  if(is.numeric(x)==TRUE) { 
-    n <- length(x) 
-    if(sort==TRUE) { 
-      x.time <- x[sort.list(x,decreasing=TRUE)] 
-      x.status <- rep(1,n) 
-      id <- c(1:n) 
-    } else { 
-      x.time <- x 
-      x.status <- rep(1,n) 
-      id <- c(1:n) 
-    } 
-  } else if (is.data.frame(x)==TRUE) { 
-    n <- dim(x)[1] 
-    if(sort==TRUE) { 
-      x.time <- x[sort.list(x[,1],decreasing=TRUE),1] 
-      x.status <- x[sort.list(x[,1],decreasing=TRUE),2] 
-      id <- row.names(x)[sort.list(x[,1],decreasing=TRUE)] 
-    } else { 
-      x.time <- x[,1] 
-      x.status <- x[,2] 
-      id <- row.names(x) 
-    } 
-  } else if (is.matrix(x)==TRUE) { 
-    n <- dim(x)[1] 
-    if(sort==TRUE) { 
-      x.time <- x[sort.list(x[,1],decreasing=TRUE),1] 
-      x.status <- x[sort.list(x[,1],decreasing=TRUE),2] 
-      id <- rawnames(x)[sort.list(x[,1],decreasing=TRUE)] 
-    } else { 
-      x.time <- x[,1] 
-      x.status <- x[,2] 
-      id <- rawnames(x) 
-    } 
-  } 
-  plot(x.time,c(1:n), 
-       xlim=c(0,max(x.time)), 
-       type="n", 
-       ylab="#", 
-       xlab="Lifetime data",yaxt="n") 
-  points(x.time,c(1:n),pch=c(1,4)[x.status+1],...) 
-  for( i in c(1:n) ) { 
-    lines(c(0,x.time[i]),c(i,i),...) 
-  } 
-} 
-</code> 
- 
 上に与えたデータの図を描くには次の１行を実行する。 上に与えたデータの図を描くには次の１行を実行する。
  
 <code> <code>
-draw.reliability.data.diagram(X)+plot.reliability.data.diagram(X)
 </code> </code>
  
行 157: 行 116:
 </code> </code>
  
-これを実行しておかないと、すぐ次のコマンドの実行時に赤い文字でエラーが表示される。+これを実行しておかないと、すぐ次のコマンドの実行時に青い文字でエラーが表示されるかもしれない。
  
 === 生存曲線と故障率曲線のノンパラメトリック推定 === === 生存曲線と故障率曲線のノンパラメトリック推定 ===
行 298: 行 257:
 === 関数 === === 関数 ===
  
-完全データについて、これら三つの確率分布を推定し、ヒストグラムと比較するための関数を用意した。 +完全データについて、これら三つの確率分布を推定し、ヒストグラムと比較するための関数compare.histgram.and.densitiesを用意した。 
-<code> +
-compare.histgram.and.densities <- function(X,hist=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(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)) +
-} +
-</code>+
 オプションのhistをhist=FALSEとすると、ヒストグラムは描かず、密度関数のみを描くようになる。 オプションのhistをhist=FALSEとすると、ヒストグラムは描かず、密度関数のみを描くようになる。
  
-またカプラン・マイヤー法で推定した生存関数と共に描くための関数も用意した。 +またカプラン・マイヤー法で推定した生存関数と共に描くための関数compare.KM.and.distributionsも用意した。
-<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) +
-} +
-</code> +
 これで、どの分布に最も近いかを検討できる。 これで、どの分布に最も近いかを検討できる。
  
行 398: 行 296:
 </code> </code>
  
-これらをカプラン・マイヤー曲線に重ねて描く関数も用意した。 +これらをカプラン・マイヤー曲線に重ねて描く関数compare.KM.and.distributions.2 も用意した。
- +
-<code> +
-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>+
  
 <code> <code>
 compare.KM.and.distributions.2(X) compare.KM.and.distributions.2(X)
 </code> </code>