Rの関数:wald.test {aod}

関数 wald.test とは

wald.test は、統計モデルにおいて推定されたパラメータに関する「ワルド検定（Wald Test）」を実行するための関数です。

ワルド検定は、最尤法などで得られたパラメータの推定値が、帰無仮説で想定される値（通常はゼロ）から統計的に有意に離れているかどうかを評価する、尤度比検定（Likelihood Ratio Test）やスコア検定（Score Test）と並ぶ、漸近検定の一つです。

本関数は、モデルの適合結果オブジェクト（glm や lm などの出力）から直接計算を行うのではなく、ユーザーが「パラメータの推定値ベクトル」と「分散共分散行列」を明示的に渡す仕様となっています。

当該仕様により、推定値と分散さえわかれば柔軟に仮説検定を適用することができます。

関数 wald.test の引数

library(aod)
args(wald.test)

function (Sigma, b, Terms = NULL, L = NULL, H0 = NULL, df = NULL, 
    verbose = FALSE) 
NULL

1. Sigma
   • パラメータ推定値の分散共分散行列です。
   • 通常、モデルの適合結果に対して vcov() 関数を適用することで取得できます。行数および列数は、後述するパラメータベクトル b の長さと一致している必要があります。
2. b
   • パラメータの推定値ベクトルです。
   • モデルの適合結果に対して coef() 関数を適用することで取得できます。
3. Terms
   • 検定対象となるパラメータのインデックス（位置）を指定する整数ベクトルです。
   • 例えば、モデルに5つのパラメータがあり、2番目と3番目のパラメータが同時にゼロであるか（結合仮説）を検定したい場合、Terms = c(2, 3) または Terms = 2:3 と指定します。後述の L とは排他的であり、どちらか一方のみを使用します。
4. L
   • パラメータの線形結合を指定する対比行列（Contrast matrix）です。
   • パラメータ間の複雑な関係性を検定する際に用います。行列の列数はパラメータベクトル b の長さと等しく、行数は検定する仮説の数となります。例えば「パラメータ1とパラメータ2の値が等しいか（パラメータ1 \(-\) パラメータ2 \(=\) 0）」を検定したい場合、該当列に 1 と -1 を配置した行列を指定します。
5. H0
   • 帰無仮説において想定される値のベクトルです。
   • デフォルト: NULL。自動的にすべての要素が 0 のベクトルとして扱われます。
6. df
   • F検定を行うための、分母の自由度を指定する数値です。
   • デフォルト: NULL。
   • 値を指定しない場合、ワルド統計量はカイ二乗分布に従うとしてP値が計算されます。線形回帰などで残差自由度が明確な場合、当該数値を指定することで、より小標本に適したF検定のP値も併せて出力されます。
7. verbose
   • 計算の途中経過や詳細な情報をコンソールに出力するかどうかの論理値です。

シミュレーションコード

以下に、wald.test 関数の機能を確認するためのシミュレーションコードを示します。

本シミュレーションでは、2つの説明変数を持つロジスティック回帰モデルを構築します。

データ生成の段階で、2つの説明変数が応答変数に与える影響（真の係数）を意図的に「全く同じ値」に設定します。

当該データに対してモデルを適合させた後、wald.test を用いて以下の2つの仮説を検定します。

1. 結合仮説の検定 (Terms の利用): 2つの変数の係数が、同時にゼロであるかどうか。
2. 線形結合の検定 (L 行列の利用): 2つの変数の係数が、互いに等しいかどうか。

なお、有意水準は5%とします。

サンプルデータの生成

# パッケージの読み込み
library(ggplot2)

# 乱数シードの固定
seed <- 20260425
set.seed(seed)

# 1. サンプルデータの生成
n_obs <- 500

# 説明変数の生成 (標準正規分布)
var_X1 <- rnorm(n_obs, mean = 0, sd = 1)
var_X2 <- rnorm(n_obs, mean = 0, sd = 1)

# 真のパラメータ設定
# 切片を -0.5、X1の係数を 0.8、X2の係数も 0.8 とします
true_beta0 <- -0.5
true_beta1 <- 0.8
true_beta2 <- 0.8

# ロジスティック回帰のデータ生成プロセス
linear_predictor <- true_beta0 + true_beta1 * var_X1 + true_beta2 * var_X2
probabilities <- 1 / (1 + exp(-linear_predictor))
response_Y <- rbinom(n_obs, size = 1, prob = probabilities)

# データフレームの作成
df_sim <- data.frame(Y = response_Y, X1 = var_X1, X2 = var_X2)

cat("--- データ概要 ---\n")
cat(sprintf("観測数: %d\n", n_obs))
cat(sprintf("真の係数: 切片=%.1f, X1=%.1f, X2=%.1f\n", true_beta0, true_beta1, true_beta2))

--- データ概要 ---
観測数: 500
真の係数: 切片=-0.5, X1=0.8, X2=0.8

# 2. ロジスティック回帰モデルの適合
model_fit <- glm(Y ~ X1 + X2, data = df_sim, family = binomial(link = "logit"))

# パラメータ推定値 (b) と 分散共分散行列 (Sigma) を抽出
est_b <- coef(model_fit)
est_Sigma <- vcov(model_fit)

cat("--- モデルの推定結果 ---\n")
print(round(est_b, 4))

--- モデルの推定結果 ---
(Intercept)          X1          X2 
    -0.5859      0.7622      0.7693 

係数の可視化

# 3. 係数の可視化
# 推定値と95%信頼区間を計算
se <- sqrt(diag(est_Sigma))
df_coef <- data.frame(
  Variable = names(est_b),
  Estimate = est_b,
  Lower = est_b - 1.96 * se,
  Upper = est_b + 1.96 * se
)

p1 <- ggplot(df_coef, aes(x = Variable, y = Estimate)) +
  geom_point(size = 4, color = "blue") +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, linewidth = 1, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "ロジスティック回帰モデルの係数推定値と95%信頼区間",
    x = "説明変数",
    y = "係数の推定値"
  ) +
  theme_minimal()

print(p1)

Figure 1 から、X1とX2の推定値は近似しており、かつ、共に95%信頼区間が0を跨がないことを確認できます。

# 4. wald.test を用いた仮説検定の実行

cat("=== 検定1: X1とX2の係数は同時にゼロであるか？ (Terms の利用) ===\n")
# 切片が1番目、X1が2番目、X2が3番目のパラメータです。
# 2番目と3番目を同時に検定するため、Terms = 2:3 を指定します。
test_joint <- wald.test(Sigma = est_Sigma, b = est_b, Terms = 2:3)
print(test_joint)

=== 検定1: X1とX2の係数は同時にゼロであるか？ (Terms の利用) ===
Wald test:
----------

Chi-squared test:
X2 = 80.6, df = 2, P(> X2) = 0.0

cat("=== 検定2: X1の係数とX2の係数は等しいか？ (L 行列の利用) ===\n")
# 帰無仮説: beta1 = beta2  =>  0*beta0 + 1*beta1 - 1*beta2 = 0
# 行列 L を (0, 1, -1) として定義します。
contrast_matrix <- matrix(c(0, 1, -1), nrow = 1)
test_equality <- wald.test(Sigma = est_Sigma, b = est_b, L = contrast_matrix)
print(test_equality)

=== 検定2: X1の係数とX2の係数は等しいか？ (L 行列の利用) ===
Wald test:
----------

Chi-squared test:
X2 = 0.0023, df = 1, P(> X2) = 0.96