生存時間分析の基礎３（生存曲線の検定）

前回の記事を up してからかなりの期間をあけてしまいましたが，ようやく重い腰をあげて続編を up しました．
ペースの遅い投稿者ですみません（汗．このような感じで遅々としてやっていくと思いますが，気がむいた際にでも読んでいただけますと幸いです．

ちなみに，続編を執筆せずに何をやっていたのかといいますと，昨年の 11 月から 今年 1 月上旬にかけて開催されていた ProbSpace の「浮世絵作者予測」や 昨年末から今年 3 月上旬まで開催されていた Kaggle の 「Bengali.AI Handwritten Grapheme Classification 」というコンペに参加していました.

「浮世絵作者予測」は名前の通り，色々な浮世絵の作者（全10人）を 224 x 224 ピクセルの RGB 画像から推定するというものでした．結果は 1 位の方（ solution ）には遠く及ばず 2 位という結果に終わりましたが，浮世絵は普段まじまじと眺めることもありませんので，新鮮な体験をすることができ良かったと思っています．
私の solution 概要と SpeakerDeck の資料は以下となりますので，よろしければご覧ください．

ProbSpace Discussion: 

2nd place solution | ProbSpace

Speaker Deck:

※・・・リーダーボードや投稿された solution などをみるには ProbSpace への登録と log in が必要です．興味のある方は是非登録してみてください．

また，「Bengali.AI Handwritten Grapheme Classification 」は，ベンガル文字という日本人にとってはあまり馴染みのない文字を， 137 x 236 ピクセルのモノクロ画像から分類するタスクでした．簡単ではありますが，solution を以下のリンク先に up していますので，興味のある方はご覧ください．

Kaggle discussion: 

Bengali.AI Handwritten Grapheme Classification | Kaggle

Speaker Deck (High resolution model pipeline figure):

さて，また前置き（投稿遅延の言い訳）が長くなってしまいましたが，今回は前回の Kaplan-Mier 生存曲線 に続いて，生存曲線の検定を確認していきたいと思います．

1.   生存曲線の検定

生存曲線の検定を考える場合，もっとも代表的なものはログランク検定です．ログランク検定とは，2 群もしくは複数群の KM 曲線が統計学的に差がないかを判断する 際に使用される検定です．
今回は，2 群のログランク検定に関して説明をします．3 群以上のログランク検定 や ログランク検定以外の検定 もありますが，興味のある方は「エモリー大学クラインバウム教授の生存時間解析」を勉強してみることをお勧めします．

エモリー大学クラインバウム教授の生存時間解析

では，2 群のログランク検定をみていきましょう．

2.   2 群のログランク検定のための前準備

ログランク検定は，「KM 曲線全体」を比較する統計量を計算することで行います．比較する 2 つの生存曲線に統計学的に差がないとは，2 つの曲線を全体的に比較する検定（今回はログランク検定）において，それぞれの真の生存曲線間に差がある，という証拠がないことを意味します．検定に使用する統計量は，サンプル数が多い大標本の仮定のもとで，χ ^2 分布（カイ二乗分布）に従うことになります．これからその計算の仕方を追ってみますが，「独立性の χ ^ 2 検定」を知っている方は理解がしやすいでしょう．
（実際にログランク検定を行う場合は，以下で計算する手順を踏まなくても，少数のコードで簡単に検定結果を得ることができます．最後の方で紹介をします．）

さて，前回同様に gehan データを使用してみていきますが，まず最初に KM 曲線を計算するところまで行います．
（Step 1~4 の内容は前回とほぼ変わりませんので，今回は説明は省略しますが，忘れてしまっていたり確認をしておきたい場合は，下の前回記事へのリンクをたどってみてください．）

library(MASS)
library(survival)
library(survminer)
library(tidyverse)

# Step 1
data(gehan)
gehan <- gehan %>% 
   mutate(treat = ifelse(treat == "6-MP", "6MP", "control"))

# Step 2
surv.obj <- Surv(time = gehan$time, event = gehan$cens)

# Step 3
ge.sf <- survfit(surv.obj ~ treat, data = gehan)

# Step 4
ge.sf.df <- surv_summary(ge.sf, data = gehan)
View(ge.sf.df)

次に，ログランク検定の計算を自分で行うために，上のコード内の Step 4 で得られた ge.sf.df データフレームを使い，時系列順に介入群（投薬有）と対照群（投薬無）及び 全体 の各時点における生存数とイベント数を表すデータフレームを以下の Step 5 で作成します．

