Как функция expt.lm () вычисляет доверительный интервал и интервал прогнозирования?

Я провел регрессию:

CopierDataRegression <- lm(V1~V2, data=CopierData1)

и моей задачей было получить

90%доверительный интервал для среднего ответаV2=6 а также90%интервал прогнозирования когдаV2=6.

Я использовал следующий код:

X6 <- data.frame(V2=6)
predict(CopierDataRegression, X6, se.fit=TRUE, interval="confidence", level=0.90)
predict(CopierDataRegression, X6, se.fit=TRUE, interval="prediction", level=0.90)

и я получил(87.3, 91.9) а также(74.5, 104.8) что кажется правильным, так как ИП должен быть шире.

Выход для обоих также включенse.fit = 1.39 который был таким же.Я не понимаю, что это за стандартная ошибка. Разве стандартная ошибка не должна быть больше для PI против CI? Как мне найти эти две разные стандартные ошибки в R?

Данные:

CopierData1 <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L, 
          4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L, 
          66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L, 
          90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L, 
          61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L, 
          10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L, 
          2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L, 
          2L, 4L, 5L)), .Names = c("V1", "V2"),
          class = "data.frame", row.names = c(NA, -45L))
 Gregor29 июн. 2016 г., 22:50
Смотря на?predict.lm, это говорит:"se.fit: стандартная ошибка прогнозируемых средств », «Предсказанное означает» звучит так, как будто оно относится только к доверительному интервалу. Если вы не хотите его видеть, установитеse.fit = FALSE.
 Mitty29 июн. 2016 г., 22:58
Спасибо. Я думаю, что я спрашиваю, как я могу вычислить две ошибки STD на картинке? Так что я могу проверить вычисления и узнать, как они получены.

Ответы на вопрос(2)

есть ли быстрый способ извлечь стандартную ошибку для интервала прогнозирования, но вы всегда можете отложить интервалы для SE (даже если это не супер элегантный подход):

m <- lm(V1 ~ V2, data = d)                                                                                                                                                                                                                

newdat <- data.frame(V2=6)                                                                                                                                                                                                                
tcrit <- qt(0.95, m$df.residual)                                                                                                                                                                                                          

a <- predict(m, newdat, interval="confidence", level=0.90)                                                                                                                                                                                
cat("CI SE", (a[1, "upr"] - a[1, "fit"]) / tcrit, "\n")                                                                                                                                                                                   

b <- predict(m, newdat, interval="prediction", level=0.90)                                                                                                                                                                                
cat("PI SE", (b[1, "upr"] - b[1, "fit"]) / tcrit, "\n") 

Обратите внимание, что CI SE является тем же значением изse.fit.

 Mitty30 июн. 2016 г., 18:49
Это сработало. Я вернулся к SE, используя 89,63 + - t (0,95,43) xSE = нижняя граница, где нижняя граница была 87,28 для CI и 74,46 для PI. SE CI был 1,39, а SE PI - 9,02. Таким образом, SE для интервала прогнозирования больше доверительного интервала. Но я все еще не понимаю, почему выходные данные в R для интервала предсказания перечисляют se.fit = 1,39. Почему в списке нет 9? Спасибо!!!
Решение Вопроса

interval а такжеlevel аргумент,predict.lm может вернуть доверительный интервал (CI) или интервал прогнозирования (PI). Этот ответ показывает, как получить CI и PI без установки этих аргументов. Есть два способа:

использовать результат средней стадии изpredict.lm;делать все с нуля.

Знание того, как работать с обоими способами, даст вам полное понимание процедуры прогнозирования.

Обратите внимание, что мы рассмотрим толькоtype = "response" (по умолчанию) чехол дляpredict.lm, Обсуждениеtype = "terms" выходит за рамки этого ответа.

Настроить

Я собираю ваш код здесь, чтобы помочь другим читателям копировать, вставлять и запускать. Я также изменяю имена переменных, чтобы они имели более четкое значение. Кроме того, я расширяюnewdat включить более одной строки, чтобы показать, что наши вычисления "векторизованы".

dat <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L, 
          4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L, 
          66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L, 
          90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L, 
          61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L, 
          10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L, 
          2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L, 
          2L, 4L, 5L)), .Names = c("V1", "V2"),
          class = "data.frame", row.names = c(NA, -45L))

lmObject <- lm(V1 ~ V2, data = dat)

newdat <- data.frame(V2 = c(6, 7))

Ниже приведены результатыpredict.lm, чтобы сравнить с нашими ручными вычислениями позже.

predict(lmObject, newdat, se.fit = TRUE, interval = "confidence", level = 0.90)
#$fit
#        fit       lwr      upr
#1  89.63133  87.28387  91.9788
#2 104.66658 101.95686 107.3763
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508

predict(lmObject, newdat, se.fit = TRUE, interval = "prediction", level = 0.90)
#$fit
#        fit      lwr      upr
#1  89.63133 74.46433 104.7983
#2 104.66658 89.43930 119.8939
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508
Используйте промежуточный результат изpredict.lm
## use `se.fit = TRUE`
z <- predict(lmObject, newdat, se.fit = TRUE)
#$fit
#        1         2 
# 89.63133 104.66658 
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508

Что такоеse.fit?

