ワイブル確率紙にデータセットを、プロットして最小二乗法で直線を当てはめる関数。 複数のデータセットを同時にプロットすることができる。

weibull.probability.plot <- function(data=data.df, 
                                     log.buff=c(0.2, 0.2), 
                                     points.pch=NULL, 
                                     points.color=NULL,
                                     lines.color=NULL, 
                                     color="grey",
                                     trimming=NULL)
{

  x <- data[,1]
  x.condition <- data[,2]
  log.x.min <- min(log10(x))
  log.x.max <- max(log10(x))
  data.condition <- unique(x.condition)

  weibull.prob.transform <- function(p) {   # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更

    return(log10(log10(1/(1-p))))
  }
  log.axis <- function(z)                         # 対数軸を描く関数
  {
    z <- floor(log10(z))                    # 対数にしたときの整数部
    log.min <- min(z)                       # 最小値
    z2 <- 1:10*10^log.min                   # 値の範囲をカバーするように
    n <- max(z)-log.min                     # 10 倍しながら順次，右の位置に目盛りを描く
    z2 <- rep(z2, n+1)*10^rep(0:n, each=10) # 対数目盛り位置の数値
    log.z2 <- log10(z2)                     # 目盛りを描く位置
    axis(1, at=log.z2, labels=z2)           # log.z2 の位置に，z2 という数値を描く
    abline(v=log.z2, col=color)             # 垂直格子線を描く
  }

  max.n <- 0
  for( i in data.condition ) {
    log.x <- log10(sort(x[x.condition==i]))
    n <- length(log.x)
    max.n <- max(max.n, n)
  }

  log.x <- log10(sort(x))
  log.x <- append(min(log.x)-log.buff[1], log.x)
  log.x <- append(log.x, max(log.x)+log.buff[2])

#  print(max(log.x))
#  print(min(log.x))
#  print(max.n)

  probs <- c(10^(-10:0), 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999)/100
  probs <- probs[probs > 0.1/max.n]
  probs.transformed <- weibull.prob.transform(probs)

  plot(c(log.x[1], log.x[length(log.x)]), 
       c(probs.transformed[1], probs.transformed[length(probs.transformed)]), 
       type="n", xaxt="n", yaxt="n",
       xlab="Observed Value", 
       ylab="Cumulative Percent",
       main="Weibull Probability Paper")

  abline(h=probs.transformed, col="grey")         # 水平の格子線
  log.axis(x)                                 # 横軸を描く
  axis(2, at=probs.transformed, labels=probs)     # 縦軸を描く

  for( i in data.condition) {
    log.x <- log10(sort(x[x.condition==i]))
    n <- length(log.x)
    y <- 1:n/(n+1)
    k <- match(i, data.condition)
    if(!is.null(points.pch)&&!is.null(lines.color)) {
      points(log.x, weibull.prob.transform(y), pch=points.pch[k], col=points.color[k])
    } else {
      points(log.x, weibull.prob.transform(y), pch=20)
    }
    if( (!is.null(trimming)) & (length(trimming)==2) ) {
      i.min <- round(n*trimming[1])
      i.max <- round(n*trimming[2])
    } else {
      i.min <- 1
      i.max <- n
    }
    y.to.fit <- y[i.min:i.max]
    log.x.to.fit <- log.x[i.min:i.max]
    lm.fit <- lm(weibull.prob.transform(y.to.fit)~log.x.to.fit)
    print(summary(lm.fit)$coefficients)
    print(summary(lm.fit)$r.squared)
    if(!is.null(lines.color)) {
      abline(lm.fit, col=lines.color[k])
    } else {
      abline(lm.fit)
    }
  }
}

使い方は

wb.data <- data.frame(time=c(1,2,3,4,5,6,7,8,9,10), condition=c("A","A","B","A","B","B","A","B","A","B"))
weibull.probability.plot(data=wb.data,
                                     log.buff=c(0.2, 0.2),
                                     points.pch=NULL,
                                     points.color=NULL,
                                     lines.color=NULL,
                                     color="grey",
                                     trimming=NULL)

