廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R による打ち切りデータのヒストグラム

こちらを参照

以下に本文中のコードを
#あまり信用しないでください

library(survival)
cdh <- function(DT,DF, digits =0, table=FALSE){
  #Scott's choice
  bw <- zapsmall((3.5*sd(DT))/length(DT)^(1/3), digits=digits) 
  if(bw==0){bw=0.5} #(>_<)
  nb <- ceiling((max(DT)-min(DT))/bw)
  bins <- (1:nb)*bw
  S1 <- survfit(Surv(DT,DF)~1)
  o1 <- order(DT)
  DT <- DT[o1]
  DF <- DF[o1]
  tab1 <- data.frame(DT = S1$time ,DF = DF[!duplicated(DT)], surv =S1$surv)
  TP <- rev(diff(c(sapply(length(bins):1, function(i) min(tab1[tab1$DT <= bins[i],3])),1)))
  barplot(TP/bw, space=0, col="#c8c8cb")
  axis(side=1,labels=c(0, bins),at=0:nb)
  if(table){lis1 <- list(nb, bw, tab1, data.frame(bins,TP))
    names(lis1) <- c("number","width","table1","table2") 
    return(lis1)
  }
}


例)

#png("tru1.png")
truehist(gehan$time)
#dev.off()
library(survival)
library(MASS)
#png("cdh1.png")
cdh(gehan$time, gehan$cens)
#dev.off()

truehist
f:id:abrahamcow:20121125233825p:plain

cdh
f:id:abrahamcow:20121125233835p:plain


参考文献:Huzurbazar (2003)

Huzurbazar (2003) 打ち切りデータのヒストグラム(R による拡張版) - 廿TT も参照されたし