Rの関数qmvt {mvtnorm}を利用したガブリエル比較区間

# 3. クリティカル・バリューの計算 (qmvtの利用)

# (1) 通常のt分布のクリティカル値 (95%両側)
# Bonferroni補正なしの単純なCI用
qt_crit <- qt(0.975, df = total_df)

# (2) ガブリエル区間用のクリティカル値 (SMM分布)
# SMM分布は、相関行列が単位行列である多変量t分布と等価です。
# tail = "both.tails" で絶対値の最大値を考えます。
corr_mat <- diag(k) # 単位行列 (グループ間は独立とみなす)

qmvt_res <- qmvt(
  p = 0.95,
  df = total_df,
  corr = corr_mat,
  tail = "both.tails",
  seed = seed
)

qsmm_crit <- qmvt_res$quantile
cat("--- クリティカル値 ---\n")
cat(sprintf("通常のt値 (95%%): %.4f\n", qt_crit))
cat(sprintf("SMM値 (Gabriel用): %.4f\n\n", qsmm_crit))


# 4. 区間幅の計算
# ガブリエル区間の半径 = (SMM値 * SE) / sqrt(2)
# sqrt(2) で割ることで、区間の重なり判定が検定と一致するように調整される

plot_df <- stats_df %>%
  mutate(
    # 通常の95%信頼区間
    CI_Lower = Mean - qt_crit * SE,
    CI_Upper = Mean + qt_crit * SE,

    # ガブリエル比較区間
    Gab_Width = (qsmm_crit * SE) / sqrt(2),
    Gab_Lower = Mean - Gab_Width,
    Gab_Upper = Mean + Gab_Width
  )

cat("--- 算出された区間幅 ---\n")
print(plot_df %>% select(Group, CI_Lower, CI_Upper, Gab_Lower, Gab_Upper))


# 5. 可視化 (エラーバーの重なり比較)
# 見やすくするために、データを縦持ちに変換して描画
plot_data_long <- bind_rows(
  plot_df %>% mutate(Type = "通常の95%信頼区間", Lower = CI_Lower, Upper = CI_Upper),
  plot_df %>% mutate(Type = "ガブリエル比較区間", Lower = Gab_Lower, Upper = Gab_Upper)
)

p1 <- ggplot(plot_data_long, aes(x = Group, y = Mean, color = Type)) +
  # エラーバー
  geom_errorbar(aes(ymin = Lower, ymax = Upper),
    width = 0.2, size = 1.2, position = position_dodge(width = 0.5)
  ) +
  # 平均値の点
  geom_point(size = 4, position = position_dodge(width = 0.5)) +

  # AとBの境界線（比較用）
  geom_segment(
    aes(
      x = 0.8, xend = 2.2,
      y = plot_df$CI_Upper[plot_df$Group == "A"],
      yend = plot_df$CI_Upper[plot_df$Group == "A"]
    ),
    linetype = "dotted", color = "blue", alpha = 0.5
  ) +
  geom_segment(
    aes(
      x = 0.8, xend = 2.2,
      y = plot_df$Gab_Upper[plot_df$Group == "A"],
      yend = plot_df$Gab_Upper[plot_df$Group == "A"]
    ),
    linetype = "dotted", color = "red", alpha = 0.5
  ) +
  labs(
    title = "通常の信頼区間 vs ガブリエル比較区間",
    y = "値",
    color = "区間の種類"
  ) +
  scale_color_manual(values = c(
    "通常の95%信頼区間" = "steelblue",
    "ガブリエル比較区間" = "firebrick"
  )) +
  theme_minimal()

print(p1)