Как функция 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))
2 ответа
При указании 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 : как получить прогнозную дисперсию суммы прогнозируемых значений.
Я не знаю, есть ли быстрый способ извлечь стандартную ошибку для интервала прогнозирования, но вы всегда можете отложить интервалы для 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,