Loading [MathJax]/extensions/tex2jax.js

четверг, 25 октября 2018 г.

Как построить доверительный интервал для прямой регрессии? (Линейная регрессия)

Чтобы определить область которой принадлежит регрессионная прямая \(E[Y]=\beta_0+\beta_1 X\) для последней можно построить доверительный интервал.

\(1-\alpha\) доверительные границы Уоркинга-Хотеллинга для регрессионной прямой в случае простой линейной модели с нормальной ошибкой в точке \(x_k\) имеют вид:\[\hat{y}_h\pm W\cdot s[\hat{y}_h],\]где \(W^2=2F_{1-\alpha}(2,n-2)\), \(\hat{y}_h\) и \(s[\hat{y}_h]\) подробно рассмотрены в статье о построении интервальной оценки для математического ожидания отклика.

Рассмотрим данные:
x <- c(80, 30, 50, 90, 70, 60, 120, 80, 100, 50, 40, 70, 90, 20, 110, 100, 30, 50, 90, 110, 30, 90, 40, 80, 70)
y <- c(399, 121, 221, 376, 361, 224, 546, 352, 353, 157, 160, 252, 389, 113, 435, 420, 212, 268, 377, 421, 273, 468, 244, 342, 323)
Построим \(90\%\) доверительный интервал для регрессионной прямой при \(x=100\)
> new <- data.frame(x = 100)
> fit <- lm(y ~ x)
> CI <- predict(fit,
+               new,
+               se.fit = TRUE,
+               interval = "confidence",
+               level = 0.9)
> W <- sqrt(2 * qf(0.9, length(fit$coefficients), fit$df.residual))
> Band <- cbind(CI$fit[1] - W * CI$se.fit, CI$fit[1] + W * CI$se.fit)
> Band

##          [,1]    [,2]
## [1,] 387.1591 451.613
Таким образом, с доверительной вероятностью \(0.9\) можно заключить, что регрессионная прямая при значении \(x=100\) находится в границах \([387.1591, 451.613]\).

Аналогичным образом можно рассчитать доверительные границы для других значений переменной \(x\) и вместе с регрессионной прямой нанести их на диаграмму рассеивания (рис. 1):
> library(ggplot2)
> CI_ALL <- predict(fit,
+                   se.fit = TRUE,
+                   interval = "confidence",
+                   level = 0.9)
> W <- sqrt(2 * qf(0.9, length(fit$coefficients), fit$df.residual))
> CI_low <- CI_ALL$fit[, 1] - W * CI_ALL$se.fit
> Y_hat <- CI_ALL$fit[, 1]
> CI_up <- CI_ALL$fit[, 1] + W * CI_ALL$se.fit
> b_0 <- fit$coefficients[1]
> b_1 <- fit$coefficients[2]
> mydata <- data.frame(x, y, Y_hat, CI_low, CI_up)
> ggplot(data = mydata) + geom_point(aes(x = x, y = y), size = 1.5) +
+   geom_abline(aes(intercept = b_0,  slope = b_1), size = .5) +
+   geom_errorbar(
+     aes(
+       x = x,
+       y = Y_hat,
+       ymin = CI_low,
+       ymax = CI_up
+     ),
+     color = "red",
+     width = 3,
+     size = 1
+   )
Рис 1
Из рисунка видно, что регрессионная прямая оценена довольно точно, доверительные интервалы становятся шире для крайних (максимальных и минимальных) значений \(x\), регрессионная прямая имеет явный положительный наклон.

Комментариев нет:

Отправить комментарий