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

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

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

и моя задача состояла в том, чтобы получить

  • 90% доверительный интервал для среднего отклика при V2=6 и
  • 90% -ный интервал outlookирования, когда 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) что кажется правильным, поскольку PI должен быть шире.

Выход для обоих также включал 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)) 

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

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

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

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


Настроить

Я собираю ваш код здесь, чтобы помочь другим читателям копировать, вставлять и запускать. Я также изменяю имена переменных, чтобы они имели более четкие значения. Кроме того, я расширяю 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 - стандартная ошибка outlookируемого среднего z$fit , используемая для построения CI для 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 шире, чем 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 Вы не устанавливаете newdata при использовании type = "prediction" в predict.lm ), потому что веса известны (вы должны были предоставить его с помощью аргумента weight при использовании lm ) , Для outlookирования вне выборки (т. newdata Вы передаете newdata для predict.lm ), predict.lm ожидает, что вы скажете ему, как следует учитывать взвешенную дисперсию. Вам нужно либо использовать аргумент pred.var либо pred.var в predict.lm , иначе вы получите предупреждение от predict.lm на недостаточную информацию для построения PI. Ниже приводятся следующие данные из. ?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. 

Обратите внимание, что тип CI не зависит от типа регрессии.


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

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

Прогнозируемое среднее значение может быть вычислено умножением матриц-векторов Xp %*% b , где Xp - matrix линейного предсказателя, а 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 - matrix дисперсии-ковариации b которую можно вычислить по формуле

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

Полная matrix дисперсии-ковариации yh не требуется для вычисления точечного CI или PI. Нам нужна только его основная диагональ. Поэтому вместо выполнения 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 функция, имитирующая predict.lm

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

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

Interesting Posts

MATLAB: применить фильтр низких частот или верхних частот к массиву

std :: пара ссылок

Как получить временную метку ближайшего ключевого кадра перед заданной меткой времени с помощью FFmpeg?

Spring AOP против AspectJ

Как предоставить файлы данных для тестов на андроид

В чем разница между «1L» и «1»?

Как получить текстовый узел элемента?

LibreOffice Calc: Как получить общее количество для HH: MM: SS-ячейки

Масштабировать пользовательский интерфейс для нескольких разрешений / различных устройств

Компрессор Flash и изображений

Как это сделать, я получаю вывод команды оболочки, выполненной с использованием переменной из Jenkinsfile (groovy)?

union ‘punning’ structs w / “common initial sequence”: Почему C (99+), но не C ++, предусматривает «видимое объявление типа объединения»?

Bash 4 ассоциативных массива: ошибка “объявить: -A: недействительный вариант”

Как получить доступные точки доступа Wi-Fi и их уровень сигнала в .net?

Найти общие подстроки между двумя символьными переменными

Давайте будем гением компьютера.