［R］例3.3 配達時間 (1),(2)（「線形回帰分析」、pp.106-113）

式(3.62)のモデルで計算している。表3.2（p.112）の対数尤度などは完全に一致しない。書籍では、途中で値を切り捨てて計算をしているからだと思われる。

> dtf <- read.csv("table3_1.csv", header = TRUE)
> n <- nrow(dtf)
> xx2i <- as.double(dtf$CASE)
> xx3i <- as.double(dtf$DIS) ^ 2 / 1.e5
> yyi <- as.double(dtf$DVT)
> k <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1., n), xx2i, xx3i), ncol = 3)
> mxbh <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxe <- mxy - mxyh
> ssy2 <- sum(t(mxy) %*% mxy) - sum(yyi) ^ 2 / n
> ssyh2 <- sum(t(mxyh) %*% mxyh) - sum(yyi) ^ 2 / n
> sse2 <- as.vector(t(mxe) %*% mxe)
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> mxvv <- s2 * solve(t(mxxx) %*% mxxx)
> mxs <- sqrt(diag(mxvv))
> rr2 <- ssyh2 / ssy2
> rrb2 <- 1 - (n - 1) / (n - k) * (1 - rr2)
> tbh <- as.vector(mxbh / mxs)
> sh2 <- sse2 / n
> logll <- -n / 2 * log(2 * pi) - n / 2 * log(sh2) - sse2 / (2 * sh2)
> p <- k + 1
> aic <- -2 * logll + 2 * p
> sbic <- -2 * logll + p * log(n)
> gcv <- -2 * logll - 2 * n * log(1 - p / n)
> hq <- -2 * logll + 2 * p * log(log(n))
> yyhi <- as.vector(mxyh)
> ei <- yyi - yyhi
> ri <- ei / yyi * 100
> ai2 <- ei ^ 2 / sum(ei ^ 2) * 100
> #
> cat("求めた回帰係数 β\n")
求めた回帰係数 β
> print(mxbh)
         [,1]
[1,] 6.191398
[2,] 1.410930
[3,] 1.424706
> cat("回帰係数の標準誤差 s\n")
回帰係数の標準誤差 s
> print(mxs)
[1] 0.8673945 0.1276345 0.1992415
> cat("仮説を検定するためのt値\n")
仮説を検定するためのt値
> print(tbh)
[1]  7.137927 11.054459  7.150648
> cat(sprintf("R^2   = %8.4f\n", rr2))
R^2   =   0.9791
> cat(sprintf("R~^2  = %8.4f\n", rrb2))
R~^2  =   0.9772
> cat(sprintf("log L = %8.4f\n", logll))
log L = -55.1821
> cat(sprintf("AIC   = %8.4f\n", aic))
AIC   = 118.3641
> cat(sprintf("SBIC  = %8.4f\n", sbic))
SBIC  = 123.2396
> cat(sprintf("GCV   = %8.4f\n", gcv))
GCV   = 119.0818
> cat(sprintf("HQ    = %8.4f\n", hq))
HQ    = 119.7164
> st <- sprintf("%5.2f %5.2f %5.2f %6.2f %5.2f", yyi, yyhi, ei, ri, ai2)
> for (i in 1:n) cat(sprintf("%2d %s\n", i, st[i]))
 1 16.68 20.54 -3.86 -23.12 12.29
 2 11.50 11.11  0.39   3.36  0.12
 3 12.03 12.07 -0.04  -0.34  0.00
 4 14.88 11.93  2.95  19.85  7.21
 5 13.75 14.98 -1.23  -8.93  1.25
 6 18.11 17.62  0.49   2.71  0.20
 7  8.00  9.19 -1.19 -14.82  1.16
 8 17.83 16.70  1.13   6.36  1.06
 9 79.24 78.89  0.35   0.44  0.10
10 21.50 18.46  3.04  14.14  7.64
11 40.33 35.51  4.82  11.95 19.20
12 21.00 20.96  0.04   0.19  0.00
13 13.50 12.76  0.74   5.47  0.45
14 19.75 17.70  2.05  10.39  3.48
15 24.00 21.75  2.25   9.38  4.19
16 29.00 28.88  0.12   0.41  0.01
17 15.35 15.23  0.12   0.80  0.01
18 19.00 16.32  2.68  14.13  5.95
19  9.50 10.44 -0.94  -9.92  0.73
20 35.10 38.62 -3.52 -10.04 10.27
21 17.90 20.58 -2.68 -14.97  5.94
22 52.32 52.22  0.10   0.19  0.01
23 18.75 21.77 -3.02 -16.13  7.56
24 19.83 23.22 -3.39 -17.11  9.52
25 10.75 12.16 -1.41 -13.08  1.63