Syleir’s note

2020.4.1より統計検定やE資格の勉強の進捗を報告しています。統計検定準1級、E資格、G検定取得しました！当ブログへのリンクはご自由にどうぞ。

MENU

【生存時間解析】ログランク検定の実践【Part6】【統計検定準1級•1級】

統計検定

はじめに

ついに想定している最終回です。
前回まではカプランマイヤー曲線の作成を行ったのでした。
syleir.hatenablog.com

今回は

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

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

Amazon
生存時間解析入門 原書第2版

生存時間解析入門 原書第2版

Amazon

これらをベースに咀嚼してから執筆しています。ぜひ原著にも当たってください。

目次

Kaplan-Meier法で生存曲線が書けるようになった

カプランマイヤー法で生存曲線を推定して描出したものがカプランマイヤー曲線でした。
カプランマイヤー曲線は、実際のデータから生存関数を推定したいというモチベーションで作成されます。治療群と対照群で生存曲線を作った時に、それらの生存曲線は全く同じ比でイベントや打ち切りが起きていない限り、当然違うものになります。
この生存曲線の違いが、統計的に有意なものかどうかを知りたいというのが、自然な次のモチベーションです。

ログランク検定

ログランク検定(Log-rank test)は「生存曲線に差があるか」をテーマとした検定です。
「生存曲線に差がある」とはどういうことでしょうか。

2本の生存曲線に差があるか？

明らかにこれらの生存曲線は当然違います。ここでいう、生存曲線に差があるかの「差」というのは、「2本の曲線がどれくらい違うか大域的に示すような値」のことであって視覚的な差のことではありません。
つまり、もう少し補うと、生存曲線に「統計的な」差があるかです。この統計処理にもたくさんの種類がありますが、その一つがログランク検定になります。


ログランク検定は生存曲線に差があるか？という問いに答える検定です。
検定というからには帰無仮説と対立仮説があります。
ログランク検定における帰無仮説は、「全時点で、2つの生存関数が等しい」です。
この帰無仮説に
・打ち切りは群と独立に起こること
・観測された総イベント数、期待されるイベント数が大きいこと
を仮定するとログランク検定が実現されます。

ログランク統計量

ログランク検定における、「2本の曲線がどれくらい違うか大域的に示すような値」のことを「ログランク統計量」といいます。
このログランク統計量は自由度1の\chi ^2分布に従うことが知られていますので、ログランク統計量の導き方さえわかれば、あとは統計検定2級、あるいは統計検定準1級前半の知識で解くことができます。

分割表によるカイ二乗検定を復習する

ワークブック28章にありますが、変数の値の組み合わせごとの頻度を表にしてまとめたものを分割表と言います。一旦脱線しますが、今後の解説のために一度復習しておきます。
既知の方はスキップして構いません。

生存 死亡 合計
新薬 16 4 20
プラセボ 12 8 20
合計 28 12 40

この図は2×2の分割表で、

\dfrac{(観測度数-期待度数)^2}{期待度数}
をすべての群で足し合わせたものが、
が、漸近的に自由度1の \chi^2分布に従うことが知られています。

期待度数とは？

期待度数というのは、先ほどのこの表を

生存 死亡 合計
新薬 20
プラセボ 20
合計 28 12 40

このように空欄にしたとき、空欄をもっともらしく埋めた値のことです。
例えば、新薬群のうち、生存したものは、40*(20/40)*(28/40)とするのが期待度数です。
同様に、他のものも埋めていくと、このようになります。

生存 死亡 合計
新薬 40\times \dfrac{20}{40} \times \dfrac{28}{40} 40\times \dfrac{20}{40} \times \dfrac{12}{40} 20
プラセボ 40\times \dfrac{20}{40} \times \dfrac{28}{40} 40\times \dfrac{20}{40} \times \dfrac{12}{40} 20
合計 28 12 40

となり、計算すると
期待度数としては、

生存 死亡 合計
新薬 14 6 20
プラセボ 14 6 20
合計 28 12 40

となります。
これから、 \chi^2 = \dfrac{(観測度数-期待度数)^2}{期待度数}
を計算すると、
\chi^2 = \dfrac{(16-14)^2}{14} +\dfrac{(12-14)^2}{14} + \dfrac{(4-6)^2}{6} + \dfrac{(8 - 6)^2}{6}
= \dfrac{4}{14} +\dfrac{4}{14} + \dfrac{4}{6} + \dfrac{4}{6}
= \dfrac{4}{14} +\dfrac{4}{14} + \dfrac{4}{6} + \dfrac{4}{6}
= \dfrac{4}{7} +\dfrac{4}{3}
= \dfrac{4}{7} +\dfrac{4}{3}
= \dfrac{40}{21} \fallingdotseq 1.90

