Rで確率･統計:線形回帰における回帰係数の信頼区間、母回帰式の区間推定、目的変数の予測区間、回帰直線の同時信頼区間

Rでデータサイエンス

線形回帰における回帰係数の信頼区間、母回帰式の区間推定、目的変数の予測区間、回帰直線の同時信頼区間

共通準備

  1. 参考引用資料 1,2,3,7
  2. 次の1次1変数式を母回帰式とし、有意水準は全て5%とする

\[y=\alpha x + \beta + \epsilon\]

library(ggplot2)
library(tidyr)
library(dplyr)
library(grid)
# サンプルを作成
# サイズと説明変数の設定
n <- 50
x <- seq(n)
# 誤差の設定
# 今回設定した標準偏差の数値に意味はありません。
# 母回帰式との違いや信頼区間が明瞭となるよう、傾きに対して標準偏差を大きめにしています。
# よって、傾き及び切片が有意とならない場合があります。
error <- rnorm(n = n, mean = 0, sd = 15)
# 目的変数の設定
# 今回設定した傾きと切片の数値に意味はありません。
Y <- 2 * x + 10
y <- Y + error
# サンプルを利用して散布図作成そして説明変数(x)と目的変数(y)の線形回帰を取ります
# 有意水準
significant_level <- 0.05
# 自由度n-2、有意水準0.05のt値
t_value <- qt(1 - significant_level/2, n - 2)
# 説明変数xの偏差平方和
Sxx <- sum((x - mean(x))^2)
# 説明変数xと目的変数yの積和
Sxy <- sum((x - mean(x)) * (y - mean(y)))
# 回帰直線の傾きの推定量
a <- Sxy/Sxx
# 回帰直線の切片の推定量
b <- mean(y) - mean(x) * a
# 推定されたy
estimated_y <- a * x + b
# 残差
residuals <- y - estimated_y
# 残差標準偏差
sigma <- sqrt(sum((residuals - mean(residuals))^2)/(n - 2))
# 土台とする散布図作成
g1 <- ggplot(mapping = aes(x = x)) + geom_point(mapping = aes(y = y), size = 1) + theme_minimal() + theme(legend.position = "top", legend.title = element_blank(), legend.key.width = unit(1, "cm"))
# 設定した傾きと切片に基づく直線を引きます。
g2 <- g1 + geom_line(mapping = aes(x = x, y = Y, col = "母回帰式"), size = 0.7)
# 推定された傾きと切片に基づく直線を引きます。
g <- g2 + geom_line(mapping = aes(x = x, y = a * x + b, col = "推定回帰式"), size = 0.7) + scale_color_manual(values = c(推定回帰式 = "blue", 母回帰式 = "red"))
g
Figure 1 
# 推定された係数(傾き)と切片
list(a = a, b = b)
$a
[1] 2.036153

$b
[1] 8.158429
# なお推定された係数と切片は以下の方法でも取り出せます。
result_lm <- lm(y ~ x)
# result_lm$coefficients[2] # a
# result_lm$coefficients[1] # b
result_lm %>%
    summary()

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.400  -7.844  -0.032   6.578  47.665 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   8.1584     3.9687   2.056              0.0453 *  
x             2.0362     0.1355  15.032 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.82 on 48 degrees of freedom
Multiple R-squared:  0.8248,    Adjusted R-squared:  0.8212 
F-statistic:   226 on 1 and 48 DF,  p-value: < 0.00000000000000022

回帰係数の信頼区間

  1. 推定された係数 alpha と 切片 beta それぞれの信頼区間
  2. 参考資料 3,7
  3. 次の数式が信頼区間を表します

\[\hat{\alpha}\pm t_{\left(n-2,\dfrac{0.05}{2}\right)}\times\dfrac{\hat{\sigma}}{\displaystyle\sqrt{S_{xx}}},\quad \hat{\beta}\pm t_{\left(n-2,\dfrac{0.05}{2}\right)}\times\hat{\sigma}\times{\displaystyle\sqrt{\dfrac{1}{n}+\dfrac{\bar{x}^{2}}{S_{xx}}}}\]

# 係数(傾き)
ci <- t_value * sigma/sqrt(Sxx)
list(lower = a - ci, upper = a + ci)
$lower
[1] 1.763812

$upper
[1] 2.308493
# 切片
ci <- t_value * sigma * sqrt(1/n + mean(x)^2/Sxx)
list(lower = b - ci, upper = b + ci)
$lower
[1] 0.1787964

$upper
[1] 16.13806
# なお関数confint {stats}でも95%信頼区間が直接求められます
confint(result_lm, level = 0.95)
                2.5 %    97.5 %
(Intercept) 0.1787964 16.138061
x           1.7638118  2.308493

母回帰式の区間推定

  1. 母回帰式の存在範囲推定 = 母回帰式の区間推定
  2. なおggplotでse = Tとした場合に表示される信頼区間は、この母回帰式の区間推定になります。
  3. 参考資料 4
  4. 次の数式が区間推定を表します。

