生存時間分析をStanで実行してみた

#はじめに

• 生存時間分析について説明変数なしモデルとありモデルをStanで作成。
• ハザード(瞬間死亡率)の対数値についてランダムウォークモデルを仮定。
• 参考文献 Rと生存時間分析（１）
• 参考文献 Rと生存時間分析（２）
• 参考文献 生存時間分析 - ハザード関数に時間相関の制約を入れる

#確率分布の定義

• 生存時間の確率分布$f_t(t=1,\cdots\infty)$
• 死亡関数(時刻$t$までに死亡している確率) $F_t=\sum_{i=1}^tf_i$
• 生存関数(時刻$t$まで生存している確率) $S_t=1-F_t=\sum_{i=t+1}^\infty f_i$
• 瞬間死亡率(ハザード) $h_t=\frac{f_t}{S_{t-1}}\quad(0\le h_t\le1, S_0=1)$

\begin{align}
f_t &= S_{t-1} h_t = F_t-F_{t-1}= S_{t-1}-S_t\\
S_t &= S_{t-1} (1-h_t)\\
S_t &= \prod_{i=1}^t(1-h_i),\quad
f_t = h_t\prod_{i=1}^{t-1}(1-h_i)
\end{align}

• 時刻$t$でイベントが発生する確率は$f_t$となる。
• 時刻$t$で打ち切りの場合、イベントは時刻$t+1$以降に発生するので、確率は$S_t$となる。

#説明変数なしモデル

##Stanプログラム

transformed parametersブロックで瞬間死亡率$h_t$から生存関数$S_t$と生存時間の確率分布$f_t$を計算する。

survival_0.stan

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1, upper=T> Time[N];
  int<lower=0, upper=1> Cens[N];
}

parameters {
  real<lower=0> sigma_log_hazard;
  vector<upper=0>[T] log_hazard;
}

transformed parameters {
  vector[T] log_S;
  vector[T] log_f;
  log_S[1] = log1m_exp(log_hazard[1]);
  log_f[1] = log_hazard[1];
  for (t in 2:T) {
    log_S[t] = log_S[t-1] + log1m_exp(log_hazard[t]);
    log_f[t] = log_S[t-1] + log_hazard[t];
  }
}

model {
  sigma_log_hazard ~ cauchy(0, 1);
  log_hazard ~ cauchy(0, 10);
  for (t in 2:T) {
    log_hazard[t] ~ normal(log_hazard[t-1], sigma_log_hazard);
  }
  for (n in 1:N) {
    if (Cens[n] == 1) {
      target += log_f[Time[n]];
    }
    else {
      target += log_S[Time[n]];
    }
  }
}

generated quantities {
  vector[T] hazard;
  vector[T] S;
  hazard = exp(log_hazard);
  S = exp(log_S);
}

##分析対象データ

パッケージMASSに含まれているgehanという生存時間データを分析する。

library(ggplot2)
library(scales)
library(ggfortify)
library(gridExtra)
library(rstan)
library(survival)
library(MASS)
library(muhaz)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

> head(gehan)
  pair time cens   treat
1    1    1    1 control
2    1   10    1    6-MP
3    2   22    1 control
4    2    7    1    6-MP
5    3    3    1 control
6    3   32    0    6-MP
> nrow(gehan)
[1] 42

##Stanの実行

treatが6-MPとcontrolについて別々に実行する。

model_stan = stan_model("survival_0.stan")
tstart = proc.time()
treats = c("6-MP","control")
result = lapply(seq_along(treats), function(i){
  ge = gehan[gehan$treat == treats[i],]
  data_stan = list(N=nrow(ge), T=max(ge$time), Time=ge$time, Cens=ge$cens)
  stan_fit = sampling(model_stan, data=data_stan,
    chains=4, iter=4000, warmup=2000, thin=1, seed=3,
    control=list(adapt_delta=0.8, max_treedepth=10)
  )
  stan_fit
})
proc.time() - tstart # 23sec

##実行結果の抽出

生存関数と瞬間死亡率の95%信用区間を抽出する。

my_extract = function(x,name){
  tb = summary(x, probs=c(0.025,0.975),pars=name)$summary
  df = data.frame(tb[,c("mean","2.5%","97.5%"), drop=F])
  colnames(df) = c("mean","lower","upper")
  df
}

