Weibull – non-censor and censor

non-censored censored

library(patchwork)
library(ggplot2)

n1    <- 40
shape <- 2
scale <- 5

graphs <- 10

TotXaxis <- c()
fx       <- c()
grpnr    <- c()

for (j in 1:graphs) {                  # Creates graphs.
  
  raknare <- paste0(j,"(10)")
  
  time  <- rweibull(n1, shape, scale)
  
#  timeny <- ifelse(time > limit, NA, time)
#  timeLine <- ifelse(time > limit, limit, time)
  
  yaxel <- c(1:n1)
  
 # nyYaxel <- yaxel[timeLine == limit]
 # grpnr   <- nyYaxel
  
  # x   <- c(rep(limit, length(nyYaxel)), rep(16, length(nyYaxel)))
  # y   <- c(nyYaxel, nyYaxel)               
  # grp <- c(nyYaxel, nyYaxel)
  # restLines <- data.frame(x, y, grp)
  # 
  df <- data.frame(time, yaxel) 
  
  diagram1 <- ggplot(df, aes(x=time, y=yaxel)) + geom_point() + geom_segment(x=0, y=yaxel, xend=time, yend=yaxel) + xlim(0, 16) + ylim(0, 42)
  diagram1 <- diagram1 + annotate("text", x = 15, y = Inf, label=raknare, vjust = 1.5, size = 4, colour = "black" ) + xlab("time") + ylab("item number")
  diagram1 <- diagram1 + annotate("text", x =  1, y = Inf, label="Non-censored times", hjust=0, vjust =  1.5, size = 4, colour = "red" ) 
  
  suppressWarnings(print(diagram1))

# -------------------------------- Graph 1 above, graph 2 below this line
  
n     <- 150
shape <- 2
scale <- 5

graphs  <- 10
maxHist <- 0.23

Xend  <- 16
Xaxis <- seq(0, Xend, by = 0.1)      # For plotting lines

time       <- rweibull(n, shape, scale)
  fitWeibull <- survreg(Surv(time) ~ 1, dist="weibull")
  
  estScale <- exp(fitWeibull$icoef[1])
  estShape <- 1/exp(fitWeibull$icoef[2])
  
  nyttfx   <- c(dweibull(Xaxis, estShape, estScale))
  
  TotXaxis <- c(TotXaxis, Xaxis)
  fx       <- c(fx, nyttfx)
  grpnr    <- c(grpnr, rep(j, length(Xaxis)))
  
  greylines <- data.frame(TotXaxis, fx, grpnr)
  
  raknare <- paste0(j,"(10)")
  shapetext <- paste0("Est. shape = ", round(estShape, 1))
  scaletext <- paste0("Est. scale = ", round(estScale, 1))
  
  data   <- data.frame(time)
  rubrik <- "Weibull distribution (non-censured data) (2, 5)"
  
  truefunktion <- function(x){dweibull(x, shape, scale)}
  
  histW <- ggplot(data, aes(time)) + geom_histogram(binwidth=1, color='black', alpha = 0.5, boundary = 0, fill = "yellow",
                                                    aes(y = after_stat(density)), position = 'identity')
  
  histW <- histW + geom_line(data=greylines, aes(x=TotXaxis, y=fx, group=grpnr), color="grey")
  
  histW <- histW + stat_function(fun=truefunktion, geom="line", color='red', linewidth=1.0)
  histW <- histW + theme(legend.position = "none") + xlim(0, 16) + ylim(0, maxHist)
  histW <- histW + theme(axis.text.x  = element_text(size = 12), axis.text.y = element_blank(),  
                         axis.title.x = element_blank(), axis.title.y = element_blank())
  histW <- histW + annotate("text", x = 5, y = Inf, label=rubrik, vjust = 1.5, size = 4.5, colour = "blue" )
  histW <- histW + annotate("text", x = 10, y = Inf, label=shapetext, hjust=0, vjust =  7, size = 4, colour = "red" ) 
  histW <- histW + annotate("text", x = 10, y = Inf, label=scaletext, hjust=0, vjust = 10, size = 4, colour = "red" ) 
  histW <- histW + annotate("text", x = 15, y = Inf, label=raknare, vjust = 1.5, size = 4, colour = "black" ) 
  
  print(histW)
  
  tvanonCensgrafer <- diagram1 / histW    # 
  
  print(tvanonCensgrafer)
  
  filnamntext = paste0("tvanonCensgrafer", j)
  filnamn <- paste("~/Desktop/Rgrafer/", filnamntext,".jpeg", sep="")
  ggsave(filnamn, width = 7.5, height = 6.5)                                     # Units i tum
  
}

library(patchwork)
library(ggplot2)

