Rで重み付けカッパ係数

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

# シミュレーションのパラメータ設定
seed <- 20250806
set.seed(seed) # 結果の再現性をための乱数シード
num_movies <- 100 # 映画の本数
categories <- 1:5 # 評価カテゴリ（1から5）
sd_A <- 0.5 # Aさん（ベテラン）の評価のばらつき（標準偏差）
# Bさん（新人）の評価のばらつきを変化させる
sd_B_levels <- seq(2.0, 0.5, by = -0.1)
num_simulations <- 100 # 各ばらつきレベルでのシミュレーション回数

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

# Bさんの各ばらつきレベルでシミュレーションを実行
for (sd_B in sd_B_levels) {
  for (i in 1:num_simulations) {
    # 1. 真の評価を生成（平均3、標準偏差1.5の正規分布に従うとする）
    true_rating_raw <- rnorm(num_movies, mean = 3, sd = 1.5)
    # 評価を1〜5の整数に丸めてクリッピング
    true_rating <- round(pmax(min(categories), pmin(max(categories), true_rating_raw)))

    # 2. Aさんの評価をシミュレート
    rating_A_raw <- true_rating + rnorm(num_movies, mean = 0, sd = sd_A)
    rating_A <- round(pmax(min(categories), pmin(max(categories), rating_A_raw)))

    # 3. Bさんの評価をシミュレート
    rating_B_raw <- true_rating + rnorm(num_movies, mean = 0, sd = sd_B)
    rating_B <- round(pmax(min(categories), pmin(max(categories), rating_B_raw)))

    # 4. 2種類のカッパ係数を計算
    ratings_df <- data.frame(A = rating_A, B = rating_B)

    # 通常のカッパ係数（重みなし）
    kappa_unweighted <- kappa2(ratings_df, weight = "unweighted")$value

    # 重み付けカッパ係数（2乗重み）
    kappa_weighted <- kappa2(ratings_df, weight = "squared")$value

    # 結果を一時的に保存
    temp_result <- data.frame(
      sd_B = sd_B,
      simulation_run = i,
      unweighted_kappa = kappa_unweighted,
      weighted_kappa = kappa_weighted
    )

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

# シミュレーション結果の要約
summary_results <- simulation_results %>%
  group_by(sd_B) %>%
  summarise(
    mean_unweighted_kappa = mean(unweighted_kappa),
    mean_weighted_kappa = mean(weighted_kappa),
    .groups = "drop"
  )

# プロット用にデータを縦長形式に変換
plot_data <- summary_results %>%
  pivot_longer(
    cols = c("mean_unweighted_kappa", "mean_weighted_kappa"),
    names_to = "kappa_type",
    values_to = "kappa_value"
  ) %>%
  mutate(kappa_type_jp = case_when(
    kappa_type == "mean_unweighted_kappa" ~ "通常カッパ",
    kappa_type == "mean_weighted_kappa" ~ "重み付けカッパ"
  ))


# 結果の出力
print(summary_results)

# ggplot2による結果の可視化
p <- ggplot(plot_data, aes(x = sd_B, y = kappa_value, color = kappa_type_jp)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2.5) +
  scale_x_reverse() + # 横軸を反転させ、右に行くほどスキル向上するように見せる
  scale_color_manual(values = c("通常カッパ" = "darkorange", "重み付けカッパ" = "dodgerblue")) +
  labs(
    title = "新人評論家のスキル向上と2種類のカッパ係数の比較",
    subtitle = "評価の「ズレの大きさ」を考慮する重み付けカッパ係数の効果",
    x = "新人Bさんの評価のばらつき（標準偏差）※小さいほど高スキル",
    y = "カッパ係数（平均値）",
    color = "カッパ係数の種類",
    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"),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    plot.caption = element_text(color = "grey40")
  )

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