本書第 5 章「GLM の尤度比検定と検定の非対称性」をもとにした。
この章では、どのような統計モデルでも利用可能な尤度比検定（likelihood ratio test）について説明します。これは前のモデル選択の章に登場した逸脱度の差に注目する考えかたです。P.94
本章の検定と、前章の「AIC モデル選択」との目的の違いを、「5.6 検定とモデル選択、そして推定された統計モデルの解釈」をもとに記す
AIC によるモデル選択の目的
「予測の良さとは平均対数尤度」として「良い予測をするモデル」を選ぶのが目的。平均対数尤度を最大対数尤度とパラメーター数から推定する。
尤度比検定などの Nayman-Pearson のもとでの統計学的な検定の目的
帰無仮説の安全な棄却が目的。帰無仮説が棄却されたあとの対立仮説が、どのような意味で「良い」モデルなのかは明確ではありません。
データ分析では、この違いの理解、そして以下の認識は重要。
この章の主題は、統計学的な検定が「とにかく p < 0.05 さえ出せばいい、そうすれば何を主張してもよい」といった万能のツールではないことをあらためて確認するものでした。P.111
検定手順
Nayman-Pearson の統計学的検定のわくぐみでは、パラメーター数の少ないモデルを帰無仮説と位置づけ、帰無仮説が棄却できるかどうかの確率評価に専念する。P.110
ここでは、「GLMの目的と特徴」と同様に平均種子数のモデルを用いて、以下の帰無仮説と対立仮説で検定する。
- 帰無仮説：「一定モデル」、種子数の平均 λi が定数（体サイズに非依存）モデル（傾き β2 = 0, パラメーター数 = 1）
- 対立仮説：「x モデル」、種子数の平均 λi が体サイズに依存するモデル（傾き β2 ≠ 0, パラメーター数 = 2）
ここでの 尤度比（likelihood ratio）は
一定モデルの最大尤度 / x モデルの最大尤度
ΔD1,2 = (-2 log L1) - (-2 log L2) = -2 × (log L1 - log L2)
同じデータに対してパラメーター数の多いモデルの方が逸脱度は常に小さくなる。
> fit.2 <- glm(y ~ x,data=d,family = poisson(link = "log"))
> logLik(fit.1);logLik(fit.2)
'log Lik.' -237.6432 (df=1)
'log Lik.' -235.3863 (df=2)
> -2*as.numeric(logLik(fit.1) - logLik(fit.2))
[1] 4.513941
ΔD1,2 は約 4.5 で、今回の尤度比検定では
検定統計量であるこの尤度比の差（ΔD1,2）が「4.5 ぐらいでは改善されていない」と言ってよいのかを調べます。P.99
検定方法は次のふた通りで
- パラメトリックブートスラップ法
- χ2 分布を使った近似計算法
観測データ数が少ない場合は、前者の「パラメトリックブートスラップ法」を用いること。
パラメトリックブートスラップ（parametric bootstrap, PB）法
いかなるめんどうな状況でも必ず p 値が計算できる PB 法 P.102
この方法を、検定の流れで説明する。
> fit.1$coefficients[1]
(Intercept)
2.057963
> exp(2.058); mean(d$y)
[1] 7.830294
[1] 7.83
先の帰無仮説の推定結果 exp(2.058) と 標本平均 mean(d$y) はほぼ等しい。
この標本平均で乱数生成（シミュレーションデータ作成）
> d$y.rnd <- rpois(100,lambda = mean(d$y))
作成したデータで、帰無仮説と対立仮説のモデルで推定して、逸脱度の差を求める。
> fit.12<-glm(y.rnd ~ 1,data=d,family=poisson)
> fit.22<-glm(y.rnd ~ x,data=d,family=poisson)
> fit.12$deviance - fit.22$deviance
[1] 0.2984923
体サイズ x と無意味な平均値一定のポアソン乱数データに対しても、帰無仮説である「一定モデル」の方があてはまりが良い fit.12$deviance > fit.22$deviance となっている。
このステップを、今回は 1,000 回実施する（pb 関数は末尾参照）
> setwd("/***/kubobook2012/chapter05/")
> source("pb.R")
> dd12<-pb(d,n.bootstrap = 1000)
> summary(dd12)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000001 0.091790 0.428500 0.964200 1.231000 11.510000
結果のヒストグラム（本書図 5.4 に相当）
> ggplot(data.frame(ΔD=dd12),aes(x=ΔD))+geom_histogram(binwidth = 0.2,fill="white",color="black")+geom_vline(xintercept = 4.5, color="red")+annotate("text",label="ΔD = 4.5",x=5.5,y=250,color="red")
逸脱度の差 4.5 より大きい割合は
> sum(dd12>=4.5)/length(dd12)
[1] 0.03
95% の範囲にある「逸脱度の差」は
> quantile(dd12, 0.95)
95%
3.78107
よって、ΔD1,2 ≤ 3.78107 は、有意水準 5 %の統計学的検定のわくぐみでは「よくある差」とみなされる。
結論は、逸脱度の差 4.5 の p 値 0.03 は、有意水準 0.05 より小さいので、有意差があり（significantly different）、帰無仮説の「一定モデル」は棄却され「x モデルを採択」となる。
χ2 分布を使った近似計算法
帰無仮説と対立仮説の推定結果から、anova 関数で
> fit.1 <- glm(y ~ 1,data=d,family = poisson(link = "log"))
> fit.2 <- glm(y ~ x,data=d,family = poisson(link = "log"))
> anova(fit.1, fit.2, test="Chisq")
Analysis of Deviance Table
Model 1: y ~ 1
Model 2: y ~ x
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 99 89.507
2 98 84.993 1 4.5139 0.03362 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
結果から、逸脱度の差 ΔD1,2 が 4.5139 になる p 値は 0.03362 で、帰無仮説は棄却できる。
ただし、先の PB 法での p 値 0.03 とは一致していない。これは、
χ2 分布近似はサンプルサイズが大きい場合に有効な近似計算であり、この例題で調べた植物の個体数は 100 にすぎないので、このように "Chisq" 指定によって近似的に得られた p 値はあまり正確でない可能性があります。P.107
「ロジスティック回帰分析」に続く。
pb 関数
function(d, n.bootstrap)
{
n.sample <- nrow(d)
y.mean <- mean(d$y)
cat("# ")
v.d.dev12 <- sapply(
1:n.bootstrap,
function(i) {
cat(".")
if (i %% 50 == 0) cat("\n# ")
d$y.rnd <- rpois(n.sample, lambda = y.mean)
fit1 <- glm(y.rnd ~ 1, data = d, family = poisson)
fit2 <- glm(y.rnd ~ x, data = d, family = poisson)
fit1$deviance - fit2$deviance
}
)
cat("\n")
v.d.dev12
}
0 件のコメント:
コメントを投稿