df_S =lapply(result,function(x){
  df = my_extract(x,"S")
  df$time = 1:nrow(df)
  df
})

df_h =lapply(result,function(x){
  df = my_extract(x,"hazard")
  df$time = 1:nrow(df)
  df
})

##グラフの作成

生存関数と瞬間死亡率をsurvfit関数による結果と比較する。
グラフの色を合わせて表示したいのでggplot2が使用する色コードを取得。

# ggplot2のデフォルトの色を知りたい
# https://qiita.com/hoxo_b/items/c569da6dbf568032e04a
ggColorHue <- function(n, l=65) {
  hues <- seq(15, 375, length=n+1)
  hcl(h=hues, l=l, c=100)[1:n]
}
ggcols = ggColorHue(n=2)

kphaz = kphaz.fit(time=gehan$time, status=gehan$cens, strata=gehan$treat)
df_kphaz = lapply(1:2,function(i){
  f = (kphaz$strata==i)
  data.frame(time=kphaz$time[f],haz=kphaz$haz[f])
})

fit_sv = survfit(Surv(time, cens) ~ treat, data = gehan)
summary(fit_sv)

plot_S = autoplot(fit_sv) +
  geom_line(data=df_S[[1]],aes(x=time,y=mean),size=0.5,linetype=2,color=ggcols[1]) +
  geom_ribbon(data=df_S[[1]],aes(x=time,ymin=lower,ymax=upper),alpha=0.3,fill=ggcols[1],inherit.aes=F) +
  geom_line(data=df_S[[2]],aes(x=time,y=mean),size=0.5,linetype=2,color=ggcols[2]) +
  geom_ribbon(data=df_S[[2]],aes(x=time,ymin=lower,ymax=upper),alpha=0.3,fill=ggcols[2],inherit.aes=F)

plot_haz = ggplot() +
  geom_line(data=df_kphaz[[1]],aes(x=time,y=haz),size=0.5,linetype=1,color=ggcols[1]) +
  geom_line(data=df_kphaz[[2]],aes(x=time,y=haz),size=0.5,linetype=1,color=ggcols[2]) +
  geom_line(data=df_h[[1]],aes(x=time,y=mean),size=0.5,linetype=2,color=ggcols[1]) +
  geom_ribbon(data=df_h[[1]],aes(x=time,ymin=lower,ymax=upper),alpha=0.3,fill=ggcols[1]) +
  geom_line(data=df_h[[2]],aes(x=time,y=mean),size=0.5,linetype=2,color=ggcols[2]) +
  geom_ribbon(data=df_h[[2]],aes(x=time,ymin=lower,ymax=upper),alpha=0.3,fill=ggcols[2]) +
  scale_y_continuous(labels = scales::percent)

grid.arrange(plot_haz,plot_S,ncol=2,widths = c(0.8,1))

実線がsurvfit関数で推定した結果、破線がStanで推定した結果。

#説明変数ありモデル

• 対象とする説明変数は0または1の値を取るものとする。
• Cox比例ハザードモデルを適用する。説明変数の値=1の効果を$\beta$とすると、$\beta$は時刻に依存せず、瞬間死亡率を$\beta$倍にすると仮定する。
• データが長期間に渡る場合、全ての時刻で瞬間死亡率を推定すると計算時間がかかるため、イベント発生時刻・打ち切り時刻の間では、瞬間死亡率は一定と仮定する。
• $S_l=S_k(1-h_l)^{l-k}, f_l=S_k(1-h_l)^{l-k-1}h_l$
• ここで$h_l$は、時刻$k+1$から$l$までの瞬間死亡率とする。
• 計算時間を短縮するため説明変数の値の組合せでユニークなものを予めグループとしてまとめておく。データには対応するグループの番号を指定する。

##Stanプログラム

transformed parametersブロックで瞬間死亡率$h_t$からグループ毎の生存関数$S_t$と生存時間の確率分布$f_t$を計算する。
generated quantitiesブロックでベースラインを計算する。

survival_4.stan

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1> G;
  int<lower=1> K;
  int<lower=1, upper=T> TimeID[N];
  int<lower=0, upper=1> Cens[N];
  int<lower=0, upper=G> GroupID[N];
  int<lower=1> Time[T];
  matrix<lower=0, upper=1>[G, K] Group;
}

