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

Рассмотрим вопрос прогноза величины отклика $y_{h(new)}$ для заданного значения $x_h$ переменного $x$. Значение $y_{h(new)}$ следует рассматривать как результат нового испытания, независимого от тех испытаний по которым построена регрессионная модель. Следует также предположить, что условия в которых построена модель остаются неизменными при расчета величины $y_{h(new)}$.

Отметим также разницу между двумя оценками $\hat{y}_h$ и $y_{h(new)}$. В первом случае оценивается математическое ожидание распределения $Y_h$, во втором случае оценивается индивидуальный исход, извлеченный из этого распределения.

Чтобы раскрыть природу предсказательного интервала для нового наблюдения $y_{h(new)}$ предположим вначале, что все параметры регрессии известны: $\beta_0=0.10$, $\beta_1=0.95$, $E[Y]=0.10+0.95x$, $\sigma=0.12$ и $x_h=3.5$. 

Тогда в случае простой линейной регрессионной модели с нормальной ошибкой случайная величина $Y_h$ будет иметь нормальное распределение с математическим ожиданием $E(Y_h)=0.10+0.95(3.5)=3.425$ и стандартным отклонением $\sigma=0.12$. Предположим, что требуется построить прогноз для отклика в пределах $E[Y_h]\pm3\sigma$. В нашем случае получим $3.425\pm3(0.12)=[3.065, 3.785]$. Так как вероятность попасть в симметричный относительно среднего интервал $\pm3\sigma$ для нормального распределения равна $0.9973$, то с такой же вероятностью рассчитанный выше интервал дает верный прогноз для отклика при $x_h=3.5$.

Таким образом, идея построения предсказательного интервала состоит в том, чтобы выбрать ту область значений распределения $Y$ которой принадлежит большая часть наблюдений и, следовательно, можно утверждать, что новое наблюдение также попадет в данный промежуток.

В общем случае, когда параметры простой линейной регрессионной модели с нормальной ошибкой известны, $1-\alpha$ предсказательный интервал для $y_{h(new)}$ имеет вид: $$E[Y_h]\pm z_{1-\alpha/2}\sigma.$$ Центрирование границ относительно $E[Y_h]$ позволяет получить наиболее узкий интервал при заданной вероятности.

Когда параметры регрессии не известны, они должны быть оценены: вместо математического ожидания $E[Y_h]$ естественно принять оценку $\hat{y}_h$, вместо дисперсии $var[Y_i]$ оценку MSE. Однако, точное значение $E[Y_h]$ не известно и лишь оценивается с помощью доверительного интервала. Поэтому при построении предсказательного интервала $y_{h(new)}$ следует учесть неопределенность связанную с центром распределения $Y_h$ и неопределенность связанную с вариацией непосредственно $Y_h$.

Предсказательный интервал для нового наблюдения $y_{h(new)}$ в точке $x_h$ может быть получен из следующей теоремы: стандартизированная статистика $$\frac{y_{h(new)}-\hat{y}_h}{s_{pred}}$$ имеет $t$-распределение с $n-2$ степенями свободы. Откуда $1-\alpha$ предсказательный интервал для нового наблюдения $y_{h(new)}$ имеет вид: $$\hat{y}_h\pm t_{1-\alpha/2}(n-2)\cdot s_{pred}.$$

Заметим, что числитель в рассматриваемой статистике показывает насколько новое наблюдение $y_{h(new)}$ отклоняется от оценки $\hat{y}_h$ неизвестного математического ожидания $E[Y_h]$, рассчитанной по включенным в модель $n$ наблюдениям. Эту разницу можно рассматривать в качестве ошибки предсказания, если считать $\hat{y}_h$ лучшей оценкой для $y_{h(new)}$. 

Дисперсию ошибки предсказания можно получить используя независимость нового наблюдения $y_{h(new)}$ и рассчитанной по $n$ наблюдениям оценки $\hat{y}_h$: $$\sigma^2_{pred}=var[y_{h(new)}-\hat{y}_h]=var[y_{h(new)}]+var[\hat{y}_h]=\sigma^2+var[\hat{y}_h].$$ Заметим, что дисперсия $\sigma^2_{pred}$ состоит из двух компонент: дисперсии $\sigma^2$ случайной величины $Y$ в точке $x_h$ и дисперсии $var[\hat{y}_h]$ выборочного распределения $\hat{y}_h$. Несмещенной оценкой $\sigma^2_{pred}$ является: $$s^2_{pred}=MSE+s^2[\hat{y}_h]=MSE\left[1+\frac{1}{n}+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right].$$

Рассмотрим данные:

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_1=100$

> new <- data.frame(x = 100)

> fit <- lm(y ~ x)

> predict(fit,

+         new,

+         se.fit = FALSE,

+         interval = "prediction",

+         level = 0.9)

##      fit      lwr      upr

## 419.3861 332.2072 506.5649

Таким образом, с доверительной вероятностью 0.9 можно предсказать, что значение отклика в точке $x=100$ находится в границах $[332.2072, 506.5649]$.

Обычно более широкий предсказательный интервал используется для контроля соответствия исследуемого процесса исходным параметрам.

Иногда требуется построить предсказательный интервал для среднего арифметического $\bar{y}_{h(new)}$ из $m$ откликов при заданном $x_h$ переменной $x$. 

Можно показать, что в предположении независимости новых наблюдений $Y$, соответствующие $1-\alpha$ границы имеют вид: $$\hat{y}_h\pm t_{1-\alpha/2}(n-2)\cdot s_{pred-mean},$$ где $$s^2_{pred-mean}=\frac{MSE}{m}+s^2[\hat{y}_h]=MSE\left[\frac{1}{m}+\frac{1}{n}+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right].$$ Заметим, что оценка $s^2_{pred-mean}$ состоит из двух компонент: дисперсии среднего арифметического $m$ наблюдений из распределения $Y$ в точке $x_h$ и дисперсии $var[\hat{y}_h]$ выборочного распределения $\hat{y}_h$.

По тем же данным построим $90\%$ предсказательный интервал для среднего арифметического трех независимых откликов при $x_1=100$

> predict(fit,

+         new,

+         se.fit = TRUE,

+         interval = "confidence",

+         level = 0.9);

## $`fit`

##      fit      lwr     upr

## 419.3861 394.9251 443.847

## 

## $se.fit

## [1] 14.27233

## 

## $df

## [1] 23

## 

## $residual.scale

## [1] 48.82331

> MSE <- 48.82331^2;

> s <- sqrt(MSE/3+14.27233^2);

> qt <- qt(p = 0.95, df = 23);

> c(419.3861-s*qt,419.3861+s*qt)

## [1] 365.2356 473.5366

Таким образом, с доверительной вероятностью $0.9$ можно предсказать, что среднее арифметическое значение трех отклика в точке $x=100$ находится в границах $[365.2356, 473.5366]$.

Отметим, что полученный предсказательный интервал существенно уже рассчитанного выше аналогичного предсказательного интервала для единственного значения отклика.