n1    <- 40
  shape <- 2
  scale <- 5
  limit <- 6
  
  graphs <- 10
  
  TotXaxis <- c()
  fx       <- c()
  grpnr    <- c()
  grpnr2   <- c()
  
  areaRight <- pweibull(limit, shape, scale, lower.tail = FALSE)
  
  for (j in 1:graphs) {                  # Creates graphs.
    
  raknare <- paste0(j,"(10)")
  
  time  <- rweibull(n1, shape, scale)
  
  timeny <- ifelse(time > limit, NA, time)
  
  timeLine <- ifelse(time > limit, limit, time)
  
  yaxel <- c(1:n1)
  
  nyYaxel <- yaxel[timeLine == limit]
  grpnr   <- nyYaxel
  
  x   <- c(rep(limit, length(nyYaxel)), rep(16, length(nyYaxel)))
  y   <- c(nyYaxel, nyYaxel)               
  grp <- c(nyYaxel, nyYaxel)
  restLines <- data.frame(x, y, grp)
  
  df <- data.frame(timeny, yaxel) 
  
  diagram <- ggplot(df, aes(x=timeny, y=yaxel)) + geom_point() + geom_segment(x=0, y=yaxel, xend=time, yend=yaxel) + xlim(0, 16) + geom_vline(xintercept = limit, color="red")
  diagram <- ggplot(df, aes(x=timeny, y=yaxel)) + geom_point() + geom_segment(x=0, y=yaxel, xend=timeLine, yend=yaxel) + xlim(0, 16) +
                       geom_line(data=restLines, aes(x=x, y=y, group=grp), color="grey") + ylim(0, 42) + xlab("time") + ylab("item number")
  diagram <- diagram + annotate("text", x = 15, y = Inf, label=raknare, vjust = 1.5, size = 4, colour = "black" ) 
  diagram <- diagram + annotate("text", x =  7, y = Inf, label="Censored times (right)", hjust=0, vjust =  1.5, size = 4, colour = "red" ) 
  diagram <- diagram + annotate("text", x =  0, y = Inf, label="Non-censored times (left)", hjust=0, vjust =  1.5, size = 4, colour = "red" ) 
  diagram <- diagram + geom_vline(xintercept = limit, color="red")
  
    suppressWarnings(print(diagram))

# -------------------------------- Graph 1 above, graph 2 below this line
    
  n2    <- 200
  shape <- 2
  scale <- 5
  limit <- 6
  
  Xend  <- 16
  Xaxis <- seq(0, Xend, by = 0.1)      # For plotting lines
    
  time   <- rweibull(n2, shape, scale)
  status <- ifelse(time > limit, 0, 1)       # 1 = 'no fail', 0 = 'fail at time'
  time   <- ifelse(time > limit, limit, time)
  
 timeStatus <- data.frame(time, status)
 
 fitWeibull <- survreg(Surv(time, status) ~ 1, dist="weibull")
 
 estScale <- exp(fitWeibull$icoef[1])
 estShape <- 1/exp(fitWeibull$icoef[2])
 
 nyttfx   <- c(dweibull(Xaxis, estShape, estScale))
 
 TotXaxis <- c(TotXaxis, Xaxis)
 fx       <- c(fx, nyttfx)
 grpnr2   <- c(grpnr2, rep(j, length(Xaxis)))
 
 greylines <- data.frame(TotXaxis, fx, grpnr2)
 
 raknare <- paste0(j,"(10)")
 shapetext <- paste0("Est. shape = ", round(estShape, 1))
 scaletext <- paste0("Est. scale = ", round(estScale, 1))
 
 data   <- data.frame(time)
 rubrik <- "Weibull distribution (censured data) (2, 5)"
 
 truefunktion <- function(x){dweibull(x, shape, scale)}
 
 histCens <- ggplot(data, aes(time)) + geom_histogram(binwidth=1, color='black', alpha = 0.5, boundary = 0, fill = "yellow",
                                                   aes(y = after_stat(density)), position = 'identity')
 
 histCens <- histCens + geom_line(data=greylines, aes(x=TotXaxis, y=fx, group=grpnr2), color="grey")
 
 histCens <- histCens + stat_function(fun=truefunktion, geom="line", color='red', linewidth=1.0)
 histCens <- histCens + theme(legend.position = "none") + xlim(0, 16) + ylim(0, 0.45) + xlab("time")
 histCens <- histCens + theme(axis.text.x  = element_text(size = 12), axis.text.y = element_blank(),  
                         axis.title.y = element_blank())
 histCens <- histCens + annotate("text", x =  5, y = Inf, label=rubrik, vjust = 1.5, size = 4.5, colour = "blue" )
 histCens <- histCens + annotate("text", x = 10, y = Inf, label=shapetext, hjust=0, vjust =  7, size = 4, colour = "red" ) 
 histCens <- histCens + annotate("text", x = 10, y = Inf, label=scaletext, hjust=0, vjust = 10, size = 4, colour = "red" ) 
 histCens <- histCens + annotate("text", x = 15, y = Inf, label=raknare, vjust = 1.5, size = 4, colour = "black" ) 
 histCens <- histCens + geom_segment(x = limit, y = 0, xend = limit, yend = 0.40, color = "red")
 
 suppressWarnings(print(histCens))
 
 tvaCensgrafer <- diagram / histCens    # 
 
 filnamntext = paste0("tvaCensgrafer", j)
 filnamn <- paste("~/Desktop/Rgrafer/", filnamntext,".jpeg", sep="")
 ggsave(filnamn, width = 7.5, height = 6.5)                                     # Units i tum
 
}

Show/hide comments

Show/hide R-code

Weibull – analysis of time data