文書の過去の版を表示しています。
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),...) } } 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)) } 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) }