parameters {
  real<lower=0> sigma_log_hazard;
  vector<upper=0>[T] log_hazard;
  vector[K] beta;
}

transformed parameters {
  matrix[G, T] log_S;
  matrix[G, T] log_f;
  for (g in 1:G) {
    vector[T] new_hazard = log_hazard + dot_product(Group[g], beta);
    log_S[g, 1] = log1m_exp(new_hazard[1]) * (Time[1]);
    log_f[g, 1] = log1m_exp(new_hazard[1]) * (Time[1] - 1) + new_hazard[1];
    for (t in 2:T) {
      log_S[g, t] = log_S[g, t-1] + log1m_exp(new_hazard[t]) * (Time[t] - Time[t-1]);
      log_f[g, t] = log_S[g, t-1] + log1m_exp(new_hazard[t]) * (Time[t] - Time[t-1] - 1) + new_hazard[t];
    }
  }
}

model {
  sigma_log_hazard ~ cauchy(0, 1);
  log_hazard ~ cauchy(0, 10);
  beta ~ cauchy(0, 10);
  for (t in 2:T) {
    log_hazard[t] ~ normal(log_hazard[t-1], sigma_log_hazard);
  }
  for (n in 1:N) {
    if (Cens[n] == 1) {
      target += log_f[GroupID[n], TimeID[n]];
    }
    else {
      target += log_S[GroupID[n], TimeID[n]];
    }
  }
}

generated quantities{
  vector[T] hazard;
  vector[T] log_S0;
  vector[T] S0;
  vector[K] exp_beta;
  hazard = exp(log_hazard);
  log_S0[1] = log1m_exp(log_hazard[1]) * (Time[1]);
  for (t in 2:T) {
    log_S0[t] = log_S0[t-1] + log1m_exp(log_hazard[t]) * (Time[t] - Time[t-1]);
  }
  S0 = exp(log_S0);
  exp_beta = exp(beta);
}

##分析対象データ

パッケージsurvivalに含まれているkidneyという生存時間データを分析する。
説明変数としてsexとdiseaseを使用する。

> head(kidney)
  id time status age sex disease frail
1  1    8      1  28   1   Other   2.3
2  1   16      1  28   1   Other   2.3
3  2   23      1  48   2      GN   1.9
4  2   13      0  48   2      GN   1.9
5  3   22      1  32   1   Other   1.2
6  3   28      1  32   1   Other   1.2
> nrow(kidney)
[1] 76

##時刻情報とグループ情報の抽出

df_kidney = kidney
df_kidney$sex = c("M","F")[kidney$sex]
df_kidney$sex = as.factor(df_kidney$sex)

dum = cbind(
  model.matrix(~df_kidney$sex)[,-1,drop=F],
  model.matrix(~df_kidney$disease)[,-1]
)
colnames(dum) = c("sexM","diseaseGN","diseaseAN","diseasePKD")

times = sort(unique(df_kidney$time))
time_id = match(df_kidney$time,times)

group = unique(dum)
rownames(group) = 1:nrow(group)
group_char = apply(group,1,paste,collapse="")
group_id = match(apply(dum,1,paste,collapse=""),group_char)

> times
 [1]   2   4   5   6   7   8   9  12  13  15  16  17  22  23  24  25  26  27  28  30  34  38  39
[24]  40  43  46  53  54  58  63  66  70  78  96 108 113 114 119 130 132 141 149 152 154 156 159
[47] 177 185 190 196 201 245 292 318 333 402 447 511 536 562
> time_id
 [1]  6 11 14  9 13 19 57 54 20  8 15 52  5  7 58 20 27 50 10 44  5 55 41  6 34 22 42 32 59 16 12
[32]  2 48 47 53 37 13 46 10 35 43 60 56 15  9 31 23 26  8 24 36 51 40 45 21 20  1 16 39 17 18 29
[63]  3 25 43 20 49  3 38  6 28 11  4 33 30  6
> group
  sexM diseaseGN diseaseAN diseasePKD
1    1         0         0          0
2    0         1         0          0
3    0         0         0          0
4    1         1         0          0
5    0         0         1          0
6    1         0         1          0
7    1         0         0          1
8    0         0         0          1
> group_id
 [1] 1 1 2 2 1 1 3 3 1 1 3 3 4 4 2 2 5 5 4 4 5 5 3 3 5 5 5 5 3 3 6 6 3 3 3 3 2 2 3 3 7 7 3 3 5 5 5