となります。

観測値 B_1 B_2 合計
A_1 x_{11} x_{12} x_{1\cdot}
A_2 x_{21} x_{22} x_{2\cdot}
合計 x_{\cdot 1} x_{\cdot 2} x_{\cdot\cdot}

と抽象化した時に、
\chi^2 = \dfrac{x_{\cdot\cdot} (x_{11} x_{22} - x_{12} x_{21})^2}{x_{1\cdot}x_{2\cdot}x_{\cdot 1}x_{\cdot 2}}としても同じになることも知られています。

これらの \chi^2をピアソンのカイ２乗検定統計量と言います。
これが帰無仮説のもとでカイ２乗分布に分布収束することが示されています。
詳細は、

現代数理統計学の基礎 共立講座 数学の魅力

現代数理統計学の基礎 共立講座 数学の魅力

Amazon
の7.4 適合度検定をご参照ください。証明もあります。

この例では、
\chi^2 = \dfrac{40 \times (16 \times 8 - 4 \times 12)^ 2}{20 \times 20 \times 28 \times 12}
\fallingdotseq 1.90
として同じ値が出てきます。

より正確には、フィッシャーの正確検定により検定を行います。

ログランク検定とは一体何をやっているのか？

ログランク検定はノンパラメトリックな検定である

まず最初に理解する必要があるのは、ログランク検定は生存曲線の関数の形を仮定する必要のないノンパラメトリックな検定ということです。実際のデータをもとに、ログランク統計量という、「2本の曲線がどれくらい違うか大域的に示すような値」を算出しますが、ここに曲線がどのような関数であるかを指定する余地はありません。

ログランク統計量の計算方法

上に書いたイベント数が多いこと、という過程を早速無視して、手計算可能なレベルで以下の例からログランク検定をやってみましょう。

t 新薬 死亡数 打ち切り リスクセット プラセボ 死亡数 打ち切り リスクセット
1 0 0 6 2 0 6
2 2 0 6 1 0 4
3 0 0 4 1 0 3
4 0 1 4 0 0 2
5 1 0 3 1 0 2
6 1 0 2 0 0 1
7 0 0 1 0 0 1

この例を使ってログランク統計量を計算していきましょう。

①どちらかの群のイベント発生数と期待イベント発生数の差をイベント発生時点ごとに計算

イベント発生ごとに分割表を作成します。

t = 1

⚪︎イベント発生数

生存 死亡 合計
新薬 6 0 6
プラセボ 4 2 6
合計 10 2 12

⚪︎期待イベント発生数

生存 死亡 合計
新薬 6
プラセボ 6
合計 10 2 12

これを埋めて、

生存 死亡 合計
新薬 5 1 6
プラセボ 5 1 6
合計 10 2 12

となります。t = 1での新薬群の実際のイベント数と期待イベントの差は0 - 1 = -1です。

上の分割表のピアソンのカイ２乗検定統計量と違って、4箇所すべてを計算する必要はありませんが、常に同じところを計算し、実際のイベント数と期待発生数の引く順番も常に揃える必要があります。
今回は（死亡, 新薬）かつ実際のイベント数ー期待発生数で統一して計算します。
揃えてさえいれば、どの群、どちらの引き算で計算しても、正負が逆転するだけで、同じ絶対値の点数になるのでどちらでも構いません。

t = 2

⚪︎イベント発生数

生存 死亡 合計
新薬 4 2 6
プラセボ 3 1 4
合計 7 3 10

⚪︎期待イベント発生数

生存 死亡 合計
新薬 5.2 1.8 6
プラセボ 2.8 1.2 4
合計 7 3 10

t = 2での死亡群の実際のイベント数と期待イベントとの差は2 - 1.8 = 0.2です。

t = 3

⚪︎イベント発生数

生存 死亡 合計
新薬 4 0 4
プラセボ 2 1 3
合計 6 1 7

⚪︎期待イベント発生数

生存 死亡 合計
新薬 3.43 0.57 4
プラセボ 2.57 0.43 3
合計 6 1 7

t = 3での死亡群の実際のイベント数と期待イベントとの差は0 - 0.57 = -0.57です。

