尤度比検定

Rのlmtestパッケージにあるlrtest関数で尤度比検定をやってみます。

パルマーペンギンのデータを用います。なんでも、R 4.5.0からデフォルトのデータセットに組み込まれるとのことです。

library(palmerpenguins)
library(ggplot2)
library(lmtest)
## 要求されたパッケージ zoo をロード中です
## 
## 次のパッケージを付け加えます: 'zoo'
## 以下のオブジェクトは 'package:base' からマスクされています:
## 
##     as.Date, as.Date.numeric

アデリーペンギンのデータを抜き出します。欠測は完全にランダムに発生していると仮定して、欠測値のある行は除きます。

生息地の島ごとに色を変えて嘴高と嘴長をプロットします。

Adelie <- penguins |>
  dplyr::filter(species == "Adelie",
                !is.na(penguins$bill_depth_mm),
                !is.na(penguins$bill_length_mm))
ggplot(Adelie) +
  geom_point(aes(x = bill_depth_mm, y = bill_length_mm,
                 colour = island))

モデル

嘴長を目的変数とします。モデル1は説明変数が嘴高のみの線形回帰のモデルです。モデル2はさらに島の因子を加え、島により線形回帰の切片が異なるとする共分散分析のモデルです。

モデル1: bill_length_mm ~ bill_depth_mm
モデル2: bill_length_mm ~ bill_depth_mm + island

model1 <- lm(bill_length_mm ~ bill_depth_mm,
             data = Adelie)
model2 <- lm(bill_length_mm ~ bill_depth_mm + island,
             data = Adelie)

lrtest関数で尤度比検定をおこないます。

lrtest(model1, model2)
## Likelihood ratio test
## 
## Model 1: bill_length_mm ~ bill_depth_mm
## Model 2: bill_length_mm ~ bill_depth_mm + island
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   3 -349.12                     
## 2   5 -348.78  2 0.6744     0.7138

p > 0.05 なので、帰無仮説は棄却できない、すなわち島によって切片が異なるとはいえない、ということになりました。

baseの関数だけでおこなう場合

baseの関数だけでも以下のようにしてできます。

モデル1の対数尤度を求めます。

(logLik1 <- logLik(model1))
## 'log Lik.' -349.1176 (df=3)

モデル2の対数尤度を求めます。

(logLik2 <- logLik(model2))
## 'log Lik.' -348.7804 (df=5)

両モデルの自由度と対数尤度の差×-2を求めます。

df1 <- attr(logLik1, "df")
df2 <- attr(logLik2, "df")
(delta <- -2 * (as.numeric(logLik1) - as.numeric(logLik2)))
## [1] 0.6744218

両モデルの自由度の差を自由度とし、対数尤度の差×-2を統計量としてカイ2乗分布の確率を求めます。

pchisq(delta, df = df2 - df1, lower.tail = FALSE)
## [1] 0.7137583

lrtestと同じ値が得られました。