廿TT

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

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

生存関数のプロットは便利だけど密度関数のプロットと比べると分布の形状を把握しにくい。

そこで打ち切りデータのヒストグラムというのが提案されている(Huzurbazar, A. V. (2005). A Censored Data Histogram. Communications in Statistics - Simulation and Computation, 34 : pp. 113-120. http://www.tandfonline.com/doi/abs/10.1081/SAC-200047089

なんらかの方法で生存率が求まれば適当な幅 bw で生存率の差分を取り、ビンの幅で割ることで密度関数のノンパラメトリックな推定量を得ることができる。

それを関数化しました(make_df4hist · GitHub)。

シミュレーションで右打ち切りのデータを作り、密度関数のカーブと合わせてプロットすると、だいたい一致していることがわかる。

library(survival)
library(tidyverse)
source("https://gist.githubusercontent.com/abikoushi/b9c56028929dc0720e27e36d784a027e/raw/069dfd04ea76885b1186c93eefb75f928cdf6270/make_df4hist")
set.seed(1); x <- rweibull(10000,2,2)
d <- x<2
x2 <- ifelse(d,x,2)

sf_test <- survfit(Surv(x2,d)~1)
bw <- 0.2
dfhisttest <- make_df4hist(sf_test,bw)

p_test <- ggplot(dfhisttest,aes(x=midtime,y=density))+
  geom_col(fill="white",colour="black",width = bw)+
  stat_function(fun = dweibull, args = list(shape=2,scale=2))

print(p_test)

f:id:abrahamcow:20180124103039p:plain

パラメトリックモデルの当てはまりを見るのなんかにも使えるかもしれない。

lung2 <- lung %>% 
  filter(ph.ecog!=3) %>% 
  mutate(sex=sex-1)

sf <-survfit(Surv(time,status)~sex+ph.ecog,data=lung2)
sr <-survreg(Surv(time,status)~sex+ph.ecog,data=lung2)

bw <- 120
df4hist <- make_df4hist(sf,bw)

dfparam <- df4hist %>% 
  mutate(sex=as.integer(sex),ph.ecog=as.integer(ph.ecog)) %>% 
  group_by(strata) %>% 
  summarise(m=1/sr$scale,eta=first(exp(sr$coefficients[1]+sr$coefficients[2]*sex+sr$coefficients[3]*ph.ecog)))

xv <- seq(0, max(sf$time), len=100)
dfdens <- do.call(rbind,lapply(1:nrow(dfparam),function(i){
  with(dfparam[i,],data.frame(strata, x=xv, y=dweibull(xv,shape = m, scale = eta),stringsAsFactors = FALSE))}))

p2<-ggplot(df4hist)+
  geom_col(aes(x=midtime,y=density),width = bw,fill="white",colour="black")+
  geom_line(data = dfdens,aes(x=x,y=y),colour="royalblue",size=1)+
  facet_wrap(~strata)+
  xlab("elapsed time")
print(p2)

f:id:abrahamcow:20180124103020p:plain

abrahamcow.hatenablog.com

abrahamcow.hatenablog.com