\[\left(\hat{\alpha}x_{0}+\hat{\beta}\right)\pm t_{\left(n-2,\dfrac{0.05}{2}\right)}\times\displaystyle\sqrt{\dfrac{1}{n}+\dfrac{\left(x-\bar{x}\right)^{2}}{S_{xx}}}\times \hat{\sigma}\]

ci <- t_value * sqrt(1/n + ((x - mean(x))^2)/Sxx) * sigma
lower1 <- estimated_y - ci
upper1 <- estimated_y + ci
g3 <- g2 + geom_line(mapping = aes(y = estimated_y, col = "推定回帰式"), size = 0.7) + scale_color_manual(values = c(推定回帰式 = "blue", 母回帰式 = "red"))
gg1 <- g3 + geom_line(mapping = aes(y = upper1), size = 0.7, col = "blue", linetype = "dashed") + geom_line(mapping = aes(y = lower1), size = 0.7, col = "blue", linetype = "dashed")
# geom_smoothで回帰直線と信頼区間を描きますと重なります。
gg2 <- gg1 + geom_smooth(mapping = aes(x = x, y = y), method = "lm", se = T, inherit.aes = T)
gridExtra::grid.arrange(gg1, gg2, layout_matrix = c(1, 2) %>%
    matrix(nrow = 1), top = textGrob("母回帰式の区間推定", gp = gpar(font = 1)))
Figure 2

目的変数の予測区間

  1. 説明変数 x の値をx0 に指定したときの目的変数 y の信頼区間 = 予測区間
  2. 参考資料 4
  3. 次の数式が予測区間を表します。

\[\left(\hat{\alpha}x_{0}+\hat{\beta}\right)\pm t_{\left(n-2,\dfrac{0.05}{2}\right)}\times\displaystyle\sqrt{1+\dfrac{1}{n}+\dfrac{\left(x-\bar{x}\right)^{2}}{S_{xx}}}\times \hat{\sigma}\]

# グレーの範囲が予測区間です。
ci <- t_value * sqrt(1 + 1/n + ((x - mean(x))^2)/Sxx) * sigma
lower2 <- estimated_y - ci
upper2 <- estimated_y + ci
g3 <- g2 + geom_line(mapping = aes(y = estimated_y, col = "推定回帰式"), size = 0.7) + scale_color_manual(values = c(推定回帰式 = "blue", 母回帰式 = "red"))
g <- g3 + geom_line(mapping = aes(y = upper2), size = 0.7, col = "blue", linetype = "dashed") + geom_line(mapping = aes(y = lower2), size = 0.7, col = "blue", linetype = "dashed") + geom_ribbon(mapping = aes(ymin = lower2, ymax = upper2), alpha = 0.2) + labs(title = "目的変数 y の信頼区間 = 予測区間")
g
Figure 3

回帰直線の同時信頼区間

  1. x = x0 だけでなく、すべての x に対する回帰直線の同時信頼区間
  2. 参考資料 5,6
  3. 次の数式が同時信頼区間を表します。

\[\left(\hat{\alpha}x_{0}+\hat{\beta}\right)\pm \displaystyle\sqrt{2F_{\left(2,n-2,\dfrac{0.05}{2}\right)}}\times\displaystyle\sqrt{\dfrac{1}{n}+\dfrac{\left(x-\bar{x}\right)^{2}}{S_{xx}}}\times \hat{\sigma}\]

# グレーの範囲が同時信頼区間です。
ci <- sqrt(2 * qf(significant_level, 2, n - 2)) * sqrt(1/n + ((x - mean(x))^2)/Sxx) * sigma
lower3 <- estimated_y - ci
upper3 <- estimated_y + ci
g3 <- g2 + geom_line(mapping = aes(y = estimated_y, col = "推定回帰式"), size = 0.7) + scale_color_manual(values = c(推定回帰式 = "blue", 母回帰式 = "red"))
g <- g3 + geom_line(mapping = aes(y = upper3), size = 0.7, col = "blue", linetype = "dashed") + geom_line(mapping = aes(y = lower3), size = 0.7, col = "blue", linetype = "dashed") + geom_ribbon(mapping = aes(ymin = lower3, ymax = upper3), alpha = 0.2) + labs(title = "すべての x に対する回帰直線の同時信頼区間")
g
Figure 4

参考引用資料

  1. http://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/BS704_Confidence_Intervals/BS704_Confidence_Intervals_print.html
  2. http://lbm.ab.a.u-tokyo.ac.jp/~omori/kokusai11/kokusai11_1212.html
  3. http://www2.econ.osaka-u.ac.jp/~tanizaki/class/2016/basic_econome/03.pdf
  4. http://www.hil.hiroshima-u.ac.jp/sys2/a/kaikibunseki.pdf
  5. http://stat.sys.i.kyoto-u.ac.jp/titech/class/doc/lec20021114-9.pdf
  6. https://www.jstage.jst.go.jp/article/jscswabun/14/1/14_KJ00001896557/_pdf
  7. https://www2.isye.gatech.edu/~yxie77/isye2028/lecture12.pdf

最終更新

Sys.time()
[1] "2024-04-16 05:19:16 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.542'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者