尤度比検定
モデル間の尤度を比較する
尤度likelihoodとは、あるモデル（データが従う確率分布、パラメータ）を仮定した時に、手元のデータがそのモデルから生成される確率、あるいは確率密度に比例する値として定義されるのであった。この尤度を用いたパラメータ推定方法が最尤法most likelihood methodであり、その名の通り尤度をパラメータの関数としてみたときに、最も大きな尤度を実現するパラメータを推定値とする方法である。尤度は、手元のデータを実現する確率（に準ずるもの）として考えられるのだから、最も大きな「確率」を実現できるパラメータを選ぶという方法は手元のデータの数値に対して最も適合するものを求める方法として、至極自然な考えである。
さて今まである説明変数がモデルの予測の改善に役立ちそうかを判断する材料として、Wald統計量や統計的仮説検定とは異なる視点から議論したAICを紹介した。本項ではこれらとはまた異なる視点で、説明変数が有意な効果を持つか判断する材料として尤度比検定likelihood ratio testを紹介する。尤度比検定は尤度を計算する推定方法であればあまねく利用可能な便利な方法である。また、glm関数などの出力を要約して判断するうえでも重要な手法となる。
尤度比検定likelihood ratio test（LRT）
尤度比検定はその名の通り、比較したい2つのモデルの尤度を比較するときに用いられ、その検定統計量は尤度の比に比例するものとなる。単純な2仮説における比較においては、尤度比検定は一様最強検出力検定Uniformly most powerful testであると表現される。つまり、共通の危険率を持つ検定群において最も検出力が高い。ただし、あらゆるモデル間を比較できるのではなく、一般化線形モデルなどで説明変数がモデルの改善に役立ったかを尤度比検定するときは、一方のモデルは他方のモデルのネストとなっている必要がある。例えば下図なら、2つのパラメータβ0、β1をもつモデルが、1つのパラメータβ0しか持たないモデルをネストしている状態にある。
このとき、ネストしている側のモデルをフルモデルfull model、されている側のモデルを縮小モデルreduced modelと呼び、一般的な統計的仮説検定に倣って言えば、縮小モデル側が帰無仮説、フルモデル側が対立仮説に相当する。上記の例に立ち戻って考えてみると、縮小モデル側はフルモデルにおいてβ1 = 0を仮定している状態と考えられるので、縮小モデルとフルモデルの尤度比検定は、β1 = 0を帰無仮説とする統計的仮説検定である、ということである。ただし、ややこしいことを言うようで恐縮だが、尤度比検定における帰無仮説は縮小モデルとフルモデルでは尤度の差がない、ということである。
ここで尤度比検定の検定統計量を示す。以下のとおりである。
尤度比検定においては、尤度比Λそのものではなく、サンプルサイズが十分に大きければ-2 log(Λ)がフルモデルと縮小モデルの自由度の差を自由度とするカイ二乗分布に従うことを利用して検定を行う。モデル間の自由度の差とは、基本的には推定するパラメータの数の差である。あるいは後述のように近似ではなくシミュレーションによって-2 log(Λ)の分布を生成してブートストラップ法を通じて検定することも可能である。log(Λ)は整理すれば、縮小モデルの最大対数尤度からフルモデルの最大対数尤度を引いたものである。一般にパラメータの数が多いフルモデルのほうが縮小モデルよりも尤度が大きいため、log(Λ)は負の値を示す。ゆえに、-2 log(Λ)は正の値を示す統計量である。AICの項を見てもらえばわかるが、-2 log(Λ)とは逸脱度残差deviance residualsそのものである。
尤度比検定は、その検定統計量を見ればわかるが、あるパラメータを追加して推定した時に、その最大尤度がどれだけ改善されたかを判定する検定である。Wald検定のようにパラメータの効果が有意に0からずれているかというもの以上に、「よりモデルが改善したかどうかを判断する」検定であると考えてよいだろう。
尤度比が従う確率分布
尤度比検定自体は、どんなモデル間比較でも行われる可能性があるが、今回は、その代表格である一般化線形モデルを具体例として紹介しよう。
・連続説明変数が1つ
では、いくつかのシミュレーションを通じてモデル間の尤度比の性質を確かめてみよう。例えば、以下のような切片項だけを持つモデル（縮小モデル）と切片項と連続変数xの係数のパラメータを持つモデル（フルモデル）の比較を考える。今回は簡単のため、1つのパラメータで記述できるポアソン分布に従うデータとする。ここで、連続変数xに関する尤度比検定を考える。さて、尤度比検定の帰無仮説は縮小モデルとフルモデルの尤度比が同じである、つまり連続変数xにかかる係数が0であるという前提をもとにデータを生成する。
------------------------------------------------------
library(plotn)
library(car)
b0 <- 0
b1 <- 0 #連続変数xにかかる真のパラメータ
n <- 100
x <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x))
fitN <- glm(y ~ 1, family = poisson(link = "log")) #reduced model
fitF <- glm(y ~ x, family = poisson(link = "log")) #full model
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN2
##
## Call:
## glm(formula = y ~ 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.42127 -1.42127 -0.00997 0.86763 2.83106
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.00995 0.09950 0.1 0.92
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 129.54 on 99 degrees of freedom
## Residual deviance: 129.54 on 99 degrees of freedom
## AIC: 272.93
##
## Number of Fisher Scoring iterations: 5
fitF2
##
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.49352 -1.40486 -0.05967 0.77230 2.72831
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.112730 0.218471 -0.516 0.606
## x 0.007724 0.012045 0.641 0.521
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 129.54 on 99 degrees of freedom
## Residual deviance: 129.13 on 98 degrees of freedom
## AIC: 274.52
##
## Number of Fisher Scoring iterations: 5
------------------------------------------------------
さて問題は、連続変数xの追加はモデルの改善に役立ったのだろうか。それを検定するのが尤度比検定である。実際には、xにかかる真のパラメータb1は0なのだから、大体は効果なし＝役立たないだろうと予想される。
ここで以下のように推定されたパラメータから対数尤度を計算する関数を定義する。xには説明変数をまとめた行列、yには被説明変数、btは各モデルから得られた推定値をまとめた行列をとる。
------------------------------------------------------
f <- function(x, y, bt){
if(ncol(bt) > 1){
colSums(dpois(x = y, lambda = exp(X%*%bt), log = T))
#行列の積を使い、まとめて対数尤度を計算
} else {
sum(dpois(x = y, lambda = exp(X%*%bt), log = T))
}
}
------------------------------------------------------
上記を使い、以下のように縮小モデルとフルモデルのパラメータ推定値を使い、対数尤度を計算する。出力で得られるベクトルのうち、1項目は縮小モデル、2項目はフルモデルの最大対数尤度である。
------------------------------------------------------
b0eN <- coef(fitN)[1] #reduced modelの切片項推定値
b0eF <- coef(fitF)[1] #full modelの切片項推定値
b1eF <- coef(fitF)[2] #full modelのxの係数推定値
X <- cbind(1,x) #1列目は切片項を表しすべて1とする。2列目は説明変数xそのもの。
bt <- cbind(c(b0eN,0),c(b0eF,b1eF)) #reduced modelのxの係数推定値は0とする
f(X, y, bt)
## [1] -135.4644 -135.2580
------------------------------------------------------
最大対数尤度の差を-2倍したものが尤度比検定の検定統計量である。これは、上記のfitN2とfitF2に格納されている逸脱度を使っても計算できる。さらに、fitF2におけるNull devianceとResidual devianceの差でもある。
------------------------------------------------------
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #-2 log(Λ)
## [1] 0.4128908
fitN$deviance - fitF$deviance
## [1] 0.4128908
------------------------------------------------------
さて、上記のように検定統計量が計算できたが、問題はこれがどんな確率分布に従っているかである。そこで、以下のように10000回、上記と同じ設定でデータを生成して、分布を確認してみることにする。
------------------------------------------------------
b0 <- 0
b1 <- 0
n <- 100
R <- NULL
for (i in 1:10000){
x <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
X <- cbind(1,x)
bt <- cbind(c(b0eN,0),c(b0eF,b1eF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図1の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 1),
type = "l", col = "red")
)
------------------------------------------------------
図1 -2 log(Λ)の分布と自由度1のカイ二乗分布（赤線）
すると-2 log(Λ)は自由度1のカイ二乗分布とよく一致することがわかる。ここで、カイ二乗分布の自由度は、モデル間の自由度の差で決まり、今回は縮小モデルでは推定するパラメータが1個、フルモデルでは2個であるので、その差の1がカイ二乗分布の自由度となる。したがって、帰無仮説（縮小モデル）が正しい＝b1が0のとき、-2 log(Λ)はカイ二乗分布に従うことを利用して統計的仮説検定が可能である。もし、-2 log(Λ)がカイ二乗分布の右端5%の領域に存在するなら、帰無仮説が正しいという仮定に無理がある、つまり説明変数xはモデルの改善に有効であると判断できるというわけである。
今回は、サンプルサイズ100で解析を行っているため、-2 log(Λ)はカイ二乗分布によく従う。しかし、サンプルサイズが小さい場合はその限りではないので、上記のようにシミュレーションによって分布を作成し、それを帰無分布とすることでブートストラップ確率を求めて、それをp値とすることもできる。
・連続説明変数が3つ
ほかの状況も確認しておこう。例えば、連続変数が3つになった時は、この3つすべてに関する-2 log(Λ)は以下のように計算できる。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b3 <- 0
n <- 100
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
x3 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b3*x3))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x1 + x2 + x3, family = poisson(link = "log"))
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN2
##
## Call:
## glm(formula = y ~ 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.39284 -1.39284 0.03031 0.03031 2.29653
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.03046 0.10153 -0.3 0.764
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 127.96 on 99 degrees of freedom
## Residual deviance: 127.96 on 99 degrees of freedom
## AIC: 267.67
##
## Number of Fisher Scoring iterations: 5
fitF2
##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.64266 -1.32085 0.00037 0.42098 2.21261
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.300398 0.349254 -0.860 0.3897
## x1 -0.010967 0.012655 -0.867 0.3861
## x2 0.007928 0.011439 0.693 0.4883
## x3 0.020334 0.012117 1.678 0.0933 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 127.96 on 99 degrees of freedom
## Residual deviance: 122.98 on 96 degrees of freedom
## AIC: 268.69
##
## Number of Fisher Scoring iterations: 5
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0eN,0,0,0),c(b0eF,b1eF,b2eF,b3eF))
f(X, y, bt)
## [1] -132.8337 -130.3427
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #-2 log(Λ)
## [1] 4.982034
fitN$deviance - fitF$deviance
## [1] 4.982034
------------------------------------------------------
このときの-2 log(Λ)の分布は以下の通り。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b3 <- 0
n <- 100
R <- NULL
for (i in 1:10000){
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
x3 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b3*x3))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x1 + x2 + x3, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0eN,0,0,0),c(b0eF,b1eF,b2eF,b3eF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図2の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 3),
type = "l", col = "red")
)
------------------------------------------------------
図2 -2 log(Λ)の分布と自由度3のカイ二乗分布（赤線）
すると-2 log(Λ)は自由度3のカイ二乗分布とよく一致することがわかる。今回は縮小モデルでは推定するパラメータが1個、フルモデルでは4個であるので、その差の3がカイ二乗分布の自由度となる。
・連続説明変数が5つ
さらに連続変数が5つになった時、これはシミュレーションだけで済ませてしまう。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b3 <- 0
b4 <- 0
b5 <- 0
n <- 100
R <- NULL
for (i in 1:10000){
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
x3 <- runif(n = n, 0, 30)
x4 <- runif(n = n, 0, 30)
x5 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b3*x3 + b4*x4 + b5*x5))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x1 + x2 + x3 + x4 + x5, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
b4eF <- coef(fitF)[5]
b5eF <- coef(fitF)[6]
X <- cbind(1,x1,x2,x3,x4,x5)
bt <- cbind(c(b0eN,0,0,0,0,0),c(b0eF,b1eF,b2eF,b3eF,b4eF,b5eF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図3の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 5),
type = "l", col = "red")
)
------------------------------------------------------
図3 -2 log(Λ)の分布と自由度5のカイ二乗分布（赤線）
やはり、-2 log(Λ)は自由度5のカイ二乗分布とよく一致することがわかる。今回は縮小モデルでは推定するパラメータが1個、フルモデルでは6個であるので、その差の5がカイ二乗分布の自由度となる。
・連続説明変数が2つとその交互作用
では、ちょっとわかりにくくなる例を紹介しよう。連続変数が2つとし、これらの交互作用を含むモデルをフルモデルとする。縮小モデルは切片だけだ。シミュレーションによって、これらのモデル間の尤度比がどんな分布に従うか調べる。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b12 <- 0 #交互作用項の係数
n <- 100
R <- NULL
for (i in 1:10000){
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b12*x1*x2))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x1*x2, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b12eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x1*x2)
bt <- cbind(c(b0eN,0,0,0),c(b0eF,b1eF,b2eF,b12eF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図4の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 1),
type = "l", col = "red"),
points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 2),
type = "l", col = "blue"),
points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 3),
type = "l", col = "green"),
points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 4),
type = "l", col = "purple"),
points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 5),
type = "l", col = "lightblue")
)
------------------------------------------------------
図4 -2 log(Λ)の分布と自由度1～5のカイ二乗分布（赤、青、緑、紫、薄青）
交互作用項は2つの説明変数の組み合わせだが、パラメータは1つで表されている。ゆえに、フルモデルの自由度は4、縮小モデルの自由度は1である。ゆえに、-2 log(Λ)は自由度3のカイ二乗分布に従うことが期待されるが、確かにそうなっていることがわかるだろう。交互作用項はやや特殊な項だが、あくまでも重要なのはパラメータの数であることがわかる。
・カテゴリカル説明変数（3水準）が1つ
次にカテゴリカル変数が説明変数の時の-2 log(Λ)の分布を確認してみよう。まずは3水準のカテゴリカル変数が1つの場合である。
------------------------------------------------------
b0 <- 0 #切片
bB <- 0 #水準Bの係数
bC <- 0 #水準Cの係数
n <- 120
g1 <- rep(LETTERS[1:3], each = n/3) #カテゴリカル説明変数
------------------------------------------------------
ここで乱数の生成方法だが、カテゴリカル変数の場合、ダミー変数が必要になるのであった。連続変数なら、上記のように行列Xを生成する方法で十分だったが、カテゴリカル変数の数が多くなるとダミー変数を手作業で生成するのはかなりめんどくさい。そこで、model.matrix関数を使えば簡単に生成できる。
------------------------------------------------------
X <- model.matrix( ~ g1)
head(X)
## (Intercept) g1B g1C
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 0 0
## 6 1 0 0
b <- c(b0,bB,bC)
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log")) #reduced model
fitF <- glm(y ~ g1, family = poisson(link = "log")) #full model
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ x, family = poisson(link = "log"))
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN2
##
## Call:
## glm(formula = y ~ 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.40831 -1.40831 0.00836 0.88848 2.26730
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.008368 0.091670 -0.091 0.927
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 133.04 on 119 degrees of freedom
## Residual deviance: 133.04 on 119 degrees of freedom
## AIC: 310.86
##
## Number of Fisher Scoring iterations: 5
fitF2
##
## Call:
## glm(formula = y ~ g1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5000 -1.2450 -0.0732 0.7426 2.5845
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2549 0.1796 -1.419 0.156
## g1B 0.3727 0.2334 1.597 0.110
## g1C 0.3272 0.2356 1.389 0.165
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 133.04 on 119 degrees of freedom
## Residual deviance: 130.03 on 117 degrees of freedom
## AIC: 311.85
##
## Number of Fisher Scoring iterations: 5
------------------------------------------------------
さて、各カテゴリの水準の係数は0であるから、カテゴリカル変数g1はモデルの改善に役立たないだろうと予測される。
今までのように、以下のように縮小モデルとフルモデルのパラメータ推定値を使い、対数尤度、-2 log(Λ)を計算する。
------------------------------------------------------
b0eN <- coef(fitN)[1] #reduced modelの切片項推定値
b0eF <- coef(fitF)[1] #full modelの切片項推定値
bBeF <- coef(fitF)[2] #full modelの水準Bの推定値
bCeF <- coef(fitF)[3] #full modelの水準Cの推定値
X <- model.matrix( ~ g1) #ダミー変数で構築された説明変数の行列
bt <- cbind(c(b0eN,0,0),c(b0eF,bBeF,bCeF)) #reduced modelの水準BとC推定値は0とする
f(X, y, bt)
## [1] -154.4312 -152.9270
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #-2 log(Λ)
## [1] 3.008374
fitN$deviance - fitF$deviance
## [1] 3.008374
------------------------------------------------------
同様に10000回、上記と同じ設定でデータを生成して、-2 log(Λ)分布を確認してみることにする。
------------------------------------------------------
b0 <- 0 #切片
bB <- 0 #水準Bの係数
bC <- 0 #水準Cの係数
n <- 120
R <- NULL
g1 <- rep(LETTERS[1:3], each = n/3)
for (i in 1:10000){
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ g1, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bCeF <- coef(fitF)[3]
X <- model.matrix( ~ g1)
bt <- cbind(c(b0eN,0,0),c(b0eF,bBeF,bCeF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図5の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 2),
type = "l", col = "red")
)
------------------------------------------------------
図5 -2 log(Λ)の分布と自由度2のカイ二乗分布（赤線）
すると-2 log(Λ)は自由度2のカイ二乗分布とよく一致することがわかる。注意したいのは、カテゴリカル説明変数は1つだが、その中に3水準あるため、フルモデルでは推定するべきパラメータが3個ある。縮小モデルでは1つだ。ゆえに、その差の2がカイ二乗分布の自由度となる。カイ二乗分布の自由度は、説明変数の数ではなく、推定するべきパラメータの数の差であることのを覚えておこう。
・カテゴリカル説明変数（3水準と2水準）の2つ
さらに3水準と2水準の2つのカテゴリカル変数が説明変数の時を考えてみよう。
------------------------------------------------------
b0 <- 0 #切片
bB <- 0 #水準Bの係数
bb <- 0 #水準bの係数
bc <- 0 #水準cの係数
n <- 120
g1 <- rep(LETTERS[1:2], each = n/2) #2水準カテゴリカル説明変数
g2 <- rep(letters[1:3], times = n/3) #3水準カテゴリカル説明変数
X <- model.matrix( ~ g1 + g2)
head(X)
## (Intercept) g1B g2b g2c
## 1 1 0 0 0
## 2 1 0 1 0
## 3 1 0 0 1
## 4 1 0 0 0
## 5 1 0 1 0
## 6 1 0 0 1
b <- c(b0,bB,bb,bc)
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ g1 + g2, family = poisson(link = "log"))
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN2
##
## Call:
## glm(formula = y ~ 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.36015 -1.36015 0.07696 0.29939 2.35882
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07796 0.09491 -0.821 0.411
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 152.17 on 119 degrees of freedom
## Residual deviance: 152.17 on 119 degrees of freedom
## AIC: 313.19
##
## Number of Fisher Scoring iterations: 5
fitF2
##
## Call:
## glm(formula = y ~ g1 + g2, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.59278 -1.35832 -0.00898 0.42389 2.36232
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.29977 0.20209 -1.483 0.138
## g1B 0.30874 0.19210 1.607 0.108
## g2b -0.08961 0.24458 -0.366 0.714
## g2c 0.22884 0.22649 1.010 0.312
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 152.17 on 119 degrees of freedom
## Residual deviance: 147.49 on 116 degrees of freedom
## AIC: 314.51
##
## Number of Fisher Scoring iterations: 5
------------------------------------------------------
今までのように、以下のように縮小モデルとフルモデルのパラメータ推定値を使い、対数尤度、-2 log(Λ)を計算する。
------------------------------------------------------
b0eN <- coef(fitN)[1] #reduced modelの切片項推定値
b0eF <- coef(fitF)[1] #full modelの切片項推定値
bBeF <- coef(fitF)[2] #full modelの水準Bの推定値
bbeF <- coef(fitF)[3] #full modelの水準bの推定値
bceF <- coef(fitF)[4] #full modelの水準cの推定値
X <- model.matrix( ~ g1 + g2) #ダミー変数で構築された説明変数の行列
bt <- cbind(c(b0eN,0,0,0),c(b0eF,bBeF,bbeF,bceF)) #reduced modelの水準BとC推定値は0とする
f(X, y, bt)
## [1] -155.5932 -153.2530
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #-2 log(Λ)
## [1] 4.680252
fitN$deviance - fitF$deviance
## [1] 4.680252
------------------------------------------------------
同様に10000回、上記と同じ設定でデータを生成して、-2 log(Λ)分布を確認してみることにする。
------------------------------------------------------
b0 <- 0 #切片
bB <- 0 #水準Bの係数
bb <- 0 #水準bの係数
bc <- 0 #水準cの係数
n <- 120
n <- 120
R <- NULL
g1 <- rep(LETTERS[1:2], each = n/2)
g2 <- rep(letters[1:3], times = n/3)
for (i in 1:10000){
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ g1 + g2, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bbeF <- coef(fitF)[3]
bceF <- coef(fitF)[4]
X <- model.matrix( ~ g1 + g2)
bt <- cbind(c(b0eN,0,0,0),c(b0eF,bBeF,bbeF,bceF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図6の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 3),
type = "l", col = "red")
)
------------------------------------------------------
図6 -2 log(Λ)の分布と自由度3のカイ二乗分布（赤線）
やはり、-2 log(Λ)は自由度3のカイ二乗分布とよく一致することがわかる。今回は縮小モデルでは推定するパラメータが1個、フルモデルでは4個であるので、その差の3がカイ二乗分布の自由度となる。
・カテゴリカル説明変数（3水準と2水準）の2つとそれらの交互作用
最後、カテゴリカル変数（2水準と3水準）とそれらの交互作用を含むモデルの-2 log(Λ)の分布を計算しよう。途中過程も必要なら追ってみるとよい。
------------------------------------------------------
b0 <- 0 #切片
bB <- 0 #水準Bの係数
bb <- 0 #水準bの係数
bc <- 0 #水準cの係数
bBb <- 0 #水準Bとbの交互作用係数
bBc <- 0 #水準Bとcの交互作用係数
n <- 120
g1 <- rep(LETTERS[1:2], each = n/2)
g2 <- rep(letters[1:3], times = n/3)
X <- model.matrix( ~ g1 * g2)
head(X)
## (Intercept) g1B g2b g2c g1B:g2b g1B:g2c
## 1 1 0 0 0 0 0
## 2 1 0 1 0 0 0
## 3 1 0 0 1 0 0
## 4 1 0 0 0 0 0
## 5 1 0 1 0 0 0
## 6 1 0 0 1 0 0
b <- c(b0,bB,bb,bc,bBb,bBc)
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ g1 * g2, family = poisson(link = "log"))
fitN2 <- summary(fitN)
fitF2 <- summary(fitF)
fitN2
##
## Call:
## glm(formula = y ~ 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3038 -1.3038 0.1582 0.1582 2.4679
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.16252 0.09901 -1.641 0.101
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 132.65 on 119 degrees of freedom
## Residual deviance: 132.65 on 119 degrees of freedom
## AIC: 291.15
##
## Number of Fisher Scoring iterations: 5
fitF2
##
## Call:
## glm(formula = y ~ g1 * g2, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.51658 -1.30384 0.05086 0.70499 2.32394
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2231 0.2500 -0.893 0.372
## g1B -0.5754 0.4167 -1.381 0.167
## g2b 0.1719 0.3393 0.506 0.613
## g2c 0.1178 0.3436 0.343 0.732
## g1B:g2b 0.7664 0.5193 1.476 0.140
## g1B:g2c 0.5182 0.5366 0.966 0.334
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 132.65 on 119 degrees of freedom
## Residual deviance: 125.85 on 114 degrees of freedom
## AIC: 294.35
##
## Number of Fisher Scoring iterations: 5
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bbeF <- coef(fitF)[3]
bceF <- coef(fitF)[4]
bBbeF <- coef(fitF)[5]
bBceF <- coef(fitF)[6]
X <- model.matrix( ~ g1 * g2)
bt <- cbind(c(b0eN,0,0,0,0,0),c(b0eF,bBeF,bbeF,bceF,bBbeF,bBceF))
f(X, y, bt)
## [1] -144.5767 -141.1760
-2*(f(X, y, bt)[1]-f(X, y, bt)[2])
## [1] 6.801411
fitN$deviance - fitF$deviance
## [1] 6.801411
b0 <- 0
bB <- 0
bb <- 0
bc <- 0
bBb <- 0
bBc <- 0
n <- 120
R <- NULL
g1 <- rep(LETTERS[1:2], each = n/2)
g2 <- rep(letters[1:3], times = n/3)
for (i in 1:10000){
y <- rpois(n = n, lambda = exp(X%*%b))
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fitF <- glm(y ~ g1 * g2, family = poisson(link = "log"))
b0eN <- coef(fitN)[1]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bbeF <- coef(fitF)[3]
bceF <- coef(fitF)[4]
bBbeF <- coef(fitF)[5]
bBceF <- coef(fitF)[6]
X <- model.matrix( ~ g1 * g2)
bt <- cbind(c(b0eN,0,0,0,0,0),c(b0eF,bBeF,bbeF,bceF,bBbeF,bBceF))
R <- c(R, -2*(f(X, y, bt)[1]-f(X, y, bt)[2]))
}
histn(R, freq = F) #図7の図示
overdraw(points(seq(0, 30, length = 200),
dchisq(seq(0, 30, length = 200), df = 5),
type = "l", col = "red")
)
------------------------------------------------------
図7 -2 log(Λ)の分布と自由度5のカイ二乗分布（赤線）
フルモデルでは6つのパラメータを推定し、縮小モデルでは1つパラメータを推定している。ゆえに、-2 log(Λ)は自由度5のカイ二乗分布に従う。
以上、徹底的に尤度比の従う分布に関して議論してきた。ここからはより具体的に、尤度比検定の手続きについて紹介していく。
Rで行う尤度比検定
・連続説明変数が3つのとき
では、Rで尤度比検定を行う。解析自体は非常に簡単で、GLMならcarパッケージ内のAnova関数を使うことで簡単に実行できる。デフォルトで実装されているanova関数は基本的に線形回帰限定でしたつかえないので注意が必要だ。
以下のような、いずれの説明変数もモデルの改善に役立たなそうな状況を考える。Anova関数に必要なのはフルモデルだけだ。フルモデルのglmオブジェクトをAnova関数に通すことで、いわゆる分散分析表のようなものが出力される。これが尤度比検定の結果となる。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b3 <- 0
n <- 100
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
x3 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b3*x3))
fitF <- glm(y ~ x1 + x2 + x3, family = poisson(link = "log")) #full model
fitF2 <- summary(fitF)
fitF2 #Wald検定
##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.57812 -1.33710 -0.06744 0.28616 2.91073
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.118618 0.311459 0.381 0.703
## x1 0.003359 0.011349 0.296 0.767
## x2 0.003070 0.011787 0.260 0.794
## x3 -0.012567 0.011878 -1.058 0.290
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 119.24 on 99 degrees of freedom
## Residual deviance: 117.87 on 96 degrees of freedom
## AIC: 274.14
##
## Number of Fisher Scoring iterations: 5
Anova(fitF) #尤度比検定
## Analysis of Deviance Table (Type II tests)
##
## Response: y
## LR Chisq Df Pr(>Chisq)
## x1 0.08750 1 0.7674
## x2 0.06791 1 0.7944
## x3 1.11672 1 0.2906
------------------------------------------------------
いずれの説明変数もp値が0.05より大きいから、有意な効果はないと判断できる。LR Chisqが今まで求めてきた尤度比に相当する-2 log(Λ)、Dfは説明変数の自由度、そしてPr(>Chisq)がp値である。Chisqという表記からカイ二乗分布を使っていることがわかる。
では、各値を求めていこう。x1の-2 log(Λ)とp値は以下の通り。フルモデルからx1だけ説明変数を除いて、再度GLMを行い、その係数を使って計算する。x1が除かれたモデルが縮小モデルである。この時、フルモデルと縮小モデルはパラメータ数の差が1だから、カイ二乗分布の自由度は1である。
------------------------------------------------------
fit <- glm(y ~ x2 + x3, family = poisson(link = "log")) #x1に関するreduced model
b0e <- coef(fit)[1]
b2e <- coef(fit)[2]
b3e <- coef(fit)[3]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0e,0,b2e,b3e),c(b0eF,b1eF,b2eF,b3eF))
f(X, y, bt)
## [1] -133.1118 -133.0680
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #x1に関する-2 log(Λ)
## [1] 0.08750475
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #x1の尤度比検定のp値
## [1] 0.7673736
------------------------------------------------------
確かに一致することを確認しよう。
x2も同様にフルモデルからx2だけ説明変数を除いて、再度GLMを行い、その係数を使って計算する。
------------------------------------------------------
fit <- glm(y ~ x1 + x3, family = poisson(link = "log")) #x2に関するreduced model
b0e <- coef(fit)[1]
b1e <- coef(fit)[2]
b3e <- coef(fit)[3]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0e,b1e,0,b3e),c(b0eF,b1eF,b2eF,b3eF))
f(X, y, bt)
## [1] -133.102 -133.068
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #x2に関する-2 log(Λ)
## [1] 0.06790598
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #x2の尤度比検定のp値
## [1] 0.7944104
------------------------------------------------------
最後、x3も同様である。
------------------------------------------------------
fit <- glm(y ~ x1 + x2, family = poisson(link = "log")) #x3に関するreduced model
b0e <- coef(fit)[1]
b1e <- coef(fit)[2]
b2e <- coef(fit)[3]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b3eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0e,b1e,b2e,0),c(b0eF,b1eF,b2eF,b3eF))
f(X, y, bt)
## [1] -133.6264 -133.0680
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #x3に関する-2 log(Λ)
## [1] 1.116715
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #x3の尤度比検定のp値
## [1] 0.2906268
------------------------------------------------------
注意したいのはヌルモデル（切片だけ）に説明変数を加えたモデルとの尤度比検定ではないという点である。つまり、下記のように計算しても必要な値は求まらない。
------------------------------------------------------
#以下は間違い
fitN <- glm(y ~ 1, family = poisson(link = "log"))
fit <- glm(y ~ x3, family = poisson(link = "log"))
b0e <- coef(fit)[1]
b3e <- coef(fit)[2]
b0eN <- coef(fitN)[1]
X <- cbind(1,x1,x2,x3)
bt <- cbind(c(b0eN,0,0,0),c(b0e,0,0,b3e))
f(X, y, bt)
## [1] -133.7509 -133.1574
-2*(f(X, y, bt)[1]-f(X, y, bt)[2])
## [1] 1.186934
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F)
## [1] 0.2759489
------------------------------------------------------
下記のようにx1の係数だけが0でないとき、ちゃんと尤度比検定でも検出できる（Wald検定でも検出できている）。
------------------------------------------------------
b0 <- 0
b1 <- 0.007
b2 <- 0
b3 <- 0
n <- 100
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
x3 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b3*x3))
fitF <- glm(y ~ x1 + x2 + x3, family = poisson(link = "log")) #full model
fitF2 <- summary(fitF)
fitF2 #Wald検定
##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8971 -1.2968 -0.1198 0.7388 1.7329
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.057824 0.297029 0.195 0.8456
## x1 0.019859 0.010094 1.967 0.0491 *
## x2 -0.017579 0.010689 -1.645 0.1000
## x3 0.006202 0.010884 0.570 0.5688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 115.18 on 99 degrees of freedom
## Residual deviance: 107.59 on 96 degrees of freedom
## AIC: 285.23
##
## Number of Fisher Scoring iterations: 5
Anova(fitF) #尤度比検定
## Analysis of Deviance Table (Type II tests)
##
## Response: y
## LR Chisq Df Pr(>Chisq)
## x1 3.9140 1 0.04789 *
## x2 2.7437 1 0.09764 .
## x3 0.3251 1 0.56856
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
・連続説明変数が2つとそれらの交互作用があるとき
次は、交互作用を含む場合を考えよう。交互作用を含む時は、尤度比の計算を間違えやすいので注意が必要だ。
------------------------------------------------------
b0 <- 0
b1 <- 0
b2 <- 0
b12 <- 0
x1 <- runif(n = n, 0, 30)
x2 <- runif(n = n, 0, 30)
y <- rpois(n = n, lambda = exp(b0 + b1*x1 + b2*x2 + b12*x1*x2))
fitF <- glm(y ~ x1*x2, family = poisson(link = "log"))
fitF2 <- summary(fitF)
fitF2 #Wald検定
##
## Call:
## glm(formula = y ~ x1 * x2, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6635 -1.4253 -0.1001 0.7389 2.7984
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3697290 0.3316054 1.115 0.265
## x1 -0.0193161 0.0206942 -0.933 0.351
## x2 -0.0088707 0.0189731 -0.468 0.640
## x1:x2 0.0005534 0.0011772 0.470 0.638
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 120.40 on 99 degrees of freedom
## Residual deviance: 119.12 on 96 degrees of freedom
## AIC: 280.45
##
## Number of Fisher Scoring iterations: 5
Anova(fitF) #尤度比検定
## Analysis of Deviance Table (Type II tests)
##
## Response: y
## LR Chisq Df Pr(>Chisq)
## x1 1.05580 1 0.3042
## x2 0.02007 1 0.8873
## x1:x2 0.22013 1 0.6389
------------------------------------------------------
まず真っ先に、交互作用項の尤度比を求めよう。以下のように行う。
------------------------------------------------------
fit <- glm(y ~ x1 + x2, family = poisson(link = "log")) #x1:x2に関するreduced model
b0e <- coef(fit)[1]
b1e <- coef(fit)[2]
b2e <- coef(fit)[3]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b12eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x1*x2)
bt <- cbind(c(b0e,b1e,b2e,0),c(b0eF,b1eF,b2eF,b12eF))
f(X, y, bt)
## [1] -136.3344 -136.2244
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #x1:x2に関する-2 log(Λ)
## [1] 0.2201261
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #x1:x2の尤度比検定のp値
## [1] 0.6389439
------------------------------------------------------
交互作用項もパラメータは1つだから、従うカイ二乗分布は自由度1で大丈夫だ。
問題はここからで、x1やx2の尤度比検定はひと手間必要である。例えば、x1の尤度比検定を行うとき、フルモデルは交互作用のないモデル、縮小モデルは交互作用のないモデルのうち、x1を除いたモデルとなる。つまり、交互作用項を残したままでは検定できないということだ。これが、交互作用のないモデルとの違いとなる。
------------------------------------------------------
fitFr <- glm(y ~ x1 + x2, family = poisson(link = "log"))
#交互作用のないモデルを構築、これがx1とx2に関するfull modelとなる。
fit <- glm(y ~ x2, family = poisson(link = "log")) #x1に関するreduced model
b0e <- coef(fit)[1]
b2e <- coef(fit)[2]
b0eF <- coef(fitFr)[1]
b1eF <- coef(fitFr)[2]
b2eF <- coef(fitFr)[3]
X <- cbind(1,x1,x2)
bt <- cbind(c(b0e,0,b2e),c(b0eF,b1eF,b2eF))
f(X, y, bt)
## [1] -136.8623 -136.3344
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #x1に関する-2 log(Λ)
## [1] 1.055804
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #x1の尤度比検定のp値
## [1] 0.3041742
------------------------------------------------------
確かに一致した。x2に関しても同様の手順で求まる。以下に間違いの例を紹介しよう。いずれもフルモデルから交互作用項を除いていないことでうまく解析ができなくなっている。
------------------------------------------------------
#以下は間違い
fit <- glm(y ~ x2 + x1:x2, family = poisson(link = "log"))
b0e <- coef(fit)[1]
b2e <- coef(fit)[2]
b12e <- coef(fit)[3]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b12eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x1*x2)
bt <- cbind(c(b0e,0,b2e,b12e),c(b0eF,b1eF,b2eF,b12eF))
f(X, y, bt)
## [1] -136.6627 -136.2244
-2*(f(X, y, bt)[1]-f(X, y, bt)[2])
## [1] 0.8767295
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F)
## [1] 0.349099
fit <- glm(y ~ x2, family = poisson(link = "log"))
b0e <- coef(fit)[1]
b2e <- coef(fit)[2]
b0eF <- coef(fitF)[1]
b1eF <- coef(fitF)[2]
b2eF <- coef(fitF)[3]
b12eF <- coef(fitF)[4]
X <- cbind(1,x1,x2,x1*x2)
bt <- cbind(c(b0e,0,b2e,0),c(b0eF,b1eF,b2eF,b12eF))
f(X, y, bt)
## [1] -136.8623 -136.2244
-2*(f(X, y, bt)[1]-f(X, y, bt)[2])
## [1] 1.27593
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F)
## [1] 0.2586572
------------------------------------------------------
・カテゴリカル説明変数（2水準と3水準）が2つのとき
最後、2水準と3水準のカテゴリカル変数があるときについて確認しよう。さらに、3水準のうち、1つだけが有意な効果を持つ場合だ。
------------------------------------------------------
b0 <- 0
bB <- 0
bb <- 0
bc <- 0.8
n <- 120
g1 <- rep(LETTERS[1:2], each = n/2)
g2 <- rep(letters[1:3], times = n/3)
X <- model.matrix( ~ g1 + g2)
b <- c(b0,bB,bb,bc)
y <- rpois(n = n, lambda = exp(X%*%b))
fitF <- glm(y ~ g1 + g2, family = poisson(link = "log")) #full model
fitF2 <- summary(fitF)
fitF2 #Wald検定
##
## Call:
## glm(formula = y ~ g1 + g2, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1508 -1.3191 -0.2107 0.4497 2.0161
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2288 0.1600 1.430 0.152647
## g1B -0.0113 0.1503 -0.075 0.940084
## g2b -0.3567 0.2204 -1.618 0.105578
## g2c 0.6098 0.1757 3.471 0.000519 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 177.08 on 119 degrees of freedom
## Residual deviance: 148.43 on 116 degrees of freedom
## AIC: 371.66
##
## Number of Fisher Scoring iterations: 5
Anova(fitF) #尤度比検定
## Analysis of Deviance Table (Type II tests)
##
## Response: y
## LR Chisq Df Pr(>Chisq)
## g1 0.0056 1 0.9401
## g2 28.6375 2 6.046e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
このとき、尤度比検定では、おのおののカテゴリカル変数全体で有意かどうかの判定が行われる。ちょうど、分散分析の結果と同じような感じである。どの水準で有意かどうかを判定するには、さらに別の解析、Wald検定やTukey-Krammer法などが必要になる。ただ、論文によっては尤度比検定で出力される分散分析表のような要約のほうが説明しやすい場合もあるので、便利ではある。
今回の例では、g1には効果がなく、g2には効果があると判定される。では、具体的に推移値を求めていこう。まずはg1について。
------------------------------------------------------
fit <- glm(y ~ g2, family = poisson(link = "log")) #g1に関するreduced model
b0e <- coef(fit)[1]
bbe <- coef(fit)[2]
bce <- coef(fit)[3]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bbeF <- coef(fitF)[3]
bceF <- coef(fitF)[4]
X <- model.matrix( ~ g1 + g2)
bt <- cbind(c(b0e,0,bbe,bce),c(b0eF,bBeF,bbeF,bceF))
f(X, y, bt)
## [1] -181.8336 -181.8308
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #g1に関する-2 log(Λ)
## [1] 0.005649748
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 1, lower.tail = F) #g1の尤度比検定のp値
## [1] 0.9400836
------------------------------------------------------
g1は2水準なので、フルモデルと縮小モデルでの推定するパラメータ数の差は1である。ゆえに、カイ二乗分布は自由度1のものを用いる。確かに一致した。
次にg2について。
------------------------------------------------------
fit <- glm(y ~ g1, family = poisson(link = "log")) #g2に関するreduced model
b0e <- coef(fit)[1]
bBe <- coef(fit)[2]
b0eF <- coef(fitF)[1]
bBeF <- coef(fitF)[2]
bbeF <- coef(fitF)[3]
bceF <- coef(fitF)[4]
X <- model.matrix( ~ g1 + g2)
bt <- cbind(c(b0e,bBe,0,0),c(b0eF,bBeF,bbeF,bceF))
f(X, y, bt)
## [1] -196.1495 -181.8308
-2*(f(X, y, bt)[1]-f(X, y, bt)[2]) #g2に関する-2 log(Λ)
## [1] 28.63751
pchisq(-2*(f(X, y, bt)[1]-f(X, y, bt)[2]), df = 2, lower.tail = F) #g2の尤度比検定のp値
## [1] 6.045666e-07
------------------------------------------------------
g2は3水準なので、フルモデルと縮小モデルでの推定するパラメータ数の差は2である。ゆえに、カイ二乗分布は自由度2のものを用いる。確かに一致した。
以上のように尤度比検定を行うことが可能である。例え、カテゴリカル変数で交互作用がある場合でも、推定するパラメータ数などに注意して、カイ二乗分布を選べば、問題なく解析が可能である。
※2024.05.22 黒木玄先生より2024.5.19にXにて間違いの指摘があったので序文を修正。また直すかも。