［R］回帰モデルの信頼区間と予測区間

lm関数により求めた回帰モデルの信頼区間と予測区間はpredict関数で求めることができる。以下の計算により、赤実線が求めた回帰直線。桃破線がその95%信頼区間、橙点線がその95%予測区間。

> # 説明変数xと目的変数y
> x <- c(1, 2, 3, 5, 6, 7, 8, 10, 11, 12)
> y <- c(0, 1, 1, 2, 4, 4, 6,  5,  8,  8)
> # 回帰直線を求める
> r <- lm(y ~ x)
> print(r)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x  
    -0.8567       0.7318  
> # yの推定値
> yest <- fitted(r)
> # 信頼区間と予測区間の計算(後続の手順のためdata.frame化)
> dtf <- data.frame(x)
> rcon <- predict(r, dtf, interval = "confidence", level = 0.95)
> rpre <- predict(r, dtf, interval = "prediction", level = 0.95)
> rcon <- data.frame(rcon)
> rpre <- data.frame(rpre)
> # 図の作成(黒丸:観測値,赤実線:回帰直線,桃破線:信頼区間,橙点線:予測区間)
> plot(x, y, type = "n", asp = 1.0)
> lines(x, rpre$lwr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rpre$upr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rcon$lwr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, rcon$upr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, yest, col = "red", lwd = 2.0)
> points(x, y, pch = 20, col = "black")

> # 以下はpredict関数と手計算の計算結果の比較
> ah <- 0.95
> n <- length(x)
> k <- 2
> mxxx <- matrix(c(rep(1, n), x), ncol = 2)
> mxy <- matrix(y, ncol = 1)
> # 最小二乗推定量
> mxb <- solve(t(mxxx) %*% mxxx) %*% (t(mxxx) %*% mxy)
> print(mxb)
           [,1]
[1,] -0.8567050
[2,]  0.7318008
> # yの推定値
> yest <- mxxx %*% mxb
> # 残差分散
> s <- sqrt(sum((y - yest) ^ 2) / (n - k))
> #
> s0 <- sh0 <- double(n)
> for (i in 1:n) {
+ mxx0 <- matrix(c(1, x[i]), ncol = 1)
+ s0[i] <- s * sqrt(t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ sh0[i] <- s * sqrt(1 + t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ }
> # 95%信頼区間
> ycon1 <- yest - s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ycon2 <- yest + s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # 95%予測区間
> ypre1 <- yest - sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ypre2 <- yest + sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # predict関数による信頼区間と予測区間
> print(data.frame(rcon$lwr, rcon$upr, rpre$lwr, rpre$upr))
     rcon.lwr  rcon.upr   rpre.lwr  rpre.upr
1  -1.1763909 0.9265825 -2.2315171  1.981709
2  -0.3152118 1.5290049 -1.4382140  2.652007
3   0.5349489 2.1424457 -0.6558465  3.333241
4   2.1772621 3.4273356  0.8728263  4.731771
5   2.9513451 4.1168542  1.6179064  5.450293
6   3.6831458 4.8486549  2.3497072  6.182094
7   4.3726644 5.6227379  3.0682286  6.927174
8   5.6575543 7.2650511  4.4667589  8.455846
9   6.2709951 8.1152118  5.1479929  9.238214
10  6.8734175 8.9763909  5.8182913 10.031517
> # 手計算による信頼区間と予測区間
> print(data.frame(ycon1, ycon2, ypre1, ypre2))
        ycon1     ycon2      ypre1     ypre2
1  -1.1763909 0.9265825 -2.2315171  1.981709
2  -0.3152118 1.5290049 -1.4382140  2.652007
3   0.5349489 2.1424457 -0.6558465  3.333241
4   2.1772621 3.4273356  0.8728263  4.731771
5   2.9513451 4.1168542  1.6179064  5.450293
6   3.6831458 4.8486549  2.3497072  6.182094
7   4.3726644 5.6227379  3.0682286  6.927174
8   5.6575543 7.2650511  4.4667589  8.455846
9   6.2709951 8.1152118  5.1479929  9.238214
10  6.8734175 8.9763909  5.8182913 10.031517