差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:survival [2018/12/09 21:59] – watalur:survival [2018/12/11 19:42] (現在) – watalu
行 1: 行 1:
 ==== Rによる生存時間解析 ==== ==== Rによる生存時間解析 ====
  
-[[r:survival:code]]の方に記したコードを、先にRに読み込ませておいてください。+[[r:survival:code]]の方に記したコードを、先にRに貼り付けて読み込ませておくこと。いくつかのグラフを描く関数を用意しておいたものである。参考までに、用いているRの関数を個別に詳細に説明したページも用意している。 
 + 
 +  * [[r:survival:data.frame]]は寿命データの記し方を説明している 
 +  * [[r:mass:fitdistr]]は完全データからの確率分布のパラメータの最尤法に基づいて推定する 
 +  * [[r:survival:survfit]]は打ち切りを含む寿命データから、生存関数をノンパラメトリックに推定する 
 + 
 +ただし、実験の内容を実施するためには、上の３つのページは読まずとも問題なく、下の説明に従って進めれば十分である。
  
 === データの書き方 === === データの書き方 ===
行 99: 行 105:
  
 <code> <code>
-draw.reliability.data.diagram(X)+plot.reliability.data.diagram(X)
 </code> </code>
  
行 255: 行 261:
 オプションの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> +
 これで、どの分布に最も近いかを検討できる。 これで、どの分布に最も近いかを検討できる。
  
行 315: 行 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>