MLB2018の安打数は正規分布を描いているのかQ-Qplotで考察する
今回はRを用いてQ-Qplotに関する話をしたいと思います。Q-Qplotは準一級の範囲だと思います。統計好きの人が問題にしたそうな内容ですよね。知りませんが。
Q-Qplotとは？
q-qplotとは簡単に言うと２つの分布（多くは得られた分布と理想的な分布）の類似の程度を視覚的に把握するplot図の事です。右肩上がりの直線に成る程良いとされています。使う機会を見かけるのは「この得られた分布は正規分布と言っていいのだろうか？」知りたい時とかですかね。
Q-Qplotの実例
今回はQ-Qplotの中でも得られたデータが正規分布を示しているのかを検証します。
こぼ場合はqqnormを使います。
得られたデータが正規分布している時のイメージ
step1 得られた分布を累積密度関数に
標準正規分布(mean=0,sd=1)に従う分布から100個の乱数を抽出し、得られた分布がきちんと正規分布に従っているかを確認しましょう。
まずはN(0,1)に従う分布から100個の乱数を抽出します。
それをx軸に並べ、y軸は累積密度にしたグラフが以下です。今回のサンプルサイズは100なので、1値毎に累積密度は0.01だけ上がっていくので階段状になっています。
負と正が綺麗に50-50に分かれているのが怖いですが抽出された100個の乱数を昇順に並べると以下の通りです。x=-2.36が最も小さい値で1つ目なのでy軸の累積密度関数は0.01の値をとります。x=-2.09では0.02、x=-1.96では0.03…のようにplotしていき、最大値x=2.06で累積密度関数は1をとります。
step2 比較する関数を累積密度関数に
今回は正規分布と比較するので、正規分布の累積密度関数を準備します。
見たことある図だと思います。例えば、標準正規分布表によるとZ=1.96の面積は0.475ですよね。正規分布は対称なのでx=-1.96での累積密度関数は0.025を取っています。
step3 2つの累積密度関数において各分位数（同じ累積密度関数の値)における各xの値をプロット
言葉ではわかりずらいので実際に視覚的に見てみましょう。
左から順にstep1で得た100個の乱数の累積密度分布、step2で得た累積密度分布、そしてq-qplotです。
q-qplotは同じ分位数当たりの値をカウントしていく訳です。例えば40個目の乱数、今回は乱数を100個抽出したのでつまり累積密度関数が0.4の時の各xの値は何でしょう？
左のstep1においてｙ＝0.4と交わる点のxは-0.243です（図では小さくて分かりにくいですが、上図の乱数の昇順で40個目の値と同値なので⁾。
真ん中のstep2でy=0.4と交わる点のx値といえば、標準正規分布表で面積が0.1になるZ値を負にしたものです。標準正規分布表を開くとZ=0.25で面積は0.0987, Z=0.26で面積は0.1026となっています。今回はZ=0.25を採用することにします。
右のQ-Qplotでは縦軸にはstep1で得られた値で今回は-0.24が与えられます。横軸にはstep2で得た値で今回は-0.25が与えられます。この縦軸-0.24と横軸-0.25で与えられる点をplotします。
もう1つやってみましょう。累積密度関数の値=0.05の時を考えます。step1の乱数側では昇順で5番目の値なので-1.62です。step2の正規分布の累積密度分布ではp=0.05でよく使うZ=1.64なので与えられる値は-1.64です。Q-Qplotで左から5番目の点を読み取ると座標は(-1.64,-1.62)です。
これを累積密度関数値=,0.050.4だけでなく0.01から0.02,0.03…0.98,0.99,1.00の100個分でプロットすることでQ-Qplotが完成します。
Q-Qplot for rnorm(100) and normal distribution
Q-Qplotの見方
しつこいようですがQ-Qplotでは同じ分位数(累積密度関数の値)において求められる各数値を取ってplotしたものです。
今回は同じ正規分布から得たので数値のスケールは同じですが、普段は異なる機会が多いでしょう。ですがQ-Qplotでみたいのは2つの分布が類似しているかどうかだけ。重要なのは散らばり具合の類似度な訳です。同じ分位毎にplotしているので、どちらの値も尺度は違えど全体のスケールからみて同程度の値を取っていればy=xに近づきそうですよね。またy=xからの解離・ゆがみ具合でどのように分布が歪んでいるかを把握することが出来ますし、問題として作り易そうです。
（暇な時に問題作りたい)
MLB選手の安打数は正規分布を示しているかをQ-Qplotで検証
2018年に出場したMLB選手の安打数は正規分布を示しているのかを検証したいと思います(してないと思うので失敗例として扱う事を前提にしています)。本当はNPBにしようとしたけどMLBの方が人数多いしデータ揃ってそうなので。ここから引用しました。選手の抽出や.CSVでdataを得られたり本当に便利。すごい。
step1 2018年のMLB選手の安打数を累積密度分布に落とす
まずはヒストグラムを描くと
0~20本が多くて明らかに正規分布ではないです。これを累積密度分布に落とし込むと
0が多すぎて初っ端から累積密度が0.4位から始まるという笑
どうしてこうなるかというと毎日レギュラーで活躍している人は一握りで、4割近くの人はマイナーを往ったり来たりする控え選手or投手だからでしょう。
これでいいのかと考えた私はヒット数が20本以下は切り捨てて、21本以上の人に絞り込みました（本当は規定打席の方がいいんでしょうが)。すると
正規分布ではないですが、まだ使えそうな分布が得られました。累積密度分布に変換
すると
y=√x ,log(x)のような関数になりました。せっかくなのでこの2つをそれぞれQ-Qplotで正規分布と比較していきます。
step2 正規分布を累積密度分布に
これは前のと同じなので省略
step3 2つの累積密度関数において各分位数（同じ累積密度関数の値)における各xの値をプロット
まずはヒット数で限局していない1つ目の累積密度分布と正規分布でかけてやると
直線とは程遠い曲線が得られました。3つ並べると
左がヒット数の累積密度分布、真ん中が標準正規分布の累積密度分布、右がQ-Qplot.横軸が正規分布のx値で縦軸がMLB側のx値つまりヒット数です。初めは横軸が進んでも縦軸が全く縦に伸びていないのは0付近の値が多いからです。つまり半数近くは0付近のヒット数を誇っているのでしょう。
次に、ヒット数21本以上に限局した場合のQ-Qplotでは
直線とは言い難いですが、1つ目の限局しない場合よりはマシだと思います。それでも極端でないとはいえ本数が少ない選手の方が相対的に多くなっているので最初は鈍く、途中から鋭く曲がるような曲線になってしまっています。
よってMLBの選手の安打数は正規分布していないことが分かります（この程度だとhistで分かりますが）。
Q-Qplot as for the number of hit at(BAT) 2018MLB
Q-Qplotの正確な理解(未)
上で説明したのが１番イメージしやすい考えだと思っているし、Q-Qplotってそんなに正確な情報を得るために使う物ではないと勝手に思っているので十分だと思うのですが、実際は違うみたいなんですよね～ 頭いい人なら余裕で理解できるのでしょうが
だってこの考えだと累積密度関数の値が1の時の正規分布側の値って無限になりません？ ずっと引っかかってました。wikipediaを拝見した所、
正規確率プロットを使用する際には、標準正規分布の順序統計量の期待値の尺度であるランキットを使用する。
なんとなくわかった。でも確認のために標準正規分布の順位統計量を自力で求めてみればいいのだろうけどちょっと面倒くさい…それにRでhelp(qqnorm)を入れても全部英語だし。上で説明したのも的外れ？ 今後も勉学に努めます…
青木宣親から学ぶ 野球は統計学で説明できるのか？~2項分布編~
みなさんは青木宣親選手をご存知でしょうか？
Noriこと青木宣親選手は日本だけでなくMLBでも活躍できた日本の野球史に名を残す現役の名選手です。
青木選手は1982年に宮崎県日向市で産声を上げました。
それからというものの、早稲田時代ではリーグ4連覇を果たし、ヤクルトに入団後も順調に成績を重ね首位打者や最多安打・盗塁王といったタイトルを総なめにし、アナウンサーと結婚。そしてMLBに移籍しました。福留孝介選手や松井稼頭央選手など多くの名選手はMLBの成績が不安定で大幅に落としていた中で青木選手はMLBでも着実に成績を残していきました。その成績は
2012 151試合 打率.288
2013 155試合 打率.286
2014 132試合 打率.285
2015 93試合 打率.287
2016 118試合 打率.283
2017 110試合 打率.277
通算打率.285
ヒトは青木をこう呼びます
ミスター2割8分
2割8分に収束する男
アヘアへ0.280マン
青木宣親選手はシーズン中に波はありました。
今年は無理だろうと巷で交わされていた時も数知れず。
なのに終わってみると必ず2割8分台に乗せているのです。
こんなことできるでしょうか？
多くの選手は毎年少なからず成績が変動します。成績は体調面だけでなく私生活やモチベーション、相手チームの分析など多岐に依存するからでしょう。
この異端な成績を統計学の側面から考えるとどのように評価できるのでしょうか？というのが今回の議題です。
野球と統計学
ここに野球と統計の歴史を書こうと思ったのですが（マネーボールとか）、まあ要らないかなと思うので
工事中(着工日未定)
打率を統計学で考える
今回は青木選手が毎年同じような打率で安打を重ねる事は統計学的にどういう意味を持つのかを考えます。ですが複雑に考えると大変になるので状況は簡潔にします。青木選手はどのような投手・状況・場面でも0.285の確率(打率)で安打を打つと仮定します。野球が好きな人ほど抵抗がある理論だと思いますが、打数が多くなれば案外イイ感じになるんです。
この仮定において今回は、1年間の打率を2項分布で考えることが出来ます。
この青木選手が1年間で1打席だけ立ったシーズンを1000回繰り返した(1000年間シーズンを過ごす)場合、1000シーズン分の打率はどうなるでしょうか？
このplotだけでは正確な数字は分かりませんが、1000×0.285=285年は打率が10割で715年は打率が0割と言う期待値に近いことは読み取れます。
※ table(a)で正確な数字を求めたところそれぞれ740と260でした。
では1シーズンで10打席に立った青木選手の年間打率を求める事象を1000回繰り返した場合、各シーズンの打率の値はどうなるでしょうか？
| 打率 | シーズン数 |
|---|---|
| 0 | 34 |
| 0.1 | 142 |
| 0.2 | 280 |
| 0.3 | 242 |
| 0.4 | 175 |
| 0.5 | 96 |
| 0.6 | 26 |
| 0.7 | 2 |
| 0.8 | 3 |
こんな感じです。
aoki_10bats added way to make table for html style ...
10打席しかないので打率も0から1まで0.1刻みですが流石に9割打者と10割打者は居ませんね。
次は、青木選手が1シーズンで100打席に立ったシーズンの打率を求める試行を1000回繰り返す(1000シーズン分)った時の分布は
| 打率 | 0.16 | 0.17 | 0.19 | 0.2 | 0.21 | 0.22 | 0.23 | 0.24 | 0.25 | 0.26 | 0.27 | 0.28 | 0.29 | 0.3 | 0.31 | 0.32 | 0.33 | 0.34 | 0.35 | 0.36 | 0.37 | 0.38 | 0.39 | 0.4 | 0.41 | 0.42 | 0.43 | 0.44 | 0.46 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| シーズン数 | 1 | 2 | 2 | 8 | 11 | 28 | 31 | 34 | 45 | 62 | 79 | 80 | 70 | 83 | 85 | 72 | 68 | 53 | 45 | 41 | 27 | 19 | 14 | 17 | 14 | 2 | 4 | 2 | 1 |
打率0.31が85シーズンで最頻値となっています。一応、中心極限定理で正規分布を示していますが。
aoki_100bats added ways to change rows and comns
では今回のメインです。1年間怪我無く無事に送った場合に立つ打席数は500打席程度でしょう。ですので、青木選手が1シーズンで500打席に立ったシーズンの打率を求める試行を1000回繰り返す(1000シーズン送)った時の分布は
| 0.22 | 0.222 | 0.226 | 0.23 | 0.232 | 0.234 | 0.236 | 0.238 | 0.24 | 0.242 | 0.244 | 0.246 | 0.248 | 0.25 | 0.252 | 0.254 | 0.256 | 0.258 | 0.26 | 0.262 | 0.264 | 0.266 | 0.268 | 0.27 | 0.272 | 0.274 | 0.276 | 0.278 | 0.28 | 0.282 | 0.284 | 0.286 | 0.288 | 0.29 | 0.292 | 0.294 | 0.296 | 0.298 | 0.3 | 0.302 | 0.304 | 0.306 | 0.308 | 0.31 | 0.312 | 0.314 | 0.316 | 0.318 | 0.32 | 0.322 | 0.324 | 0.326 | 0.328 | 0.33 | 0.332 | 0.334 | 0.336 | 0.34 | 0.348 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 3 | 3 | 2 | 3 | 3 | 2 | 7 | 7 | 7 | 4 | 9 | 12 | 4 | 19 | 17 | 19 | 15 | 22 | 31 | 34 | 29 | 38 | 42 | 30 | 56 | 43 | 44 | 43 | 35 | 34 | 35 | 40 | 31 | 29 | 37 | 21 | 26 | 18 | 20 | 19 | 15 | 16 | 13 | 11 | 10 | 6 | 5 | 10 | 8 | 2 | 1 | 2 | 2 | 1 | 1 |
期待値となる打率の0.282 , 0.284になるシーズン数はそれぞれ43,44シーズンで、2割8分台に収束しているシーズン数は221シーズンでした。
つまり本来の打率が0.285の選手でもせいぜい500打席程度では選手の真の能力を測るには不十分であり、ぶれが生じるという話です。試しに両側95％信頼区間を計算してみると0.244≦打率≦0.326 であります。このように500打席程度では本来の能力とはかけ離れた数値を叩き出すことも容易に起こるのです。
この結果から打率を2項分布で考えるのがナンセンスだ、確率論は野球に応用できないという旨の反論も当然あると思います。
が
一方で、特定の季節・相手投手・球場 etcで調子が悪いというのも実はこのような2項分布の偶然のぶれに起因するものであるという考えもできます。
最終年だけ活躍するおかわり君や夏に強いゴジラも偶然に起因した結果なのかもしれません。また、2軍の数十打席の数字だけで1軍に上げるかどうかを判断する事にも疑問符が呈されるでしょう。
というように、この打率の考え方では1シーズン当たりの打率にぶれは欠かさないので、より選手の能力を測れる変数を探すことになりました。その結果OPSなどが作られて現在ではそれらにスポットライトが当てられている訳なのですが、青木選手はこの自然の摂理の予想に反した結果を残している訳です。
500打席あれば本来の力の打率付近に収束する確率は高いし今の言動は大袈裟かもしれませんが、少なくとも青木選手のアヘアへ0.280マンは2項分布の野球への応用に改善の余地がある可能性を示すものとなっているでしょう。
※まだ不十分なのは承知なので、今後も暇な時に野球と統計の話を追加していきます。
2項分布の話
2項分布は数学か賭博が好きな方ならお馴染みの分布でしょう。気長に聞いてください。
ベルヌーイ分布
2項分布について説明する為に、まずベルヌーイ分布(試行)について知っておくと理解が円滑に進むでしょう。ベルヌーイ分布とは成功と失敗の2項から成る分布で、仮に成功する確率をp,失敗する確率を1-pとし、成功時の確率変数Ｘは１、失敗時の確率変数Xは０とすると、この分布の確率密度関数は
となります。仰々しいですが、これはただ成功（X=1)は確率pで、失敗(X=0)は確率1-pで起こると言っている事を数式にしただけです。これをベルヌーイ分布(試行)といいます。
期待値と分散はそれぞれ
で与えられますね。
2項分布
ではこのベルヌーイ試行をn回繰り返したとき、成功した回数(確率変数=Yとする)はどのようになるでしょうか？一般化したY=k(k=0,1,2,,,,,,n)を言い換えると、n回の施行のうちk回は成功してn-k回は失敗したときの確率を求めればよい訳です。これはつまり
で与えられます。次に期待値と分散を求めたいのですが、確率密度関数から直接求めるのは大変そうです。そこでベルヌーイ試行の話を持ってきます。
ベルヌーイ分布の確率変数Xは成功時に1,失敗時に0をとります。従ってn回試行して計で成功した回数Yは
で与えられます。Xの各試行が互いに独立だとすると期待値と分散は加法性により
となります。
中心極限定理
あ、要らないなって気づいたのも後の祭りでグラフを作ってしまったので一応書きます。2項分布もnが大きければ中心極限定理によって正規分布に近似できるって話です。
成功率pは0.3にしました。
まずはベルヌーイ試行を1回行って成功した回数を確率変数とし、この試行を1000回繰り返すと
となります。1000＊0.3＝300 なので成功した回数(a=1)は約300回になりそうだという推測通りです。
次にベルヌーイ試行を10回繰り返して成功した回数を求める試行を1000回行うと
こんな感じ。やっぱり右に長くなりがち。
次はベルヌーイ試行を100回繰り返して成功した回数を求める試行を1000回行うと
綺麗。
ベルヌーイ試行を1000回繰り返して成功した回数を求める試行を1000回行うと
ベルヌーイ試行を10000回繰り返して成功した回数を求める試行を1000回やって
正規分布っぽい。実際に理論的に導き出される正規分布(期待値3000,分散2100)を重ねると
まあなんということでしょう。初めはあんなに質素だったヒストグラムが匠の技により正規分布に従う様相を呈しています。
Rのコードは以下の通りでした。
正規分布と中心極限定理について
to 正規分布は統計を学んでいく上で最も重要な連続分布であると結論付けられるでしょう。なぜならば、中心極限定理という定理によって正規分布以外の分布でも標本数nが大きければ標本の平均値の分布は正規分布に近似できるという側面を正規分布は抱えているからです。今回はそんな正規分布について考えていきたいと思います。
正規分布ってどんなもの？
みなさんは数学の授業で1次関数や2次関数を習いませんでしたか？1次関数は直線で、2次関数は放物線を描きます。傾きや曲がり度は係数に依存しています。正規分布はどのような形を描くのでしょうか？正規分布は左右対称で山なりの曲線を描きます。また、1次関数や2次関数の形は係数に依存するように、正規分布の形は期待値と分散に依存します。傾きがp,切片がqの1次関数は
で与えられます。同様に、期待値がµ,分散がσ²の正規分布の確率密度関数は
で与えられます。うえっ。なので一般にはN(µ,σ²)で省略します。正規は英語でnormalだからNなのでしょう。
上のグラフはN(0,1),N(2,1),N(0,2)の正規分布関数をRでグラフにしたものです。
ソースコードはこんな感じ
N(0,1),N(2,1)を比較するとN(2,1)はx軸に2だけ動かしたものとなっています。正規分布の確率密度関数をみてみるとxとµは(x-µ)²で与えられている事から納得でしょう。
N(0,1),N(0,2)の違いは分散の大きさです。期待値は等しいので対称軸は同じになっています。分散と言うのは分布の散らばり具合なので、大きければ平たくなります。確率密度関数なので面積は１になることから必然的に背丈も低くなることが予想できるでしょう。
サイコロで考える正規分布と中心極限定理
これだとただの数学なので実際の例で考えていきます。
サイコロを投げる試行を何度か繰り返し、その平均を求めます。
まずはサイコロを1回投げて出た目を数える施行を1000回繰り返すと
このようになります。左から1，2，3，4，5，6が出た回数です。どれも160回前後で均一になっています。これを正規分布だとは政治家でも言えないでしょう。
コードはこんな感じ
次に、サイコロを10回振って出た目の和の平均を数える試行を1000回繰り返すと
となります。少し山なりになりましたが、まだ凸凹しています。
こちらのコードは
次に、サイコロを100回投げて出た目の平均を数える試行を1000回繰り返すと
が得られました。許容範囲内の綺麗な山なりじゃないでしょうか？
コードは
この調子でサイコロを1000回投げて出た目の平均を数える試行を1000回繰り返すと
いいですね～
コードは
最後にサイコロを1万回投げて出た目の平均を数える試行を1000回行うと
ヤバイ、美しE超えて美しFやんけ!
このように、標本サイズnが大きくなるほど正規分布に近づくのです！これが中心極限定理です。
実際に、サイコロを1000回投げて出た目の平均を1000回繰り返したヒストグラムに 期待値 3.5 , 分散 35/120000 (何故こうなるかはそのうち) の正規分布を重ね合わせると
ピッタリ当てはまりますね！
1000dice mean histgram added dnorm curve
このように、サイコロの目がそれぞれ出る確率は1/6で等しいのに標本サイズnが大きくなるほど正規分布に近似できることがわかります。
ですので、正規分布はよくみかける＆使いやすい分布なのです。
研究デザイン(コホート研究・症例対照研究・横断研究)
研究について漠然と知っている方も多いと思いますが、研究の種類は多岐に渡ります。今回はその中でも最も有名な、時間軸に沿って分ける事で区別されるコホート研究・症例対照研究・横断研究の3つについて説明していきたいと思います。
コホート研究
コホート研究では、現在集まっている被験者を２つのグループに分けます。そこで、片方には要因があり、もう片方には要因がないようにします。そして未来まで研究を継続させ、疾患の有無を各グループで比較する研究スタイルのことです。
具体例で考えます。あなたはコーヒーを飲んだら寿命は延びるのか？について研究する事にしました。まずあなたは今いる100人の被験者を2つのグループに均等に分けます。片方のグループには死ぬまで毎日コーヒー1杯を飲んでもらいます。もう一方のグループにはコーヒーを一生飲んでもらわないようにします。そして数十年後あなたは被験者が全員死んだ時に、コーヒーを飲んでいたグループと飲んでいなかったグループで寿命の長さを比較します。コーヒーを飲んでいた側は平均100歳まで生きたけど飲んでいない側は65歳だった。きっとコーヒーはいい！ と結論付けるのがコホート研究です。この場合、要因はコーヒーを飲むことで疾患は死亡です。
この研究スタイルは3つの中で最も良いデザインだとされています。何故ならコホート研究では現在から未来に向かう前向きの時間軸で時間に沿って研究を行うことができるからです。自然な流れですよね。従って因果関係が明確です。加えてバイアス、つまり偏りが少ないです。
しかしデメリットとしてお金と時間がかかります。具体例で挙げたコーヒーの場合だと、被験者が死ぬまでずっとコーヒーを与えなければいけません。慣れてくると既存の安いインスタントの味に我慢できずに「オーガニック」「イタリア製のドリップは違う」「急冷させて香りを閉じ込める」「上 島 珈 琲」なんて言われるかもしれませんね。お金と手がかかります。加えて時間がかかります。死ぬまで被験者を追いかけるので何十年かかるんでしょう、途方にくれそうです。
症例対照研究（ケースコントロール研究）
症例対照研究（ケースコントロール研究）はコホート研究とは逆で、現時点で被験者を疾患があるグループと無いグループに分けて、それぞれ過去に要因があったかどうかを過去に遡って調べる方法です。
これも具体例で考えましょう。60代で亡くなったグループと60代でまだピンピンのグループに分け、過去に酒豪であったかを聞き取ります(亡くなった場合は親族に聞き取りましょう)。60代で亡くなった方の方が酒豪の割合が高ければ飲酒と寿命は関係してそうだな、なんて察しが付きます。これが症例対照研究です。
症例対照研究はお金も時間もかかりません。それぞれのグループの被験者に「過去に酒を浴びるほど飲みましたか？」と聞くだけです。それが最大の魅力でしょう。
一方デメリットとしてバイアスが大きい、つまり信頼性・エビデンスが低いので研究者や論文の読者はそれを踏まえて処理して理解する必要があります。過去に酒をどのくらい飲みましたか？と聞いても正確な情報はなかなか得られないでしょう。
横断研究
横断研究とはこれまでの２つとは異なり、時間の流れは有りません。現在の瞬間での要因と疾患の有無を調べる方法です。
具体的に考えます。最近、歯周病は全身の疾患と関係していることが分かっています。それを踏まえてあなたは歯周病と心疾患の関係を研究することにしました。歯周病を要因、心疾患を疾患として、被験者に現在歯周病と心疾患がそれぞれあるかどうかを調べ、各グループの人数を調べることで関係が有るのか考察します。時間軸を横断しているのでこれを横断研究と呼びます 。
これもお金や時間はかかりませんよね。患者に「歯周病か心疾患を罹患してます？」って聞けばいいだけです。
一方でこのデザインでは時間軸に沿わず、瞬間のデータしか参考にしていません。ということはもし歯周病と心疾患に関係がありそうでも、歯周病によって心疾患が引き起こされるのか、心疾患によって歯周病が引き起こされるのか、はたまた双方的なのかという因果関係が分かりません。その点を留意する必要があります。
ROC曲線とは？
前回の感度と特異度の話を発展させていきます。
使う表は引き続き10人の癌マーカーの数値を使います。
tmduto.hatenablog.com
| 罹患or非罹患 | マーカーA濃度 |
|---|---|
| 罹患者 | 0.19 |
| 罹患者 | 0.17 |
| 非罹患者 | 0.17 |
| 罹患者 | 0.13 |
| 非罹患者 | 0.12 |
| 罹患者 | 0.09 |
| 非罹患者 | 0.07 |
| 非罹患者 | 0.06 |
| 非罹患者 | 0.04 |
| 非罹患者 | 0.01 |
| 診断 | 合計 | |||
|---|---|---|---|---|
| 陽性 | 陰性 | |||
| 癌 | 有(罹患) | 3(人) | 1(人) | 4(人) |
| 無(非罹患) | 2(人) | 4(人) | 6(人) | |
| 合計 | 5(人) | 5(人) | 10(人) |
前回、腫瘍マーカーの濃度は個人差が有るので誤診が生じると言いました。
また、診断は基準値より上か下かで陽性か陰性か判断しますよね。
ではその基準値(カットオフ値)を変えれば患者の診断も変わると推測できます。
基準値(カットオフ値)を上げれば陽性は増え、陰性は減りそうだし、
基準値(カットオフ値)を下げれば陽性は減り、陰性は減りそうです。
具体的に基準値(カットオフ値)を上げ下げしてみましょう。
カットオフ値を変化させると…？
基準値(カットオフ値)を0.08にしてみます。
| 罹患or非罹患 | マーカーA濃度 |
|---|---|
| 罹患 | 0.19 |
| 罹患 | 0.17 |
| 非罹患 | 0.17 |
| 罹患 | 0.13 |
| 非罹患 | 0.12 |
| 罹患 | 0.09 |
| 非罹患 | 0.07 |
| 非罹患 | 0.06 |
| 非罹患 | 0.04 |
| 非罹患 | 0.01 |
わかりづらいので2×2の表にすると
| 診断 | 合計 | |||
|---|---|---|---|---|
| 陽性 | 陰性 | |||
| 癌 | 有(罹患) | 4(人) | 0(人) | 4(人) |
| 無(非罹患) | 2(人) | 4(人) | 6(人) | |
| 合計 | 6(人) | 4(人) | 10(人) |
感度が1になりましたね！
特異度は4/6=0.6666で変化していません。よって特異度を変えずに感度が高くなったので、この検査は精度が高くなったといえるでしょう。
では逆に基準値(カットオフ値)を0.13にあげてみましょう。
再び2×2の表にすると
| 診断 | 合計 | |||
|---|---|---|---|---|
| 陽性 | 陰性 | |||
| 癌 | 有(罹患) | 3(人) | 1(人) | 4(人) |
| 無(非罹患) | 1(人) | 5(人) | 6(人) | |
| 合計 | 4(人) | 6(人) | 10(人) |
こちらも特異度が5/6=.83333で改善しました！
感度は3/4=0.75で変化してないので、こちらも検査の精度が改善したといってよいでしょう。
とはいえ、極端にとっても良くないだろうと察しますよね。
たとえばカットオフ値を0.2にすると全員陰性になります。確かに非罹患者は全員陰性になるので特異度は1になりますが、癌罹患者も全員陰性になるの感度は0です。
このように基準値(カットオフ値)の選び方が重要となっていく訳です。
ROC曲線とは？
では基準値(カットオフ値)をどこに設定すればいいのでしょうか？
そもそも上のように感度と特異度は表裏一対です。感度・特異度を両方高める事はできません。つまり、感度と特異度の両方が適度な値をとれる場所を知りたいです。では、基準値(カットオフ値)を変数として感度と特異度の値の動きが可視化できればいいなあと考えます。それがROC曲線です。
ROC曲線では横軸に偽陽性率(1ー特異度)、縦軸に陽性率(感度)をとります。そして、基準値(カットオフ値)を変えて生じる幾分かの感度と特異度のセットの点をプロットして、最近傍の点同士を結んだ線がROC曲線です。
ROC曲線を描くのは面倒なのでRを使いましょう。
莫大なerrorやwarningがでてきましたがグラフは得られたので良しとしましょう。
こんな仕上がりです。
カットオフ値が 0.1 , 0.08 , 0.1 の時の感度と特異度はそれぞれ
(1-特異度,感度)=(0.3333 , 0.75) , (0.3333 , 1 ) , (0.16666 , 0.75)
でした。御覧のように、これらはROC曲線の点(角)になっていますよね。
このように、カットオフ値を変化させた場合の(1-特異度,感度)の各点を繋いだ線がROC曲線です。
ROC曲線の目的は、一目で感度と特異度の変移がみてとれることです。
感度と特異度が互いにいい点、つまり左上に近い点を取る時のカットオフ値ほど採用すべき値だと考えられるでしょう。
感度・特異度とかの話(仮)
現在の日本の医療は西洋医学を主軸としている訳ですが、東洋医学と比較して幾つか特徴があります。例えば西洋医学では風邪を引いた患者にずっと同じ風邪薬が処方されるのに対して、東洋医学では風邪の引き始めやピーク・引き終わりは異なる状態と考えて別の薬を処方されます。その中でも西洋医学の１番の特徴と言えば検査医学でしょう。検査数値を病気の診断に利用します。非常に便利なのですがどうしても誤診は避けられません。永久機関は理論上不可能である事と似てます。
誤診とは？
様々な種類のある癌腫瘍の中には、癌細胞に特異的な物質を放出する腫瘍が有ります。この物質は血液中に溶け込みます。つまり、採血検査でその物質の濃度が異常値であれば癌であると診断できます。この物質はマーカーとして働くことから腫瘍マーカーと呼ばれます。
例えばある腫瘍マーカー（マーカーAとします)の血液濃度が0.1mg/ml 以上なら癌と診断する基準を考えます。ここで癌患者と非癌患者合わせて10人の患者に検査を行った時をみてみましょう。
| 罹患or非罹患 | マーカーA濃度 |
|---|---|
| 罹患者 | 0.19 |
| 罹患者 | 0.17 |
| 非罹患者 | 0.17 |
| 罹患者 | 0.13 |
| 非罹患者 | 0.12 |
| 罹患者 | 0.09 |
| 非罹患者 | 0.07 |
| 非罹患者 | 0.06 |
| 非罹患者 | 0.04 |
| 非罹患者 | 0.01 |
前置きとして重要なのは、人には個人差があるので癌患者なのにマーカーA濃度が低かったり、癌患者ではないのにマーカーA濃度がなぜか高い人がどうしても存在します。上の10人の例でも、上から3人目の患者は癌ではないのにマーカーA濃度が0.1mg/ml以上なので、この検査では癌を罹患していると誤診されてしまうのです。
感度・特異度とは？
上の表ではわかりづらいので、2×2のテーブル表にしましょう。
| 診断 | 合計 | |||
|---|---|---|---|---|
| 陽性 | 陰性 | |||
| 癌 | 有(罹患) | 3(人) | 1(人) | 4(人) |
| 無(非罹患) | 2(人) | 4(人) | 6(人) | |
| 合計 | 5(人) | 5(人) | 10(人) |
がんと診断する基準(カットオフ値といいます)を0.1mg/mlにした場合、表１はこのようにまとめられます。
この診断が使えるかどうかを調べる尺度として最も最適な指標は何でしょう？やはり検出力が知りたいですよね。癌患者の何％がきちんと癌と診断されるのかと非癌患者の何％がきちんと陰性として診断されるのかが気になります。これが感度と特異度です。
| 診断 | |||
|---|---|---|---|
| 陽性 | 陰性 | ||
| 癌 | 有(罹患) | 真陽性 | 偽陰性 |
| 無(非罹患) | 偽陽性 | 真陰性 |
感度 =真陽性/(真陽性+偽陰性)
| 感度 | = | 真陽性 |
| (真陽性+偽陰性) |
特異度=真陰性/(真陰性+偽陽性)
| 特異度 | = | 真陰性 |
| (真陰性+偽陽性) |
今回の例で考えてみましょう。感度と特異度はそれぞれ
感度 ：3/4=0.75
特異度：4/6=0.6666
となります。
つまりこの検査では
- 癌患者の75％は陽性と診断される =25%は陰性と誤診される
- 非癌患者の67％は陰性と診断される=23%は陽性と誤診される
だろうと考えられます。これをどう解釈するかは研究者に一任されますが、さすがにこの検査はよくないんじゃないかと考えられます。患者10連ガチャ引いても2.5人は見過ごされますからね。
陽性的中度・陰性的中度とは？
感度と尤度比は説明しました。でも私達は自分が癌であるかどうかなんて知らないし、それを知るために検査を受けています。どちらかというと、陽性と診断された際にどの程度の確率で本当に癌なのか、又は陰性と診断された際にどの程度の確率で本当に癌で無いのかを知りたいと思います。それが陽性的中度と陰性的中度です。
ですので
陽性的中度=真陽性/(真陽性+偽陽性)
陰性的中度=真陰性/(真陰性+偽陰性)
で示されます。今回の例で考えると
陽性的中度=3/5=0.6(60%)
陰性的中度=4/5=0.8(80%)
となります。
もしあなたがこの検査を受けると…
- 陽性がでた場合、60%の確率で本当に癌を患っている=40%は実際に癌でない誤診
- 陰性がでた場合、80%の確率で本当に癌を患っている=20%は実際に癌である誤診
ということになるでしょう。
これだけでみると使えませんね笑
※(仮)というのは、付け加えたい図が頭にあったり表に色を付けたかったり数式(分数)を正式な形で表記したい思いが有るのですが、技術力不足で叶えられていないので今後ブログ力をつけたら適時改良編集していきたいと思います。