[48] 5 6 6 5 5 2 2 5 5 4 4 2 2 5 5 5 5 8 8 2 2 3 3 3 3 8 8 7 7

##Stanの実行

model_stan = stan_model("survival_4.stan")

data_stan = list(
  N = nrow(df_kidney),
  T = length(times),
  G = nrow(group),
  K = ncol(group),
  TimeID = time_id,
  Cens = df_kidney$status,
  GroupID = group_id,
  Time = times,
  Group = group
)

tstart = proc.time()
stan_fit = sampling(model_stan,data=data_stan,
  chains=4, iter=4000, warmup=2000, thin=1, seed=3,
  control=list(adapt_delta=0.8,max_treedepth=10)
)
proc.time() - tstart # 122sec

##実行結果の抽出

生存関数と瞬間死亡率、説明変数の効果$\beta$の95%信用区間を抽出する。
瞬間死亡率の時刻情報を推定した区間の中央に設定する。

my_extract = function(x,name){
  tb = summary(x, probs=c(0.025,0.975),pars=name)$summary
  df = data.frame(tb[,c("mean","2.5%","97.5%"), drop=F])
  colnames(df) = c("mean","lower","upper")
  df
}

df_S0 = my_extract(stan_fit,"S0")
df_S0$time = times

df_h = my_extract(stan_fit,"hazard")
df_h$time = (function(){
  n = length(times)
  (times + c(0,times[-n]) + 1)/2
})()

df_exp_beta = my_extract(stan_fit,"exp_beta")
rownames(df_exp_beta) = c("sexM","diseaseGN","diseaseAN","diseasePKD")

beta_hist = extract(stan_fit,pars="exp_beta")$exp_beta
beta_hist = data.frame(beta_hist)
colnames(beta_hist) = c("sexM","diseaseGN","diseaseAN","diseasePKD")

##グラフの作成

推定された説明変数の効果$\beta$の分布。

pvline = geom_vline(xintercept=1,color="red")
plot1 = ggplot(beta_hist,aes(x=sexM)) + geom_histogram(bins=30) + pvline
plot2 = ggplot(beta_hist,aes(x=diseaseGN)) + geom_histogram(bins=30) + pvline
plot3 = ggplot(beta_hist,aes(x=diseaseAN)) + geom_histogram(bins=30) + pvline
plot4 = ggplot(beta_hist,aes(x=diseasePKD)) + geom_histogram(bins=30) + pvline

grid.arrange(plot1,plot2,plot3,plot4,ncol=2,nrow=2)

coxph関数による効果$\beta$の推定値との比較。

kidney.cox = coxph(Surv(time, status) ~ sex+disease, data=df_kidney)
summary(kidney.cox)
df_exp_beta

> kidney.cox = coxph(Surv(time, status) ~ sex+disease, data=df_kidney)
> summary(kidney.cox)
Call:
coxph(formula = Surv(time, status) ~ sex + disease, data = df_kidney)

  n= 76, number of events= 58 

              coef exp(coef) se(coef)      z Pr(>|z|)    
sexM        1.4774    4.3815   0.3569  4.140 3.48e-05 ***
diseaseGN   0.1392    1.1494   0.3635  0.383   0.7017    
diseaseAN   0.4132    1.5116   0.3360  1.230   0.2188    
diseasePKD -1.3671    0.2549   0.5889 -2.321   0.0203 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

           exp(coef) exp(-coef) lower .95 upper .95
sexM          4.3815     0.2282   2.17689    8.8188
diseaseGN     1.1494     0.8700   0.56368    2.3437
diseaseAN     1.5116     0.6616   0.78245    2.9202
diseasePKD    0.2549     3.9238   0.08035    0.8084

Concordance= 0.696  (se = 0.045 )
Rsquare= 0.206   (max possible= 0.993 )
Likelihood ratio test= 17.56  on 4 df,   p=0.001501
Wald test            = 19.77  on 4 df,   p=0.0005533
Score (logrank) test = 19.97  on 4 df,   p=0.0005069

> df_exp_beta
                mean      lower     upper
M          4.3978555 2.06594992 8.0631262
diseaseGN  0.9693798 0.45369567 1.7495952
diseaseAN  1.3627789 0.68891947 2.4198428
diseasePKD 0.2758943 0.08059133 0.6748802

