［R］回帰分析における対数尤度を簡単に求める

回帰分析をlm関数で行い、その戻り値を用いてlogLik関数を使う。以下は「情報量規準」（朝倉書店）に掲載の計算例（pp.59-61）を再現した結果。

観測値はpp.59-60に掲載されており、以下のとおり。これをtable_p059.csvと保存する。

 i,    x,     y
 1, 0.00, 0.854
 2, 0.05, 0.786
 3, 0.10, 0.706
 4, 0.15, 0.763
 5, 0.20, 0.772
 6, 0.25, 0.693
 7, 0.30, 0.805
 8, 0.35, 0.739
 9, 0.40, 0.760
10, 0.45, 0.764
11, 0.50, 0.810
12, 0.55, 0.791
13, 0.60, 0.798
14, 0.65, 0.841
15, 0.70, 0.882
16, 0.75, 0.879
17, 0.80, 0.863
18, 0.85, 0.934
19, 0.90, 0.971
20, 0.95, 0.985

以下、計算結果。当該書籍では、上記の観測値を用いた多項式回帰モデルの残差分散と対数尤度の一覧を表にまとめて掲載している（p.61）。

> dtf <- read.csv("table_p059.csv", header = TRUE)
> n <- nrow(dtf)
> for (k in 1:9) {
+     r <- lm(y ~ poly(x, k), data = dtf)
+     rss <- sum(r$re ^ 2)
+     cat(sprintf("次数:%d, 残差分散:%8.6f, 対数尤度:%5.2f\n", k, rss / n, logLik(r)))
+ }
次数:1, 残差分散:0.002587, 対数尤度:31.19
次数:2, 残差分散:0.000922, 対数尤度:41.51
次数:3, 残差分散:0.000833, 対数尤度:42.52
次数:4, 残差分散:0.000737, 対数尤度:43.75
次数:5, 残差分散:0.000688, 対数尤度:44.44
次数:6, 残差分散:0.000650, 対数尤度:45.00
次数:7, 残差分散:0.000622, 対数尤度:45.44
次数:8, 残差分散:0.000607, 対数尤度:45.69
次数:9, 残差分散:0.000599, 対数尤度:45.83