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

Рассмотрим вопрос прогноза величины отклика \(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]\).

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