廿TT

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

rstan で生存関数の推定(離散時間, 2群)

いろいろやり方はあると思うけど, かんたんに書けるので, とりあえずこれでいいかなと思った.

生存関数を S(t) で表し, 離散時間でハザードを考えると,

p_1=1-S(t_1),
p_2=\{S(t_1)-S(t_2)\}/S(t_1),
p_3=\{S(t_2)-S(t_3)\}/S(t_2),...

生存関数を考えると,

S(t_1)=1-p_1,
S(t_2)=(1-p_2)S(t_1),
S(t_3)=(1-p_3)S(t_2),...

各時間区切りごとのイベント数を y, リスクセット(その時間の直前までイベントの生起していない被験者数)を R で表すと, ハザード p に関する尤度は,

 L \propto p^{y}(1-p)^{R-y}

と二項分布で書ける.

さらにハザード p に対してロジットリンクで回帰構造を入れる.

切片項を  \alpha_j とし, 時間ごとに異なるとすると, ベースラインのハザードはだんだん大きくなると考えるのが自然なので, \alpha_{j-1} < \alpha_j とする.

Stan のコードは以下のようになる.

data{
  int<lower=1> N;
  int<lower=1> Tim;
  int<lower=0> Y[N];
  int<lower=0> R[N];
  int<lower=1,upper=Tim> time[N];
  vector<lower=0,upper=1>[N] treat;
}
parameters{
  ordered[Tim] alpha;
  real beta;
}
model{
  for(i in 1:N){
   Y[i] ~ binomial_logit(R[i],alpha[time[i]]+beta*treat[i]); 
  }
}
generated quantities{
  vector[Tim] S1;
  vector[Tim] S2;
  S1[1] = 1-inv_logit(alpha[1]);
  S2[1] = 1-inv_logit(alpha[1]+beta);
  for(i in 2:Tim){
  S1[i] = (1-inv_logit(alpha[i]))*S1[i-1];
  S2[i] = (1-inv_logit(alpha[i]+beta))*S2[i-1];
  }
}

これをベースにして, 回帰係数  \beta が時間とともに変化するようなモデルとか, ランダム効果を入れるとか, \alpha をスムージングするとか, いろいろ考えられると思う(ちゃんと推定できるかはわからない).

R のパッケージ MASS に入っている Gehan データ(白血病寛解時間のデータだったと思う)に対して, 生存関数を推定してみた.

f:id:abrahamcow:20171011045225p:plain

MCMCの収束は大丈夫そうである.

生存関数をプロットしてみる. 帯は95%信用区間.

f:id:abrahamcow:20171011045323p:plain

6-MP(薬)投与群と, 対照群では明確に差が見て取れる.

回帰係数  \beta の95%信用区間に 0 が含まれていないことからも, 薬の効果はあると判断できる.

> summary(survfit,pars="beta")$summary
         mean     se_mean        sd     2.5%      25%      50%     75%
beta 3.184638 0.006650618 0.5730268 2.129856 2.780925 3.167513 3.56786
        97.5%    n_eff     Rhat
beta 4.360173 7423.794 1.000251
> 

以下に R のコードを貼ります.

library(tidyverse)
library(rstan)
data(gehan,package = "MASS")

gehan2 <-gehan %>%
  group_by(time,treat) %>%
  summarise(Y=sum(cens),n=n()) %>% 
  group_by(treat) %>% 
  mutate(R=rev(cumsum(rev(n)))) %>% 
  ungroup() %>% 
  mutate(time2=as.integer(factor(time)))

dat <- list(N = nrow(gehan2),
            Y = gehan2$Y,
            R = gehan2$R,
            Tim = max(gehan2$time2),
            time = gehan2$time2,
            treat=as.integer(gehan2$treat)-1)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
survfit_model <-stan_model("survfit.stan")
survfit <- sampling(survfit_model,dat,iter=5000,
                  control=list(adapt_delta=0.99))

traceplot(survfit,pars=c("alpha","beta"))
all(summary(survfit)$summary[,"Rhat"]<1.1,na.rm = TRUE)

summary(survfit,pars="beta")$summary

gehan3 <-data_frame(time=sort(unique(gehan2$time)),
        S1=get_posterior_mean(survfit,pars="S1")[,"mean-all chains"],
        S1_l = summary(survfit,pars="S1")$summary[,"2.5%"],
        S1_u = summary(survfit,pars="S1")$summary[,"97.5%"],
        S2=get_posterior_mean(survfit,pars="S2")[,"mean-all chains"],
        S2_l = summary(survfit,pars="S2")$summary[,"2.5%"],
        S2_u = summary(survfit,pars="S2")$summary[,"97.5%"])

p1 <-ggplot(gehan3,aes(x=time))+
  geom_step(aes(y=S1,colour="6-MP"),size=1)+
  geom_ribbon(aes(ymin=S1_l,ymax=S1_u,fill="6-MP",alpha=0.1))+
  geom_step(aes(y=S2,colour="control"),size=1)+
  geom_ribbon(aes(ymin=S2_l,ymax=S2_u,fill="control",alpha=0.1))+
  theme(legend.title = element_blank())+
  guides(fill=FALSE,alpha=FALSE)+
  scale_color_manual(values=c("royalblue","tomato"))+
  scale_fill_manual(values=c("royalblue","tomato"))+
  ylab("survival")
print(p1)

追記

カプラン・マイヤープロットとの比較を忘れていた. 実線がカプラン・マイヤー推定量, 点線が提案法.

f:id:abrahamcow:20171011215610p:plain

これを見ると提案手法はちょっと平滑化しすぎかなあ.

library(survival)
kmfit <-survfit(Surv(time,cens)~treat,data=gehan)
plot(kmfit,col=c("royalblue","tomato"),lwd=2)
lines(gehan3$time,gehan3$S1,col="royalblue",lty=2,lwd=2)
lines(gehan3$time,gehan3$S2,col="tomato",lty=2,lwd=2)

アマゾンアフィリエイトのコーナー

生存時間分析の本, なにかいいのあったら教えてください.

赤澤『サバイバルデータの解析』は分厚すぎず, いいんじゃないかと思う.

Cox比例ハザードモデル (医学統計学シリーズ)

Cox比例ハザードモデル (医学統計学シリーズ)

中村『Cox比例ハザードモデル』もけっこう好き.