ポアソン回帰における逸脱度検定の仕組み
カウントデータを扱う際の標準的な手法としてポアソン回帰モデルが広く利用されています。多くの統計ソフトウェアでは、モデル適合後に自動的に適合度検定の結果を提供するか、容易に実行できる機能を持っています。しかし、この検定結果を無条件に信頼することは危険であり、その限界を理解しておく必要があります。
モデルの適合度を評価する指標として「逸脱度（Deviance）」が用いられます。これは、現在のモデルと飽和モデル（各観測値ごとにパラメータを持つモデル）との対数尤度比に基づいています。理論的には、モデルが正しく指定されている場合、この逸脱度は自由度が「観測値数 - パラメータ数」のカイ二乗分布に従うとされています。
したがって、逸脱度の値に対してカイ二乗分布の右側確率を計算することで、モデルの適合度に関する p 値を得ることができます。null 仮説は「モデルが正しく指定されている」というものです。
R 言語での実装例
以下に、R 言語を用いてポアソン回帰モデルを適合し、逸脱度に基づいた適合度検定を手動で計算する手順を示します。まず、一様分布する共変量と、それに基づくポアソン分布する応答変量を模擬生成します。
set.seed(2023)
sample_size <- 1000
covariate <- runif(sample_size)
lambda_vec <- exp(covariate)
count_data <- rpois(sample_size, lambda_vec)
glm_model <- glm(count_data ~ covariate, family = poisson)
summary(glm_model)
出力結果には「Residual deviance」（残差逸脱度）と「df.residual」（残差自由度）が含まれています。例えば、残差逸脱度が 1110.3、自由度が 998 である場合、適合度検定の p 値は以下のように計算します。
p_value <- pchisq(glm_model$deviance,
df = glm_model$df.residual,
lower.tail = FALSE)
print(p_value)
この計算により得られた p 値が有意水準（例えば 0.05）を下回る場合、モデルの指定が不適切であるという証拠があると解釈されます。
シミュレーションによる検定性能の検証
この検定が理論通りに機能するかを確認するため、モンテカルロシミュレーションを実施します。既知の正しいモデルからデータを生成し、それがどの程度の頻度で「不適切」と判定されるか（第 1 種の過誤）を調べます。
まず、平均カウント数が比較的小さい状況で 10,000 回の試行を行います。
simulation_runs <- 10000
p_value_storage <- numeric(simulation_runs)
for (i in seq_len(simulation_runs)) {
n <- 1000
x_var <- runif(n)
mu <- exp(x_var)
y_var <- rpois(n, mu)
fit <- glm(y_var ~ x_var, family = poisson)
p_value_storage[i] <- pchisq(fit$deviance,
df = fit$df.residual,
lower.tail = FALSE)
}
rejection_rate <- sum(p_value_storage < 0.05) / simulation_runs
print(rejection_rate)
このシミュレーションを実行すると、 rejection_rate は約 0.94 となります。これは、モデルが正しく指定されているにもかかわらず、約 94% のケースで適合度検定がモデルを棄却してしまうことを意味します。つまり、平均値が小さい個数データにおいて、逸脱度のカイ二乗近似は精度が非常に低いです。
次に、期待されるカウント数が larger な場合で同様の実験を行います。平均値が 20 から 55 の範囲になるようにデータ生成過程を変更します。
simulation_runs <- 10000
p_value_storage <- numeric(simulation_runs)
for (i in seq_len(simulation_runs)) {
n <- 1000
x_var <- runif(n)
mu <- exp(3 + x_var)
y_var <- rpois(n, mu)
fit <- glm(y_var ~ x_var, family = poisson)
p_value_storage[i] <- pchisq(fit$deviance,
df = fit$df.residual,
lower.tail = FALSE)
}
rejection_rate <- sum(p_value_storage < 0.05) / simulation_runs
print(rejection_rate)
この条件下では、棄却率は約 0.0635 に低下し、名目上の 5% という有意水準に近づきます。これは、期待値が十分に大きければ、逸脱度のカイ二乗近似が妥当になることを示唆しています。