t = 4

t = 4は打ち切りはありますが、イベントはないので、計算を行いません。

t = 5

飽きた人はスクロールしてください。
⚪︎イベント発生数

生存 死亡 合計
新薬 2 1 3
プラセボ 1 1 2
合計 3 2 5

⚪︎期待イベント発生数

生存 死亡 合計
新薬 1.8 1.2 3
プラセボ 1.2 0.8 2
合計 3 2 5

t = 5での生存群の実際のイベント数と期待イベントとの差は1 - 1.2 = -0.2です。

t = 6

⚪︎イベント発生数

生存 死亡 合計
新薬 1 1 2
プラセボ 1 0 1
合計 2 1 3


⚪︎期待イベント発生数

生存 死亡 合計
新薬 1.33 0.67 2
プラセボ 0.67 0.33 1
合計 2 1 3

t = 6での生存群の実際のイベント数と期待イベントとの差は1 - 0.67 = 0.33です。

t = 7

イベントが発生していないので計算は起こりません。

②これらの値の総和をとる

次に、これらの値の総和を取ります。
この教科書で定義される通りに文字を定義していくと、
i = 1, 2(群：1 = 新薬, 2 = プラセボ)
f = 1, 2, 3, 4, 5(時間：イベントが起きたタイミングの昇順、今回はf=1,2,3,4,5がt=1,2,3,5,6に対応)
m_{if}(群i, 時間fでの実観測イベント数)
e_{if}(群i, 時間fでの期待イベント数)
となり、上で計算したことをまとめていくと次のようになります。

t f m_{1f} e_{1f} m_{1f} - e_{1f}
1 1 0 1 -1
2 2 2 1.8 0.2
3 3 0 0.57 -0.57
4 0 0 0
5 4 1 1.2 -0.2
6 5 1 0.67 0.33
7 0 0 0
合計 4 5.24 -1.24

さらに文字が増えてごちゃごちゃしますが、
「群iの、観測度数ー期待度数の値をすべてのイベントが起きた時間について合計したもの(-1.24のことです)」
を O_{i} - E_{i} = \sum_{f} (m_{if} - e_{if})と定義すると、
群1（新薬群）についての O_{1} - E_{1} = -1.24
となります。

③分散を計算する

ここまできたらログランク統計量まであと一息です。
定義としては、
Log-rank statistic = \dfrac{(O_1 - E_1)^2}{Var(O_1 - E_1)}
です。(O_1 - E_1) は計算済みですから、Var(O_1 - E_1)を計算すればよく、これは観察度数ー期待度数の分散のことです。

この分散は、
\hat{Var}(O_1 - E_1)
= \sum_f \dfrac{n_{1f}n_{2f}(m_{1f} + m_{2f})(n_{1f} + n_{2f} - m_{1f} - m_{2f})}{(n_{1f} + n_{2f})^2 (n_{1f} + n_{2f} - 1)}

として推定が可能です。こちらの教科書の教科書のシグマの対象の添字にはミスがありますので持っている方はご参照ください。
各イベントが発生した時間fについて、
n_{1f}：新薬群のリスクセットの人数
n_{2f}：プラセボ群のリスクセットの人数
m_{1f}：新薬群の死亡数
m_{2f}：プラセボ群の死亡数
です。

当てはめるだけなので簡単に計算できますが、手計算は人間の仕事ではなく、現実的ではないので、簡単に計算でき、出題が想定される例をご紹介します。

各時間でイベントが1人にしか起こらないとき

各fでイベントが1人にしか起こらないとは、すべてのfでm_{1f} + m_{2f} = 1のことです。

ここではm_{1f} = 1, m_{2f} = 0としましょう。
\dfrac{n_{1f}n_{2f}(m_{1f} + m_{2f})(n_{1f} + n_{2f} - m_{1f} - m_{2f})}{(n_{1f} + n_{2f})^2 (n_{1f} + n_{2f} - 1)}
= \dfrac{n_{1f}n_{2f} (n_{1f} + n_{2f} - 1)}{(n_{1f} + n_{2f})^2 (n_{1f} + n_{2f} - 1)}
= \dfrac{n_{1f}n_{2f} }{(n_{1f} + n_{2f})^2}
= \dfrac{n_{1f}}{n_{1f} + n_{2f}} \dfrac{n_{2f}}{n_{1f} + n_{2f}}

