生物医学統計を学ぶ！: 生存時間データを用いた群間比較

2群の比較

ある疾患を対象に、対照薬群と治療薬群とを比較するランダム化試験を行ったとする。

両群の生存時間分布に違いがあるかを検定したい。

```
  # https://drive.google.com/file/d/1WnbMnyzcXigGkQEzxb5jOmC5oULe9Pl9/view?usp=sharing
  # 模擬データをダウンロードします
  !curl -L -o "survdiff-test.csv" "https://drive.google.com/uc?export=download&id=1WnbMnyzcXigGkQEzxb5jOmC5oULe9Pl9"
  # データを読み込みます
  import pandas as pd
  # treatは薬の違いで分けた群、timeは生存期間（日数と仮定）、censorは打ち切り（=死亡）の有無
  data = pd.read_csv(filepath_or_buffer="survdiff-test.csv", names=["treat","time","censor"], encoding="utf-8")
  print(data.shape)
  print(data.head())
  print(data.describe())
  print(round(data[data["censor"]==1].shape[0]/data.shape[0],2)) # 23/32が死亡
  # 治療群を分ける
  t1_df = data[data["treat"]==1]
  t2_df = data[data["treat"]==2]
  print(t1_df.shape)
  print("control 打ち切り数",t1_df[t1_df["censor"]==1].shape)
  print(t2_df.shape)
  print("therapeutic 打ち切り数",t2_df[t2_df["censor"]==1].shape)
  # 生存確率を予測するモデルを定義
  # control群の予測器
  sf = sm.SurvfuncRight(t1_df["time"], t1_df["censor"])
  print(sf.summary())
  print(t1_df[t1_df["censor"]==1].shape)
  '''
  生存分布の分位数の点推定と信頼区間を取得できます。
  死亡した被験者の割合の分位数しか推定できません。
  （例えば、controlでは、13/16未満）
  '''
  print(sf.quantile(0.5)) # 50パーセンタイルの点推定値（50%の確率で22日間生存）
  print(sf.quantile_ci(0.5, alpha=0.05, method='cloglog')) # 信頼区間 
  # survival function, 生存関数のプロット
  fig = sf.plot()
  lcb, ucb = sf.simultaneous_cb(alpha=0.05, method='hw', transform='log')
  ax = fig.get_axes()[0]
  ax.fill_between(sf.surv_times, lcb, ucb, color='lightgrey')
  ax.set_ylabel("Proportion alive")
  ax.set_xlabel("Days since enrollment")
  # 群間比較
  gb = data.groupby("treat")
  ax = plt.axes()
  itr = 0
  for g in gb:
      sf = sm.SurvfuncRight(g[1]["time"], g[1]["censor"])
      sf.plot(ax)
      lcb, ucb = sf.simultaneous_cb()
      if itr == 0:
        ax.fill_between(sf.surv_times, lcb, ucb, alpha=0.2, color='lightblue')
      else :
        ax.fill_between(sf.surv_times, lcb, ucb, alpha=0.2, color='lightgray')
      itr += 1
  li = ax.get_lines()
  li[1].set_visible(False)
  li[3].set_visible(False)
  plt.figlegend((li[0], li[2]), ["control","therapeutic"], "center right")
  ax.set_ylabel("Proportion alive")
  ax.set_xlabel("Days since enrollment")
  # logrank test（原点が同じで、差が開いていくような分布の比較のときはログランクを使う）
  stat, pv = sm.duration.survdiff(data.time, data.censor, data.treat)
  print(pv)# 有意差あり。新しい治療のほうが治療成績（予後）が良い
    
  # 補足
  # CIの範囲が途中で切れるのは、その区間のサンプル数が1つしか無いため。
  print(data[data["time"]>50].shape)
  print(t1_df[t1_df["time"]>50].shape)
  print(t2_df[t2_df["time"]>62].shape)
  ```

References

http://arxiv.org/pdf/math/0409180.pdf
https://humboldt-wi.github.io/blog/research/information_systems_1920/group2_survivalanalysis/#deepsurv

Visionary Imaging Services, Inc.

Imaging CROサービスをご提供させていただきます。

お問い合わせお待ちしております。

https://www.vis-ionary.com/