Rでカッパ係数

# 必要なライブラリを読み込む
library(dplyr)
library(ggplot2)
library(irr) # カッパ係数の計算に利用

# シミュレーションのパラメータ設定
seed <- 20250805
set.seed(seed) # 結果の再現性をための乱数シード
num_patients <- 100 # 患者の数
true_positive_prop <- 0.3 # 真の陽性率
accuracy_A <- 0.95 # 医師A（ベテラン）の診断精度
# 医師B（新人）の診断精度を50%から95%まで5%刻みで変化させる
accuracy_B_levels <- seq(0.50, 0.95, by = 0.05)
num_simulations <- 100 # 各精度レベルでのシミュレーション回数

# シミュレーション結果を格納するためのデータフレームを初期化
simulation_results <- data.frame()

# 医師Bの各精度レベルでシミュレーションを実行
for (acc_B in accuracy_B_levels) {
  # 各精度レベルで複数回シミュレーションを行い、結果を安定させる
  for (i in 1:num_simulations) {
    # 1. 真の状態を生成
    # rbinomを使い、1=陽性, 0=陰性とする
    true_status <- rbinom(num_patients, 1, true_positive_prop)

    # 2. 医師Aの診断をシミュレート
    # 医師Aが間違う確率 (1 - accuracy_A)
    # 医師Aの診断 = 真の状態 XOR 間違い(1) or 正解(0)
    mistake_A <- rbinom(num_patients, 1, 1 - accuracy_A)
    rating_A <- true_status != mistake_A

    # 3. 医師Bの診断をシミュレート
    mistake_B <- rbinom(num_patients, 1, 1 - acc_B)
    rating_B <- true_status != mistake_B

    # 4. カッパ係数を計算
    # irr::kappa2()関数は、評価者を列にしたデータフレームを要求する
    ratings_df <- data.frame(A = rating_A, B = rating_B)
    # 欠損値がない完全なケースのみを使用
    kappa_result <- kappa2(ratings_df, weight = "unweighted")

    # 結果を一時的に保存
    temp_result <- data.frame(
      accuracy_B = acc_B,
      simulation_run = i,
      kappa_value = kappa_result$value
    )

    # シミュレーション結果を結合
    simulation_results <- rbind(simulation_results, temp_result)
  }
}

# シミュレーション結果の要約
# 医師Bの各精度レベルで、カッパ係数の平均値を計算
summary_results <- simulation_results %>%
  group_by(accuracy_B) %>%
  summarise(mean_kappa = mean(kappa_value), .groups = "drop")


# 結果の出力
cat("--- シミュレーション結果 ---\n\n")
print(summary_results)


# ggplot2による結果の可視化
p <- ggplot(summary_results, aes(x = accuracy_B, y = mean_kappa)) +
  geom_line(aes(group = 1), color = "dodgerblue", linewidth = 1.2) +
  geom_point(color = "dodgerblue", size = 3) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "新人医師のスキル向上とカッパ係数の関係",
    subtitle = "ベテラン医師（精度95%）とのレントゲン診断一致度",
    x = "新人医師Bの診断精度",
    y = "カッパ係数（平均値）",
    caption = "各精度レベルで100回のシミュレーションを実施した結果"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 18),
    plot.subtitle = element_text(size = 12),
    axis.title = element_text(face = "bold"),
    plot.caption = element_text(color = "grey40")
  )

# プロットを表示
print(p)