# Step 5

ge.table <- ge.sf.df %>% 
   select(time, n.risk, n.event, n.censor, treat) %>% 
   pivot_wider(names_from = treat, 
               values_from = c(n.risk, n.event, n.censor)) %>% 
   arrange(time) %>% 
   fill(starts_with("n.risk"), .direction = "up") %>% 
   fill(starts_with("n.risk"), .direction = "down") %>% 
   replace_na(list(n.event_6MP = 0, n.event_control = 0, 
                   n.censor_6MP = 0, n.censor_control = 0)) %>% 
   mutate(n.event = n.event_6MP + n.event_control,
          n.risk = n.risk_6MP + n.risk_control)
View(ge.table)

そうすると，下図 10 のようなデータフレームが得られるはずです．

各カラムの意味は以下の通りです．

time : イベントや打ち切りが発生した時点（生存時間）
n.risk_X : その時点での X 群の生存数（リスクセットとも呼ばれます）
n.event_X : その時点での X 群のイベント発生数（gehan データの場合は死亡数）
n.censor_X : その時点での X 群の打ち切り数
n.event : その時点での全体の生存数
n.risk : その時点での全体のイベント発生数

ここまでは，各時点での各群・全体の生存数・イベント発生数をまとめているだけで，特別なことは行っていませんが，次に説明する手順・考え方が重要となります．

ログランク検定では，帰無仮説は「2 つの群の生存曲線に差はない」であり，対立仮説は「2 つの群の生存曲線に差がある」となります．

帰無仮説に相当する 2 つの群の生存曲線に差がない場合，各時点でのイベントの発生確率は同じなわけですから，イベント数はそれぞれの群の生存数に比例したもの（期待イベント発生数）になると考えられます．それを式で表すと以下のようになります．

$$
e_{1f}=n_{1f} \times \frac{m_{1f} \; + \; m_{2f}}{n_{1f} \; + \; n_{2f}} \\

e_{2f}=n_{2f} \times \frac{m_{1f} \; + \; m_{2f}}{n_{1f} \; + \; n_{2f}}
$$

この数式中で使用されている記号の意味は以下となります．

$${e_{1f}}$$ , $${e_{2f}}$$: 各群（1: 介入群 2: 対照群）の時点 $${f}$$ における期待イベント発生数
$${n_{1f}}$$ , $${n_{2f}}$$: 各群の時点 $${f}$$ における生存数
$${m_{1f}}$$ , $${m_{2f}}$$: 各群の時点 $${f}$$ におけるイベント発生数

ログランク検定においては，この期待イベント発生数を使い，独立性の χ^2 検定に似た形で，検定に使用する統計量を計算していきます．
ここでは，2 種類の統計量の計算の仕方を紹介します．

3.   超幾何分布にもとづいたログランク統計量の計算方法

最初に紹介する方法は 超幾何分布 を利用した方法で，先ほどの数式の記号を使い統計量を表すと，以下のようになります．

統計量 $${T_1}$$ は 2 つの群（介入群と対照群）間に制約式が存在するため，帰無仮説の元で自由度 1 の $${\chi^2}$$ 分布に従うことになります ( $${\chi^2}$$ 検定における自由度の考え方を思い出すとよいと思います)．ここで，$${T_1}$$ の分子は「観測量の合計である $${O_i}$$ から，帰無仮説の元での期待値の合計である $${E_i}$$ を引いたものを二乗」しています．一方で，分母は「$${O_i}$$ の分散の値」となっているわけですから，中心極限定理を仮定すれば，$${T_1}$$ は正規分布に従う変数の二乗，つまり $${χ^2}$$ 分布に従う変数と考えてよいわけです．
$${O_i}$$ の分散が上式のようになっているのは，前述したように $${O_i}$$ が 超幾何分布 に従うと考えることができるからです．超幾何分布の定義はリンク先をみていただくとして，超幾何分布の期待値と分散の式は以下のようになります（wikipedia より転載）．

上式内の記号の対応関係は以下となります．代入してみると，同じ式になることが確認できると思います．

$$
N = n_{1f} + n_{2f} \\
K = n_{1f} \\
N - K = n_{2f} \\
n = m_{1f} + m_{2f} \\
$$

では，Step 5 で作成したデータフレームを利用して，以下のコードで $${T_1}$$ 統計量を算出するためのカラムをいくつか作成しましょう．

# Step 6

