Как рассчитать доверительный интервал X-перехвата в линейной регрессии?

9

Поскольку стандартная ошибка линейной регрессии обычно задается для переменной отклика, мне интересно, как получить доверительные интервалы в другом направлении - например, для x-перехвата. Я могу представить, что это может быть, но я уверен, что должен быть прямой способ сделать это. Ниже приведен пример того, как R визуализировать это:

set.seed(1)
x <- 1:10
a <- 20
b <- -2
y <- a + b*x + rnorm(length(x), mean=0, sd=1)

fit <- lm(y ~ x)
XINT <- -coef(fit)[1]/coef(fit)[2]

plot(y ~ x, xlim=c(0, XINT*1.1), ylim=c(-2,max(y)))
abline(h=0, lty=2, col=8); abline(fit, col=2)
points(XINT, 0, col=4, pch=4)
newdat <- data.frame(x=seq(-2,12,len=1000))

# CI
pred <- predict(fit, newdata=newdat, se.fit = TRUE) 
newdat$yplus <-pred$fit + 1.96*pred$se.fit 
newdat$yminus <-pred$fit - 1.96*pred$se.fit 
lines(yplus ~ x, newdat, col=2, lty=2)
lines(yminus ~ x, newdat, col=2, lty=2)

# approximate CI of XINT
lwr <- newdat$x[which.min((newdat$yminus-0)^2)]
upr <- newdat$x[which.min((newdat$yplus-0)^2)]
abline(v=c(lwr, upr), lty=3, col=4)

введите описание изображения здесь

Марк в коробке
источник
1
Вы можете самонастройки это: library(boot); sims <- boot(data.frame(x, y), function(d, i) { fit <- lm(y ~ x, data = d[i,]) -coef(fit)[1]/coef(fit)[2] }, R = 1e4); points(quantile(sims$t, c(0.025, 0.975)), c(0, 0)). Для интервалов обратного прогнозирования в файле справки chemCal:::inverse.predictприводится следующая справка, которая также может помочь при получении КИ: Massart, LM, Vandenginste, BGM, Buydens, LMC, Де Йонг, С., Леви, PJ, Смайерс-Вербеке, J. (1997 ) Справочник по хемометрике и квалиметрии: часть А, с. 200
Роланд
1
То, что вы показываете на графике, не является CI для перехвата. Вы показываете точки, где нижняя и верхняя доверительные линии прогнозов пересекают ось.
Роланд
1
Yi=α+βxi+εiwhere ε1,εni.i.d. N(0,σ2),
YxxYx
Майкл Харди
1
@AdrienRenaud - Мне кажется, что ваш ответ слишком упрощен, учитывая асимметричные аспекты, которые я упомянул, и выделен в процессе начальной загрузки, который проиллюстрировал Роланд. Если я не слишком много спрашиваю, возможно, вы могли бы расширить подход, о котором вы упомянули.
Марк в коробке

Ответы:

9

Как рассчитать доверительный интервал X-перехвата в линейной регрессии?

Asumptions

  • yi=α+βxi+εi
  • ϵ|XN(0,σ2In)
  • Подходит, используя обычный метод наименьших квадратов

3 процедуры для расчета доверительного интервала на x-intercept

Расширение Тейлора первого порядка

Y=aX+bσaσbabσab

aX+b=0X=ba.

σXX

(σXX)2=(σbb)2+(σaa)22σabab.

MIB

См. Код от Марка во вставке в разделе Как рассчитать доверительный интервал x-пересечения в линейной регрессии? ,

CAPITANI-POLLASTRI

CAPITANI-POLLASTRI предоставляет кумулятивную функцию распределения и функцию плотности для отношения двух коррелированных нормальных случайных величин. Он может использоваться для вычисления доверительного интервала x-пересечения в линейной регрессии. Эта процедура дает (почти) те же результаты, что и результаты MIB.

β^N(β,σ2(XTX)1)β^

Процедура следующая:

  • ab
  • σa,σb,σab=ρσaσb
  • abN(a,b,σa,σb,ρ)xintercept=ba
  • xintercept=ba

Сравнение 3 процедур

Процедуры сравниваются с использованием следующей конфигурации данных:

  • х <- 1:10
  • <- 20
  • б <- -2
  • y <- a + b * x + rnorm (длина (x), среднее значение = 0, sd = 1)

10000 различных образцов генерируются и анализируются с использованием 3 методов. Код (R), использованный для генерации и анализа, можно найти по адресу: https://github.com/adrienrenaud/stackExchange/blob/master/crossValidated/q221630/answer.ipynb.

  • MIB и CAPITANI-POLLASTRI дают эквивалентные результаты.
  • Разложение Тейлора первого порядка существенно отличается от двух других методов.
  • MIB и CAPITANI-POLLASTRI страдают от недостаточного покрытия. Обнаружено, что 68% (95%) ci содержат истинное значение 63% (92%) времени.
  • Расширение Тейлора первого порядка страдает от чрезмерного покрытия. Обнаружено, что 68% (95%) ci содержат истинное значение 87% (99%) времени.

Выводы

Распределение x-перехвата асимметрично. Это оправдывает асимметричный доверительный интервал. MIB и CAPITANI-POLLASTRI дают эквивалентные результаты. У CAPITANI-POLLASTRI есть хорошее теоретическое обоснование, и оно дает основания для MIB. MIB и CAPITANI-POLLASTRI страдают от умеренного недостаточного покрытия и могут использоваться для установки доверительных интервалов.

Адриен Рено
источник
Спасибо за этот хороший ответ. Означает ли этот метод, что стандартная ошибка x-пересечения является симметричной? Интервалы предсказания на моей фигуре подразумевают, что это не так, и я видел ссылку на это в другом месте.
Марк в коробке
Да, это подразумевает симметричный интервал. Если вы хотите асимметричный, вы можете использовать вероятность профиля, рассматривая параметры вашей модели как параметры неприятности. Но это больше работы :)
Адриен Рено
(σX/X)2
@fcop Это расширение Тейлора. Взгляните на en.wikipedia.org/wiki/Propagation_of_unterminty
Адриен Рено
2

Я бы порекомендовал начальную загрузку остатков:

library(boot)

set.seed(42)
sims <- boot(residuals(fit), function(r, i, d = data.frame(x, y), yhat = fitted(fit)) {

  d$y <- yhat + r[i]

  fitb <- lm(y ~ x, data = d)

  -coef(fitb)[1]/coef(fitb)[2]
}, R = 1e4)
lines(quantile(sims$t, c(0.025, 0.975)), c(0, 0), col = "blue")

результирующий сюжет

На графике показаны точки, в которых нижний / верхний предел доверительной полосы прогнозов пересекает ось. Я не думаю, что это пределы достоверности перехвата, но, возможно, это грубое приближение.

Roland
источник
Отлично - это уже выглядит более разумно, чем пример из вашего комментария. Еще раз спасибо.
Марк в коробке