など。

ワイブル確率紙に複数のデータセットを、同時にプロットして最小二乗法で直線を当てはめる関数。 複数のデータセットを同時にプロットすることができる上、打ち切りデータにも対応している。

weibull.probability.plot.2 <- function(data=data.df, 
                                       log.buff=c(0.2, 0.2), 
                                       points.pch=NULL, 
                                       points.color=NULL,
                                       lines.color=NULL, 
                                       color="grey",
                                       trimming=NULL,
                                       debugging=FALSE,
                                       xlim=NULL)
{

  x <- data[,1]
  x.condition <- data[,2]
  x.status <- data[,3]
  # x <- data$time
  # x.condition <- data$condition
  # x.status <- data$status
  log.x.min <- min(log10(x))
  log.x.max <- max(log10(x))
  data.condition <- unique(sort(x.condition))
  
  weibull.prob.transform <- function(p) {   # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更
    return(log10(log10(1/(1-p))))
    # return(log10(p))
    # return(p)
  }
  weibull.prob.transform2 <- function(p) {   # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更
    # return(log10(log10(1/(1-p))))
    return(log10(p))
  }
  
  log.axis <- function(z)                         # 対数軸を描く関数
  {
    z <- floor(log10(z))                    # 対数にしたときの整数部
    log.min <- min(z)                       # 最小値
    z2 <- 1:10*10^log.min                   # 値の範囲をカバーするように
    n <- max(z)-log.min                     # 10 倍しながら順次，右の位置に目盛りを描く
    z2 <- rep(z2, n+1)*10^rep(0:n, each=10) # 対数目盛り位置の数値
    log.z2 <- log10(z2)                     # 目盛りを描く位置
    axis(1, at=log.z2, labels=z2)           # log.z2 の位置に，z2 という数値を描く
    abline(v=log.z2, col=color)             # 垂直格子線を描く
  }
  
  max.n <- 0
  for( i in data.condition ) {
    log.x <- log10(sort(x[x.condition==i]))
    n <- length(log.x)
    max.n <- max(max.n, n)
  }
  
  log.x <- log10(sort(x))
  log.x <- append(min(log.x)-log.buff[1], log.x)
  log.x <- append(log.x, max(log.x)+log.buff[2])
  
  #  print(max(log.x))
  #  print(min(log.x))
  #  print(max.n)
  
  probs <- c(10^(-10:0), 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999)
  # probs <- c(10^(-10:0),2:10,20,30,40,50,100,200,500)
  probs <- probs[probs > 10/max.n]
  probs.transformed <- weibull.prob.transform(probs/100)
  # probs.transformed <- weibull.prob.transform2(probs/100)
  
  if(is.null(xlim)){
    plot(c(log.x[1], log.x[length(log.x)]), 
         c(probs.transformed[1], probs.transformed[length(probs.transformed)]), 
         type="n",
         xaxt="n",
         yaxt="n",
         xlab="Observed Value", 
         ylab="F(t) %",
         main="Weibull Probability Paper")
    # main="Cumulative Hazard Paper")
  }else{
    plot(c(log.x[1], log.x[length(log.x)]), 
         c(probs.transformed[1], probs.transformed[length(probs.transformed)]), 
         type="n",
         xaxt="n",
         yaxt="n",
         xlab="Observed Value", 
         ylab="F(t) %",
         main="Weibull Probability Paper",
         xlim=xlim)
    # main="Cumulative Hazard Paper")
  }
  
  abline(h=probs.transformed, col="grey")         # 水平の格子線
  log.axis(x)                                     # 横軸を描く
  axis(2, at=probs.transformed, labels=probs)     # 縦軸を描く
  
  
  #  i <- data.condition[1]
  for( i in data.condition) {
    if(debugging==TRUE) print(i);
    log.x <- log10(x[x.condition==i]) 
    st <- x.status[x.condition==i]
    x.sort <- sort.list(log.x)
    log.x <- log.x[x.sort]
    st <- st[x.sort]
    n <- length(log.x)
    
    hazard <- (st/c(length(st):1))
    # Hazard <-  cumsum(hazard) * 100
    Hazard <- cumsum(hazard) 
    # y <- { 1-exp(-Hazard/100) } 
    y <- 1-exp(-Hazard) 
    # y <- Hazard
    
    y <- y[st==1] # statusが1を抽出
    log.x <- log.x[st==1] # statusが1を抽出
    
    k <- match(i, data.condition)
    if(!is.null(points.pch)&&!is.null(lines.color)) {
      points(log.x, weibull.prob.transform(y), pch=points.pch[k], col=points.color[k])
    } else {
      points(log.x, weibull.prob.transform(y), pch=20,col=weibull.color(i))
      # points(log.x, weibull.prob.transform2(y), pch=20)
    }
    if( (!is.null(trimming)) & (length(trimming)==2) ) {
      i.min <- round(n*trimming[1])
      i.max <- round(n*trimming[2])
    } else {
      i.min <- 1
      i.max <- length(log.x)
    }
    
    y.to.fit <- y[i.min:i.max]
    log.x.to.fit <- log.x[i.min:i.max]
    lm.fit <- lm(weibull.prob.transform(y.to.fit)~log.x.to.fit) # for plot
    # lm.fit <- lm(weibull.prob.transform2(y.to.fit)~log.x.to.fit) # for plot
    # lm.fit <- lm(y.to.fit~log.x.to.fit) # for plot
    print(lm.fit)
    print(summary(lm.fit)$coefficients)
    print(summary(lm.fit)$r.squared)

    if(!is.null(lines.color)) {
      abline(lm.fit, col=lines.color[k])
    } else {
      abline(lm.fit, col=weibull.color(i))
    }
    
    # for m,eta
    y <- Hazard
    y <- y[st==1] # statusが1を抽出
    y.to.fit <- y[i.min:i.max]
    lm.fit <- lm(weibull.prob.transform2(y.to.fit)~log.x.to.fit) 
    m <- summary(lm.fit)$coefficients[2]
    # eta <- 10^((log10(100)-summary(lm.fit)$coefficients[1])/m)
    eta<- 10^(-summary(lm.fit)$coefficients[1]/m)
    print(paste("m",m,"eta",sprintf("%7.3e",eta)))
    mu <- eta * gamma(1 + 1/m)
    sigma <- eta * ( gamma(1 + 2/m) - (gamma(1 + 1/m))^2  )^(1/2)
    print(paste("mu",sprintf("%7.3e", mu),"sigma",sprintf("%7.3e",sigma)))
    cat("B10:"); cat( sprintf("%7.3e", eta*(-log(1-0.1))^(1/m) ));
    cat(" B1:"); cat( sprintf("%7.3e", eta*(-log(1-0.01))^(1/m) ));
    cat(" B0.1:"); cat( sprintf("%7.3e", eta*(-log(1-0.001))^(1/m) ));
    cat(" B0.01:"); cat( sprintf("%7.3e", eta*(-log(1-0.0001))^(1/m) )); cat("\n")
  }
}

使い方は次の通り。

wb.data.2 <- data.frame(time=c(1,3,4,4,5,5,5,2,2,4),
condition=c(1:10)/c(1:10), status=c(1,1,1,1,1,1,1,1,1,1))
weibull.probability.plot.2(data=wb.data.2,
                                     log.buff=c(0.2, 0.2),
                                     points.pch=NULL,
                                     points.color=NULL,
                                     lines.color=NULL,
                                     color="grey",
                                     trimming=NULL)

weibull.color()

weibull.color <- function(i=i)
{
  
  if(i == 1){
    return("blue")
  }
  if(i == 2){
    return("red")
  }
  if(i == 3){
    return("black")
  }
  if(i == 4){
    return("darkturquoise")
  }
  if(i == 5){
    return("darkorange")
  }
  if(i == 6){
    return("dimgrey")
  }
}