df_cox = data.frame(summary(kidney.cox)$conf.int)
colnames(df_cox) = c("mean","inv_mean","lower","upper")
df_cox$x = factor(1:4,levels=1:4,labels=c("sexM","diseaseGN","diseaseAN","diseasePKD"))
df_exp_beta$x = (1:4)+0.1

ggplot() + 
  geom_pointrange(data=df_cox,aes(x=x,y=mean,ymin=lower,ymax=upper),size=0.5) +
  geom_pointrange(data=df_exp_beta,aes(x=x,y=mean,ymin=lower,ymax=upper),size=0.5,linetype=2) +
  labs(x="",y="beta")

実線がcoxph関数で推定した結果、破線がStanで推定した結果。

推定された生存関数と瞬間死亡率のグラフ。

plot_S = ggplot() +
  geom_line(data=df_S0,aes(x=time,y=mean),size=0.5,linetype=2) +
  geom_ribbon(data=df_S0,aes(x=time,ymin=lower,ymax=upper),alpha=0.3,inherit.aes=F) +
  labs(y="surv") +
  scale_y_continuous(labels = scales::percent)

plot_h = ggplot() +
  geom_line(data=df_h,aes(x=time,y=mean),size=0.5,linetype=2) +
  geom_ribbon(data=df_h,aes(x=time,ymin=lower,ymax=upper),alpha=0.3,inherit.aes=F) +
  labs(y="hazard") +
  scale_y_continuous(labels = scales::percent)

grid.arrange(plot_h,plot_S,ncol=2)

coxph関数によるベースラインとの比較。

kidney.fit = survfit(kidney.cox, newdata = data.frame(sex="F", disease="Other"))

# rough estimation of hazard from survival function
df_kphaz = (function(){
  n = length(kidney.fit$surv)
  times = (kidney.fit$time[2:n] + kidney.fit$time[1:(n-1)])/2
  interval = (kidney.fit$time[2:n] - kidney.fit$time[1:(n-1)])
  ratio = kidney.fit$surv[2:n] / kidney.fit$surv[1:(n-1)]
  data.frame(
    time = times,
    hazard = (1 - ratio) / interval
  )
})()

plot_S = autoplot(kidney.fit) +
  geom_line(data=df_S0,aes(x=time,y=mean),size=0.5,linetype=2) +
  geom_ribbon(data=df_S0,aes(x=time,ymin=lower,ymax=upper),alpha=0.3,inherit.aes=F)

plot_h = ggplot() +
  geom_line(data=df_kphaz,aes(x=time,y=hazard),size=0.5) +
  geom_line(data=df_h,aes(x=time,y=mean),size=0.5,linetype=2) +
  geom_ribbon(data=df_h,aes(x=time,ymin=lower,ymax=upper),alpha=0.3,inherit.aes=F) +
  labs(y="hazard") +
  scale_y_continuous(labels = scales::percent)

grid.arrange(plot_h,plot_S,ncol=2)

実線がcoxph関数で推定した結果、破線がStanで推定した結果。

#おわりに

• geom_pointrange関数でx座標を上手くずらす方法を見つけられなかった。
• stan_ac関数で確認すると、betaの自己相関は急速に減少しているが、sigma_log_hazardとlp__については、なかなか減少していない。Rhatは全て1.1以下となっている。
• Stanの実行パラメータはデフォルトのままなので調整の余地はある。

v = c("sigma_log_hazard","beta","lp__")
stan_ac(stan_fit,pars=v)

print(stan_fit,digits=3,pars=v,probs=c(0.025,0.5,0.975))

Inference for Stan model: survival_4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

                     mean se_mean     sd     2.5%      50%    97.5% n_eff  Rhat
sigma_log_hazard    0.203   0.012  0.109    0.053    0.185    0.452    89 1.059
beta[1]             1.422   0.005  0.345    0.726    1.428    2.087  5849 1.001
beta[2]            -0.090   0.006  0.346   -0.790   -0.081    0.559  3299 1.001
beta[3]             0.257   0.007  0.324   -0.373    0.257    0.884  2046 1.002
beta[4]            -1.431   0.008  0.545   -2.518   -1.414   -0.393  4618 1.000
lp__             -174.200   2.727 31.170 -226.916 -177.323 -106.646   131 1.071