事後確率計算とlogistic regressionによる病気推定アルゴリズムについての考察
suppressMessages(library(tidyverse))
suppressWarnings(suppressMessages(library(epiR)))
suppressMessages(library(caret))1000人分のデータを考える。３種の質問があるとする。
質問A、質問B、質問Cはそれぞれ以下のconfusion matrixを持つ。
質問A
testA <- as.table(matrix(c(30,70,10,90), nrow = 2),
colnames = c("Dis+","Dis-"),
rownames = c("TestA +","TestA -")) %>%
epi.tests(conf.level = 0.95)
testA## Outcome + Outcome - Total
## Test + 30 10 40
## Test - 70 90 160
## Total 100 100 200
##
## Point estimates and 95 % CIs:
## ---------------------------------------------------------
## Apparent prevalence 0.20 (0.15, 0.26)
## True prevalence 0.50 (0.43, 0.57)
## Sensitivity 0.30 (0.21, 0.40)
## Specificity 0.90 (0.82, 0.95)
## Positive predictive value 0.75 (0.59, 0.87)
## Negative predictive value 0.56 (0.48, 0.64)
## Positive likelihood ratio 3.00 (1.55, 5.80)
## Negative likelihood ratio 0.78 (0.67, 0.90)
## ---------------------------------------------------------Positive likelihood ratioが陽性尤度比
Negative likelihood ratioが陰性尤度比
陽性尤度比 = sensitivity / (1 – specificity)
陰性尤度比 = (1 – sensitivity) / specificity
陽性尤度比 3 = ( 30/(30+70) ) / ( 1 – 90/(10+90))
陰性尤度比 0.77 = ( 1 – 30/(30+70) ) / ( 90/(10+60))
質問B
testB <- as.table(matrix(c(60,40,40,60), nrow = 2),
colnames = c("Dis+","Dis-"),
rownames = c("TestB +","TestB -")) %>%
epi.tests(conf.level = 0.95)
testB## Outcome + Outcome - Total
## Test + 60 40 100
## Test - 40 60 100
## Total 100 100 200
##
## Point estimates and 95 % CIs:
## ---------------------------------------------------------
## Apparent prevalence 0.50 (0.43, 0.57)
## True prevalence 0.50 (0.43, 0.57)
## Sensitivity 0.60 (0.50, 0.70)
## Specificity 0.60 (0.50, 0.70)
## Positive predictive value 0.60 (0.50, 0.70)
## Negative predictive value 0.60 (0.50, 0.70)
## Positive likelihood ratio 1.50 (1.12, 2.00)
## Negative likelihood ratio 0.67 (0.50, 0.89)
## ---------------------------------------------------------質問C
testC <-as.table(matrix(c(20,80,30,70), nrow = 2),
colnames =c("Dis+","Dis-"),
rownames = c("TestC +","TestC -")) %>%
epi.tests(conf.level = 0.95)
testC ## Outcome + Outcome - Total
## Test + 20 30 50
## Test - 80 70 150
## Total 100 100 200
##
## Point estimates and 95 % CIs:
## ---------------------------------------------------------
## Apparent prevalence 0.25 (0.19, 0.32)
## True prevalence 0.50 (0.43, 0.57)
## Sensitivity 0.20 (0.13, 0.29)
## Specificity 0.70 (0.60, 0.79)
## Positive predictive value 0.40 (0.26, 0.55)
## Negative predictive value 0.47 (0.38, 0.55)
## Positive likelihood ratio 0.67 (0.41, 1.09)
## Negative likelihood ratio 1.14 (0.97, 1.34)
## ---------------------------------------------------------まずはシュミレーションの基本の関数の動作を確認
２項分布を想定。
rbinom(シュミレーション回数、１回のシュミレーション内での試行回数、確率)
rbinom(10,1,0.7)## [1] 0 0 0 0 1 0 1 1 1 0上の例は、有病率0.7の場合に、ある１人が病気であるかどうかを１か０で返す。 それを１０回繰り返した場合の結果がでてくる。
例えばこれを１００回繰り返した場合
cc = c()
for(i in 1:100){
cc <- append(cc,sum(rbinom(10,1,0.7)))
}
tibble(cc) %>%
ggplot(aes(cc))+
geom_bar(fill="aquamarine3")+
xlab("１０人のうち、何人が感染者であるか")+
ylab("観測回数")+
theme_minimal()
だいたい７人くらいの場合が多くなる。
シュミレーション
1000人のうち、500人が病気・500人が病気ではないとする。
people = 1000
disease_group = people/2
nondisease_group = people/2病気、非病気に対して、それぞれのテーブルを作成。
列名は、病気有無、A陽性、B陽性、C陽性とする。
yesを１、noを１とする。
yes_disease <- data.frame(matrix(ncol = 4, nrow = disease_group))
colnames(yes_disease) <- c("disease","A", "B", "C")
no_disease <- data.frame(matrix(ncol = 4, nrow = nondisease_group))
colnames(no_disease) <- c("disease","A", "B", "C")病気列に、病気であるかないかを入力する。
set.seed(123)
yes_disease$disease <- c(rep(1,disease_group))
set.seed(234)
no_disease$disease <- c(rep(0,nondisease_group))Aの質問に対して、以下の割合をシュミレーションする。
病気のうち、テストAに陽性の割合 0.3 = 30 / (30 + 70)
非病気のうち、テストAに陽性の割合 0.1 = 10 / (10 + 90)
testA$tab## Outcome + Outcome - Total
## Test + 30 10 40
## Test - 70 90 160
## Total 100 100 200set.seed(345)
yes_disease$A <- rbinom(disease_group,1,0.3)
set.seed(456)
no_disease$A <- rbinom(nondisease_group,1,0.1)Bの質問に対して、以下の割合をシュミレーションする。
病気のうち、テストBに陽性の割合 0.6 = 60 / (60 + 40)
非病気のうち、テストBに陽性の割合 0.4 = 40 / (40 + 60)
testB$tab## Outcome + Outcome - Total
## Test + 60 40 100
## Test - 40 60 100
## Total 100 100 200set.seed(567)
yes_disease$B <- rbinom(disease_group,1,0.6)
set.seed(678)
no_disease$B <- rbinom(nondisease_group,1,0.4)Cの質問に対して、以下の割合をシュミレーションする。
病気のうち、テストCに陽性の割合 0.2 = 20 / (20 + 80)
非病気のうち、テストCに陽性の割合 0.3 = 30 / (30 + 70)
testC$tab## Outcome + Outcome - Total
## Test + 20 30 50
## Test - 80 70 150
## Total 100 100 200set.seed(789)
yes_disease$C <- rbinom(disease_group,1,0.2)
set.seed(890)
no_disease$C <- rbinom(disease_group,1,0.3)病気テーブルと非病気テーブルを上下方向に結合する。
df <- rbind(yes_disease,no_disease)シュミレーションがうまくできているか確認する。
paste("病気のうち、A陽性の割合:",
(df %>% filter(disease==1) %>% select(A) %>% sum())/disease_group,
" 病気のうち、B陽性の割合:",
(df %>% filter(disease==1) %>% select(B) %>% sum())/disease_group,
" 病気のうち、C陽性の割合:",
(df %>% filter(disease==1) %>% select(C) %>% sum())/disease_group)## [1] "病気のうち、A陽性の割合: 0.258 病気のうち、B陽性の割合: 0.612 病気のうち、C陽性の割合: 0.202"paste("非病気のうち、A陽性の割合:",
(df %>% filter(disease==0) %>% select(A) %>% sum())/nondisease_group,
" 非病気のうち、B陽性の割合",
(df %>% filter(disease==0) %>% select(B) %>% sum())/nondisease_group,
" 非病気のうち、C陽性の割合",
(df %>% filter(disease==0) %>% select(C) %>% sum())/nondisease_group)## [1] "非病気のうち、A陽性の割合: 0.102 非病気のうち、B陽性の割合 0.376 非病気のうち、C陽性の割合 0.33"事後oddsを計算する。この病気の有病率は50％とする。
事後オッズ ＝ 事前オッズ × 陽性尤度比 or 陰性尤度比
まず有病率50％を事前オッズに変換する。
以下の公式を使う
odds＝P/(1－P)
P＝odds/(1＋odds)
この病気の事前オッズは１となる。
信頼区間を無視して計算を行う。調べたのだが信頼区間を、事後確率に反映させる方法がわからない。信頼区間がある数字を次々に掛け算すると、最終的な事後確率の信頼区間はとても大きくなってしまうと思うのだが。
Aの質問にyesの場合、陽性尤度比 3 をかける
Aの質問にnoの場合、陰性尤度比 0.77 をかける
Bの質問にyesの場合、陽性尤度比 1.5 をかける
Bの質問にnoの場合、陰性尤度比 0.66 をかける
Cの質問にyesの場合、陽性尤度比 0.66 をかける
Cの質問にnoの場合、陰性尤度比 1.44 をかける
計算した確率は、新たな列oddsに保存する。
for(i in 1:people){
odds = 1
if(df$A[i]==1){
odds = 3
}else{
odds = 0.77
}
if(df$B[i]==1){
odds = odds*1.5
}else{
odds = odds*0.66
}
if(df$C[i]==1){
odds = odds*0.66
}else{
odds = odds*1.14
}
df$odds[i] = odds
}病気、非病気ごとの確率を集計する。
オッズから確率に変換する。
df <- df %>% mutate(prob = round(odds/(1+odds), digits = 2 ))
df %>%
group_by(prob,disease) %>%
summarise(counting = n()) %>%
mutate(subsum = sum(counting)) %>%
mutate(percentage = counting/subsum )## `summarise()` regrouping output by 'prob' (override with `.groups` argument)## # A tibble: 14 x 5
## # Groups: prob [7]
## prob disease counting subsum percentage
## <dbl> <dbl> <int> <int> <dbl>
## 1 0.25 0 100 134 0.746
## 2 0.25 1 34 134 0.254
## 3 0.37 0 181 287 0.631
## 4 0.37 1 106 287 0.369
## 5 0.43 0 51 95 0.537
## 6 0.43 1 44 95 0.463
## 7 0.570 0 128 323 0.396
## 8 0.570 1 195 323 0.604
## 9 0.69 0 20 66 0.303
## 10 0.69 1 46 66 0.697
## 11 0.75 0 3 18 0.167
## 12 0.75 1 15 18 0.833
## 13 0.84 0 17 77 0.221
## 14 0.84 1 60 77 0.779グラフ化する
df %>%
group_by(prob,disease) %>%
summarise(counting = n()) %>%
mutate(subsum = sum(counting)) %>%
mutate(percentage = counting/subsum ) %>%
ggplot(aes(x=prob, y=counting,fill=factor(disease)))+
geom_bar(stat="identity",width=0.05)+
geom_text(aes(x=prob,
y=counting,
group = factor(disease),
label = scales::percent(percentage)),
position = position_stack(vjust = .5))+
labs(fill = "Disease")+
xlim(0,1)+
scale_fill_brewer(palette="Set2")+
theme_minimal()## `summarise()` regrouping output by 'prob' (override with `.groups` argument)
logistic regression
今までの事後確率計算は、事前に全ての確率を知っていなければならない（既に学習している）。
一方今度は、機械学習を用いてデータの中にあるパターンを学習する手法を使用する。
説明性が高い機械学習モデルであるlogistic regressionを適応して病気を予測する。
predという新たな列に予測を保存する。
logi <- glm(disease~A+B+C,df,family = "binomial")
df$pred <- predict(logi,newdata = df, type = "response")集計し、グラフ化する。
df %>%
group_by(pred,disease) %>%
summarise(counting = n()) %>%
mutate(subsum = sum(counting)) %>%
mutate(percentage = counting/subsum )%>%
ggplot(aes(x=pred, y=counting,fill=factor(disease)))+
geom_bar(stat="identity")+
geom_text(aes(x = pred,
y=counting,
group = factor(disease),
label = scales::percent(percentage)),
position = position_stack(vjust = .5))+
labs(fill = "Disease")+
xlim(0,1)+
scale_fill_brewer(palette="Set2")+
theme_minimal()## `summarise()` regrouping output by 'pred' (override with `.groups` argument)
事後オッズを用いた場合のスコア計算。
確率が0.5以上の場合を病気と判断する。本来ならROCカーブを用いて、AUCが最大となる数値を模索するのだがベストだが、簡便に行うためとりあえず0.5を使う。
confusionMatrix(
factor( df$disease),
factor( (df$prob > 0.5)*1),
positive = "1"
)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 332 168
## 1 184 316
##
## Accuracy : 0.648
## 95% CI : (0.6175, 0.6776)
## No Information Rate : 0.516
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.296
##
## Mcnemar's Test P-Value : 0.424
##
## Sensitivity : 0.6529
## Specificity : 0.6434
## Pos Pred Value : 0.6320
## Neg Pred Value : 0.6640
## Prevalence : 0.4840
## Detection Rate : 0.3160
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.6482
##
## 'Positive' Class : 1
## logistic regressionを用いた場合のスコア計算
confusionMatrix(
factor( df$disease),
factor( (df$pred > 0.5)*1),
positive = "1"
)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 332 168
## 1 184 316
##
## Accuracy : 0.648
## 95% CI : (0.6175, 0.6776)
## No Information Rate : 0.516
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.296
##
## Mcnemar's Test P-Value : 0.424
##
## Sensitivity : 0.6529
## Specificity : 0.6434
## Pos Pred Value : 0.6320
## Neg Pred Value : 0.6640
## Prevalence : 0.4840
## Detection Rate : 0.3160
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.6482
##
## 'Positive' Class : 1
## どちらの計算も同じAccuracyとなった。
logistic regressionがデータから学んだ各変数の情報を解析する。
summary(logi)##
## Call:
## glm(formula = disease ~ A + B + C, family = "binomial", data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.86892 -0.97030 -0.07237 0.99676 1.65786
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5089 0.1081 -4.707 2.52e-06 ***
## A 1.1138 0.1856 6.002 1.95e-09 ***
## B 0.9498 0.1342 7.079 1.45e-12 ***
## C -0.5737 0.1533 -3.742 0.000183 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1386.3 on 999 degrees of freedom
## Residual deviance: 1274.2 on 996 degrees of freedom
## AIC: 1282.2
##
## Number of Fisher Scoring iterations: 4Coefficientsに、各変数の係数についての情報がある。
Pr(>|z|)を見てわかるように、どの変数も有意に０ではないことがわかる。
各変数のオッズ比を見てみる。
exp(cbind(coef(logi), confint(logi)))## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.6011842 0.4855083 0.7419957
## A 3.0459040 2.1294267 4.4127503
## B 2.5853198 1.9900121 3.3679857
## C 0.5634544 0.4164032 0.7598434Aのよこの数字が3となっている。
これは、Aが０（陰性）のときを基準としたときのAが１（陽性）のオッズ比である。 Aが１の場合は、Aが０のときに比べて、3のオッズ比となる。比が１より高いので、Aが１陽性の場合のほうが病気であるオッズは高いということが分かる。逆に、Cはオッズ比が１未満なのでC陽性であることは病気であるオッズを下げる。
このオッズ比が大きいということは、その変数が確率に対してより大きな影響力を持つ。つまり確率に大きな影響を与える、重要な変数となる。
各係数と各変数の合計が、logodds。
log odds = 0.6 + 3 * A + 2.6 * B – 0.6 * C
log oddsは確率に変換することができる。
確率 = 1 / ( 1 + e^(-log odds) )
例えば、以下の場合で病気になる確率を計算してみる。
テストA陽性、テストB陰性、テストC陽性
logodds <- 0.6 + 3*1 + 2.6*0 - 0.6*1
1 / (1 + exp(-logodds))## [1] 0.95257410.95の確率で、病気であると計算することができます。
ちなみに、logistic regressionでinterceptをなくした場合、得られるオッズ比は陽性尤度比とほぼ同じになる。理由は、intercept内に各変数の０の場合が含まれないため。
logi2 <- glm(disease~A+B+C-1,df,family = "binomial")
exp(cbind(coef(logi2), confint(logi2)))## Waiting for profiling to be done...## 2.5 % 97.5 %
## A 2.3521704 1.6743325 3.3503480
## B 1.7289749 1.4131847 2.1217294
## C 0.4242574 0.3227128 0.5543172先ほどの事後確率の計算：
Aの質問にyesの場合、陽性尤度比 3 をかける
Bの質問にyesの場合、陽性尤度比 1.5 をかける
Cの質問にyesの場合、陽性尤度比 0.66 をかける