ここで、 p_f = \dfrac{n_{1f}}{n_{1f} + n_{2f}}と置くと、 1 - p_f = \dfrac{n_{2f}}{n_{1f} + n_{2f}}
となり、

= p_f(1 - p_f)
となり、これは二項分布 B(1, p_f)の分散に相当しますから、
各イベント発生時点でのイベント数が、１例の場合にはその群の期待イベント発生確率をパラメータとする２項分布の分散をイベント発生時点ごとに計算し、合計したものが観察度数ー期待度数の分散になる
ことがわかります。
これくらいであれば手計算でもやる気の出る範囲内であり、これは2016年の統計検定1級で既出です。

分散を期待度数で近似するとき

現実には、各fでイベントが1人にしか起こらないことはありません。
もっとも、プログラミングをできない環境にいることはさらにありえません。手計算でログランク検定をやらなくてはいけないのは、統計検定の会場にいる時か、インフラのない無人島に漂着した時くらいです。でも残念ながら、我々はこのケースに該当してしまうのでやらなければなりません。

一般にログランク統計量
Log-rank statistic = \dfrac{(O_1 - E_1)^2}{Var(O_1 - E_1)}は、分散を使わずに、
\fallingdotseq \sum_i \dfrac{(O_i - E_i)^2}{E_i}
= \dfrac{(O_1 - E_1)^2}{E_1} + \dfrac{(O_2 - E_2)^2}{E_2}（２群の場合）

近似できることが知られています。
O_1 - E_1 = O_2 - E_2
ですし、群すべてのイベント数の平均をEとすると
E = \dfrac{E_1 + E_2}{2} = \dfrac{O_1 + O_2}{2}
が成り立ちますから、
E_2 = O_1 + O_2 - E_1
= 総イベント数 - E_1
とするともう一度計算しなくても簡単に求まります。
この例では
E_2 = 9 - 5.24 = 3.76
です。

④ログランク統計量を計算し、検定する

Log-rank statistic \fallingdotseq \sum_i \dfrac{(O_i - E_i)^2}{E_i}
= \dfrac{(O_1 - E_1)^2}{E_1} + \dfrac{(O_2 - E_2)^2}{E_2}
= \dfrac{(-1.24)^2}{5.24} + \dfrac{(-1.24)^2}{3.76}
= \dfrac{(-1.24)^2}{5.24} + \dfrac{(-1.24)^2}{3.76}
\fallingdotseq 0.29 + 0.41 = 0.60

となります。近似したものにせよ、近似しないものにせよ、漸近的に
ログランク統計量は自由度1のカイ２乗分布\chi(1)に従うことが知られており、
\chi_{0.05}(1) = 3.84 < 0.60ですから、有意水準0.05で検定した時、
「本例での2本の生存曲線に差がないという帰無仮説は棄却されない」という結果になります。

ログランク統計量とピアソンのカイ２乗統計量を見比べてみる

ログランク統計量
Log-rank statistic \fallingdotseq \sum_i \dfrac{(O_i - E_i)^2}{E_i}
と近似できるのでした。

ピアソンのカイ２乗統計量についても
\chi^2 = \dfrac{(観測度数-期待度数)^2}{期待度数}
をすべての群で足し合わせたもの
ですから、全く同じ形をしています。
実は、ログランク統計量というのは、すべてのイベント発生時について、分割表を書き、総和を取っていった後に計算してみるとピアソンのカイ２乗統計量の形に近似できるので、自由度1のカイ２乗分布に従う、というように解釈できます。

おわりに

このパートで当初想定していた、ログランク検定までの初歩の記事執筆はおしまいになります。
余力のある方は比例ハザードモデルあたりも理解していただければと思います。
今回もなんとか失踪せずに書き上げることができました！
少しでも初歩の理解につながっていれば幸いです。

そういえば、なかなかわかりやすい和書がなかったのですが、

現代疫学: 原著第4版

現代疫学: 原著第4版

  • 学術図書出版社
Amazon
最近出て話題のこの本、今日届いたのですが、生存時間解析のパートわかりやすくて良かったです。
あれ、この記事いらなくないですか？？？

ここまで読んでくださった方で、お時間ある方、
あと、フィードバックのコメントをたくさんお待ちしています。
次はこの記事書いてほしいとかのリクエストも（可能な範囲で）対応します。
コメント、励みになるのでぜひよろしくお願いします。

syleir.hatenablog.com
統計検定対策に役立つ記事のまとめです。他の記事も気合を入れて作っているのでぜひ読んでください。