ge.table <- ge.table %>% 
   mutate(
       n.event.expected_6MP = n.event * n.risk_6MP / n.risk,
       n.event.expected_control = n.event * n.risk_control / n.risk
       ) %>% 
   mutate(
       o_e_6MP = n.event_6MP - n.event.expected_6MP,
       o_e_control = n.event_control - n.event.expected_control,
       var = n.risk_6MP * n.risk_control * n.event * (n.risk - n.event) /
           n.risk ^ 2 / (n.risk - 1)
   )

各カラムの意味は以下の通りです．

n.event.expected_XXX : 各時点での期待イベント発生数
o_e_XXX : 各時点での T_1 統計量の分子にあたる値
var : 各時点での T_1 統計量の分母にあたる値

T_1 統計量の計算をする準備ができましたので，計算をしてみましょう．

# Step 7

sum(ge.table$o_e_6MP) ^ 2 / sum(ge.table$var)
# 16.79294

pchisq(16.79, df = 1, lower.tail = FALSE)
# 4.175275e-05

数式通りに計算をしてみると，~ 16.8 となりました．T_1 統計量は，自由度 1 の χ^2 分布に従う統計量でしたので，p 値を計算してみると ~ 4.2 e-5 となりました．結果的に，代表的な有意水準のレンジ（0.01 ~ 0.05）では，2 つの生存曲線に差がないという帰無仮説は棄却されることになりそうです．

4.   χ^2 検定と似たログランク統計量の計算方法

もう 1 つの統計量は，χ^2 検定とよく似た形で計算する方法です．数式で表すと以下のようになります．

χ^2 検定の統計量の計算の仕方とよく似ているのが分かると思います．こちらも自由度 1 の χ^2 分布に従うことになりますが，理論的にどうしてそのようになるかを知りたい方は，下のリンクをみてみると良いかもしれません．
（ところどころに数式の表記ミスが見受けられますが，証明の流れは勉強になると思います．なお，作者は阪大の鈴木教授で，先日出版された「機械学習の数理 100 問シリーズ」の著者でもある方です．）

タイトル未設定

いずれにしても先ほどよりも簡単に計算できそうですので，以下のコードのように計算をしてみましょう．

# Step 8

t2.1 <- sum(ge.table$o_e_6MP) ^ 2 / sum(ge.table$n.event.expected_6MP)
t2.2 <- sum(ge.table$o_e_control) ^ 2 / sum(ge.table$n.event.expected_control)
t2.1 + t2.2
# 15.23285

pchisq(15.23, df = 1, lower.tail = FALSE)
# 9.517935e-05

統計量の値が先ほどのやり方と比較して少し小さくなりましたが，p 値は十分に小さいので，帰無仮説が棄却されるという結果には変わりはなさそうです．

5.   survdiff 関数を使ったログランク検定

ここまで，敢えて KM 曲線から計算をしてきましたが，実際にログランク検定をする場合は，survdiff 関数を使うことで簡単に検定を行うことができます．コードは以下のようになります．

# Step 9

survdiff(surv.obj ~ treat, data = gehan)

たった 1 行のコードを実行すると上記のような結果が得られるはずです．各統計量の結果が先ほど計算したものと一致していることや，ここで得られる p 値が先に紹介した計算方法（超幾何分布の方法）による p 値（= 4.175275e-05）と一致していることが分かります．

いかがでしたでしょうか．
今回は，たった 1 つの関数（ survdiff ）で結果が得られるために，おざなりになりがちなログランク検定の統計量の計算の仕方を確認しました．

さて，本記事を書いている現時点（2020年5月初旬），世界的に感染拡大してしまった COVID-19 のため，私の住む東京では自宅での自粛が政府より要請されている状況です．自宅に待機し続けるのがつらいという方も多いとは思いますが，この際，生存時間分析を始めとした統計・分析の勉強に時間を費やすのも 1 つの選択肢かなと思います．こういった時間を有意義に使うことができれば，COVID-19 が収束し世の中がまた動き出した頃には，きっと大きな実力がついているのではないでしょうか．

長くなってしまいましたが，次回はいよいよモデリングということで，Cox 比例ハザードモデルの紹介をしたいと思います．
それでは！

6.   Codes

今回使用したコードは以下となります．

Greenwind1/survival-analysis

7.   References

エモリー大学クラインバウム教授の生存時間解析

ピアソンの定理の証明（阪大，鈴木教授）:
http://joe.bayesnet.org/wp-content/uploads/2017/05/2017-5-25-1.pdf

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)