信頼区間の求め方 信頼区間の求め方でのエラー。

Question

カラム’AAA’と’BBB’で構成した混合行列から得た数値を使ってPositive Likelyhood ratioを求め、その95%CIを求めます。
下記のプログラムで行うと「母比率pの95%信頼区間: nan < p < nan」と数値がnanとなっていしまいます。
どのプログラムにおかしな点があるのか、改善点などおしえていただけませんか？
p<0を期待しているinterval関数では難しいのでしょうか？
代替法があれば教えて頂けるとありがたいです。

from sklearn.metrics import confusion_matrix
import numpy as np
import scipy as sp

y_true = df1['AAA']
y_pred = df1['BBB']
cm = confusion_matrix(y_true, y_pred, labels=[0,1])
tn, fn, fp, tp = cm.flatten()

[tn, fn, fp, tp]
→[1075, 743, 129, 241]

Sensitivity　＝　tp / (tp + fp)
Specificity　＝ tn / (tn + fn)
PLR = sensitivity/(1-Specificity)

alpha = 0.95
n = 2188
p = PLR

bottom, up = sp.stats.binom.interval(alpha=alpha, n=n, p=p, loc=0)
print('母比率pの95%信頼区間: {:.2f} < p < {:.2f}'.format(bottom/n, up/n))
→　母比率pの95%信頼区間: nan < p < nan

Answer 1

spider-manさんが書かれておりますとおり、 は PLR は２郡間での割合の比 であり、分布は二項分布ではありませんので、二項分布より信頼区間を求めてもなんの意味もありません。

PLR の信頼区間の求め方としては、

http://www.snap-tck.com/room04/c01/stat/stat03/stat0304_2.html#note05

https://www.yuzashiki.tokyo/entry/2018/12/16/232333#%E3%83%AA%E3%82%B9%E3%82%AF%E6%AF%94%E3%81%AE%E4%BF%A1%E9%A0%BC%E5%8C%BA%E9%96%93

に書かております 『リスク比の信頼区間』　と同じ考え方で、 PLR の対数（ln(PLR)）をとることで、その分布が分散 (1/TP - 1/(TP+FN) + 1/FP - 1/(TN+FP)) の正規分布 に近似しますので、これを利用することで正規分布から信頼区間を求めることになります。（詳細はリンク先を見てください）

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt


TN, FP, FN, TP = [1075, 743, 129, 241]

sensitivity = TP / (TP + FN)
specificity = TN / (TN + FP)
PLR = sensitivity/(1-specificity)

V = 1/TP - 1/(TP+FN) + 1/FP - 1/(TN+FP)
bottom,up = np.exp(np.array(norm.interval(alpha, 0, np.sqrt(V)))) * PLR
print(f'PLR : {PLR:0.3f} 信頼区間 {bottom:0.3f} < p < {up:0.3f}')
# PLR : 1.594 信頼区間 1.452 < p < 1.749

あと、質問のコードですが、どう見ても FN と FP が逆に定義されているように思えます（私のサンプルでは入れ替えております。）ので確認ください。

Comments

[SyunSyun]

magichanさん。
また助けて頂き有難うございます。
そもそも二項分布より信頼区間を求めていたところから、間違いであったということで、ご指摘有難うございます。
リスク比の信頼区間とおなじく正規分布を利用できるということ理解できました。

Answer 2

二項分布の信頼区間を求める場合、alphaには信頼区間、nには試行回数、pには観測された確率を代入します。上記の例だとpが1を超える数になっているかと思われます。確率として認識されないため、上記のようになっているかと思います。

Comments

[SyunSyun]

有難うございます。
ありがとうございました。
ご指摘の通り、上記の例ではp＞１なので、そこでエラーが出てしまうのかもしれません。
Positive Likelyhood Ratioなどの1を超える値が多いときの95%CIを求めるときは sp.stats.binom.interval（alpah, n , p）は使用できないという事なのでしょうか。
この場合には、使えそうなモジュールなど思いつくものがあったならお教えいただけるとうれしいです。

[spider-man]

> SyunSyunさん
私も統計学の知識不足なので、あまりお力にはなれないかもしれませんが、一応分かる範囲でのみお答えさせて頂きます。

1. stats.binom.intervalについて
こちらは二項分布と呼ばれる確率分布における信頼区間を計算するものです。二項分布の分かりやすい例としてはコイントスが挙げられます。コインの様に表or裏といったが二値の場合において、表が出る確率とコイントスを行う試行回数が分かっているときに、表の出現頻度の信頼区間を求めるのがイメージし易いかと思います。例えば、p=0.3, n=10, alpha=0.95だとすると、表が出る確率が30％のコインを10回投げた時に、表が出る回数は95%でこの間にあるはず。といったイメージです。ですので、nが1より大きいことはそもそもの定義に沿っていないことになります。

2. 使えそうなモジュールということですが、こちらはそもそもの目的とかで変わってきそうなのと、私自身の知識不足もあり、はっきりとお伝えできなさそうです。尤度比検定とかで調べてみると役立つ知識が得られるかも知れません。（曖昧ですいません。）

[SyunSyun]

spider-manさん。
非常に分かりやすい説明有難うございました。
二項分布の信頼区間用のモジュールを使用しているのはそもそも間違いであったのですね。理解できました。
尤度比のための検定をもう少し勉強してみます。
お時間を取って頂き有難うございました。