重回帰分析におけるBenjamini-Hochberg法

Benjamini-Hochberg法の重回帰での意義

重回帰分析においてBenjamini-Hochberg法（BH法）を適用する主な意義は、多重検定問題に対処することである。重回帰分析では複数の説明変数の係数に対して同時に仮説検定を行うため、偽陽性（第一種の過誤）が増加するリスクがある。BH法は偽発見率（FDR）を制御することで、この問題に対処する。

重回帰分析でBH法を使用する利点は以下の通りである：

• 偽陽性の制御: 多数の変数を扱う際に、偶然による有意な結果の出現を抑制する。
• Bonferroni法より検出力が高い: より多くの真の関連を検出できる可能性がある。
• モデル選択の補助: 有意な変数の選択に役立つ可能性がある。

BH法を使用した場合と使用しない場合の主な違い

• 有意と判断される変数の数: BH法を適用すると、通常、有意と判断される変数の数が減少することがある。これは偽陽性を制御する効果である。
• p値の解釈: オリジナルのp値は個々の検定の有意性を示すが、BH法で調整したp値は全体の文脈での有意性を反映する。
• 変数選択への影響: BH法の適用により、モデルに含める変数の選択が変わる可能性がある。これは、より保守的な変数選択につながる。
• 結果の信頼性: BH法を適用することで、多重検定の問題に対処し、結果の信頼性が向上する。

重回帰分析をつくる

# 必要なライブラリのロード
library(AER)
library(stats)

# データのロード
data(CPS1985)

# データの概要を確認
summary(CPS1985)

# 回帰分析の実行
model <- lm(wage ~ education + experience + region + gender + union, data = CPS1985)

# 回帰係数とp値の取得
coef_summary <- summary(model)$coefficients

# p値の抽出（切片を除く）
p_values <- coef_summary[-1, "Pr(>|t|)"]

Benjamini-Hochberg法の適用

p_adjusted <- p.adjust(p_values, method = "BH")

# 結果の表示
results <- data.frame(
  Variable = rownames(coef_summary)[-1],
  Coefficient = coef_summary[-1, "Estimate"],
  Std_Error = coef_summary[-1, "Std. Error"],
  Original_P = p_values,
  Adjusted_P = p_adjusted,
  Significant_Original = p_values < 0.05,
  Significant_BH = p_adjusted < 0.05
)

print(results)

結果。

                 Variable Coefficient  Std_Error   Original_P   Adjusted_P Significant_Original Significant_BH
education       education   0.9130576 0.07907041 1.168831e-27 5.844153e-27                 TRUE           TRUE
experience     experience   0.1054496 0.01672143 6.046634e-10 1.511658e-09                 TRUE           TRUE
regionother   regionother   0.7852376 0.42667738 6.627722e-02 6.627722e-02                FALSE          FALSE
genderfemale genderfemale  -2.1699120 0.39028263 4.292549e-08 7.154248e-08                 TRUE           TRUE
unionyes         unionyes   1.3883889 0.51015930 6.713589e-03 8.391987e-03                 TRUE           TRUE

補正前と補正後のP値を比較した。

# 結果を比較しやすいようにソートして表示
results_sorted <- results[order(results$Original_P), ]

# P値の比較表を作成
comparison_table <- data.frame(
  Variable = results_sorted$Variable,
  Original_P = results_sorted$Original_P,
  Adjusted_P = results_sorted$Adjusted_P,
  Difference = results_sorted$Adjusted_P - results_sorted$Original_P,
  Original_Significant = results_sorted$Original_P < 0.05,
  BH_Significant = results_sorted$Adjusted_P < 0.05,
  Changed = results_sorted$Original_P < 0.05 & results_sorted$Adjusted_P >= 0.05
)

# 結果の表示
print(comparison_table)

大きくは変わらないようだ

      Variable   Original_P   Adjusted_P   Difference Original_Significant BH_Significant Changed
1    education 1.168831e-27 5.844153e-27 4.675323e-27                 TRUE           TRUE   FALSE
2   experience 6.046634e-10 1.511658e-09 9.069951e-10                 TRUE           TRUE   FALSE
3 genderfemale 4.292549e-08 7.154248e-08 2.861699e-08                 TRUE           TRUE   FALSE
4     unionyes 6.713589e-03 8.391987e-03 1.678397e-03                 TRUE           TRUE   FALSE
5  regionother 6.627722e-02 6.627722e-02 0.000000e+00                FALSE          FALSE   FALSE

有意性が変化した変数の数と最大の調整量を表示

# 有意性が変化した変数の数を表示
cat("\n有意性が変化した変数の数:", sum(comparison_table$Changed))

# 最大の調整量を表示
max_adjustment <- max(comparison_table$Difference)
cat("\n最大の調整量:", max_adjustment)

有意性が変化した変数の数: 0 最大の調整量: 0.001678397