z$se.fit стандартная ошибка предсказанного среднегоz$fitиспользуется для построения КИ дляz$fit, Нам также нужны квантили t-распределения со степенью свободыz$df.

alpha <- 0.90  ## 90%
Qt <- c(-1, 1) * qt((1 - alpha) / 2, z$df, lower.tail = FALSE)
#[1] -1.681071  1.681071

## 90% confidence interval
CI <- z$fit + outer(z$se.fit, Qt)
colnames(CI) <- c("lwr", "upr")
CI
#        lwr      upr
#1  87.28387  91.9788
#2 101.95686 107.3763

Мы видим, что это согласуется сpredict.lm(, interval = "confidence").

Какая стандартная ошибка для PI?

PI шире, чем CI, так как он учитывает остаточную дисперсию:

variance_of_PI = variance_of_CI + variance_of_residual

Обратите внимание, что это определено поэтапно. Для невзвешенной линейной регрессии (как в вашем примере) остаточная дисперсия везде одинакова (известна какгомоскедастичность), и этоz$residual.scale ^ 2, Таким образом, стандартная ошибка для PI

se.PI <- sqrt(z$se.fit ^ 2 + z$residual.scale ^ 2)
#       1        2 
#9.022228 9.058082 

и PI построен как

PI <- z$fit + outer(se.PI, Qt)
colnames(PI) <- c("lwr", "upr")
PI
#       lwr      upr
#1 74.46433 104.7983
#2 89.43930 119.8939

Мы видим, что это согласуется сpredict.lm(, interval = "prediction").

замечание

Все будет сложнее, если у вас есть весовая линейная регрессия, где остаточная дисперсия не везде одинакова, так чтоz$residual.scale ^ 2 должны быть взвешены. Проще построить PI для подгоночных значений (то есть вы не устанавливаетеnewdata когда используешьtype = "prediction" вpredict.lm), потому что вес известен (вы должны были предоставить его черезweight аргумент при использованииlm). Для прогнозирования вне выборки (то есть вы передаетеnewdata вpredict.lm),predict.lm ожидает, что вы скажете ему, как следует взвешивать остаточную дисперсию. Вам нужно либо использовать аргументpred.var или жеweights вpredict.lmв противном случае вы получите предупреждение отpredict.lm жалуется на недостаточную информацию для построения ИП. Следующие цитаты из?predict.lm:

 The prediction intervals are for a single observation at each case
 in ‘newdata’ (or by default, the data used for the fit) with error
 variance(s) ‘pred.var’.  This can be a multiple of ‘res.var’, the
 estimated value of sigma^2: the default is to assume that future
 observations have the same error variance as those used for
 fitting.  If ‘weights’ is supplied, the inverse of this is used as
 a scale factor.  For a weighted fit, if the prediction is for the
 original data frame, ‘weights’ defaults to the weights used for
 the model fit, with a warning since it might not be the intended
 result.  If the fit was weighted and ‘newdata’ is given, the
 default is to assume constant prediction variance, with a warning.

Обратите внимание, что на конструкцию ДИ не влияет тип регрессии.

Делай все с нуля

В основном мы хотим знать, как получитьfit, se.fit, df а такжеresidual.scale вz.

Предсказанное среднее может быть вычислено умножением матрицы на векторXp %*% b, гдеXp является матрицей линейного предиктора иb вектор коэффициента регрессии.

Xp <- model.matrix(delete.response(terms(lmObject)), newdat)
b <- coef(lmObject)
yh <- c(Xp %*% b)  ## c() reshape the single-column matrix to a vector
#[1]  89.63133 104.66658

И мы видим, что это согласуется сz$fit, Дисперсия-ковариация дляyh являетсяXp %*% V %*% t(Xp), гдеV матрица дисперсии-ковариацииb который может быть вычислен

V <- vcov(lmObject)  ## use `vcov` function in R
#             (Intercept)         V2
# (Intercept)    7.862086 -1.1927966
# V2            -1.192797  0.2333733

Полная дисперсионно-ковариационная матрицаyh не требуется для вычисления точечного КИ или ПИ. Нам нужна только его главная диагональ. Так что вместо того, чтобы делатьdiag(Xp %*% V %*% t(Xp))мы можем сделать это более эффективно через

var.fit <- rowSums((Xp %*% V) * Xp)  ## point-wise variance for predicted mean
#       1        2 
#1.949963 2.598222 

sqrt(var.fit)  ## this agrees with `z$se.fit`
#       1        2 
#1.396411 1.611900 

Остаточная степень свободы легко доступна в оснащенной модели:

dof <- df.residual(lmObject)
#[1] 43

Наконец, чтобы вычислить остаточную дисперсию, используйте оценку Пирсона:

sig2 <- c(crossprod(lmObject$residuals)) / dof
# [1] 79.45063

sqrt(sig2)  ## this agrees with `z$residual.scale`
#[1] 8.913508

замечание

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

sig2 <- c(crossprod(sqrt(lmObject$weights) * lmObject$residuals)) / dof
Приложение: самописная функция, которая имитируетpredict.lm

Код в «Делать все с нуля» был четко организован в функциюlm_predict в этом Q & A:линейная модель сlm: как получить прогнозную дисперсию суммы прогнозируемых значений.

Ваш ответ на вопрос