ロジスティック回帰(2)

個別データを使用した感度・特異度とROC曲線の求め方

パッケージのインストール

R
# ==== 準備 ====
# 必要パッケージ
#packs <- c("readr", "dplyr", "pROC")
#to_install <- setdiff(packs, rownames(installed.packages()))
#if (length(to_install) > 0) install.packages(to_install, dependencies = TRUE)
library(readr)
library(dplyr)
library(pROC)

pred には、ロジスティック回帰モデル fit2 が推定した「満足 = 1 になる確率(予測確率)」 が挿入される。これはアウトカム(満足)の発生確率をモデル化した値である。なお、ロジスティック回帰で 「暴露(介入)を受ける確率」 を推定した場合、その確率は 傾向スコア(propensity score) と呼ばれる。傾向スコアは因果推論で用いられる概念であり、本分析の「満足 = 1 の確率」とは目的が異なる。

R
dat2 <- dat2 %>%
  mutate(pred = predict(fit2, type = "response"))

感度・特異度を求める
例)カットオフ値を0.5に設定

R
cut_05 <- 0.5
tab_05<- table(
  hat = factor(as.integer(dat2$pred >= cut_05), levels = c(1, 0)),
  obs = factor(dat2$満足, levels = c(1, 0))
)
print(tab_05)


hat: 予測した満足「0と1」、obs: 実際の満足「0と1」

> print(tab_05)
   obs
hat  1  0
  1 35  9
  0 14 42

感度(Sensitivity):実際に満足だった(TP: True Positive, 真の陽性)人のうち、予測でも満足と判定できた割合。


# tab_05 は以下の形(行=予測, 列=真値)を前提:
#           obs=1   obs=0
# hat=1       TP      FP
# hat=0       FN      TN

TP: True Positive
FP: False Positive
FN: False Negative
TN: True Negative
R
TP <- sum(dat2$満足 == 1 & dat2$pred >= cut_05)
FP <- sum(dat2$満足 == 0 & dat2$pred >= cut_05)
FN <- sum(dat2$満足 == 1 & dat2$pred <  cut_05)
TN <- sum(dat2$満足 == 0 & dat2$pred <  cut_05)

sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
cat(sprintf("感度: %.3f\n特異度: %.3f\n", sensitivity, specificity))
> cat(sprintf("感度: %.3f\n特異度: %.3f\n", sensitivity, specificity))
感度: 0.714
特異度: 0.824

ROC曲線とAUC

注意:pROC::roc の response は {0,1}のベクトル、predictor はロジスティック回帰でもとめた連続変数(予測確率)です。

R
roc_obj <- roc(response = dat2$満足, predictor = dat2$pred, direction = "<")
#関数auc(): ROC曲線の下面積(Area Under the Curve = AUC)を計算する
auc_val <- auc(roc_obj)
#aucの95%CI
ci_auc  <- ci.auc(roc_obj)  
cat(sprintf("\nAUC: %.3f (95%%CI %.3f–%.3f)\n", auc_val, ci_auc[1], ci_auc[3]))

Youden index で最適カットオフ

R
best <- coords(
  roc_obj, 
  x = "best", 
  best.method = "youden", 
  ret = c("threshold", "sensitivity", "specificity", "ppv", "npv")
)
cat("\n--- Youden index での最適カットオフ ---\n")
print(best)

# 参考:最適カットオフでの混同行列
cut_best <- as.numeric(best["threshold"])
tab_best <- table(
  obs = dat2$満足,
  hat = as.integer(dat2$pred >= cut_best)
)
cat("\n最適カットオフでの混同行列:\n")
print(tab_best)

ROC曲線の描画

R
plot(roc_obj, print.auc=TRUE, col="blue")

# カットオフ0.5
p05 <- coords(roc_obj, x=0.5, input="threshold",
              ret=c("specificity","sensitivity"))
points(1 - p05["specificity"], p05["sensitivity"], pch=19, col="red")
text(1 - p05["specificity"], p05["sensitivity"] + 0.05,   # ← 上にずらす
     labels="cut=0.5", col="red")

# 最適カットオフ
best <- coords(roc_obj, x="best", best.method="youden",
               ret=c("threshold","specificity","sensitivity"))
points(1 - best["specificity"], best["sensitivity"], pch=19, col="orange")
text(1 - best["specificity"], best["sensitivity"] - 0.05,  # ← 下にずらす
     labels=paste0("best=", round(best["threshold"],3)),
     col="orange")
タイトルとURLをコピーしました