Jak obliczyć przedział ufności przecięcia X w regresji liniowej?

9

Ponieważ błąd standardowy regresji liniowej jest zwykle podawany dla zmiennej odpowiedzi, zastanawiam się, jak uzyskać przedziały ufności w innym kierunku - np. Dla przecięcia x. Jestem w stanie wyobrazić sobie, co to może być, ale jestem pewien, że musi istnieć prosty sposób, aby to zrobić. Poniżej znajduje się przykład w R na temat wizualizacji tego:

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)

wprowadź opis zdjęcia tutaj

Marc w pudełku
źródło
1
Można bootstrap to: 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)). W przypadku odwrotnych przedziałów prognoz plik pomocy chemCal:::inverse.predictzawiera następujące odniesienie, które może również pomóc w uzyskaniu CI: Massart, LM, Vandenginste, BGM, Buydens, LMC, De Jong, S., Lewi, PJ, Smeyers-Verbeke, J. (1997) ) Handbook of Chemometrics and Qualimetrics: Part A, str. 1. 200
Roland
1
To, co pokazano na wykresie, nie jest CI dla przechwytywania. Pokazujesz punkty, w których dolna i górna linia ufności prognoz przecina oś.
Roland
1
Często w regresji liniowej istnieje model, który mówi coś takiego:
Yja=α+βxja+εjagdzie ε1,εniid N.(0,σ2)),
tak aby Ys są traktowane jako losowe, a xs jako naprawione. Można to uzasadnić twierdzeniem, że szukasz rozkładu warunkowego, biorąc pod uwagęxs. W praktyce, jeśli weźmiesz nową próbkę, zwykle jest to nie tylkoYale także xzmiany te sugerują, że w niektórych okolicznościach należy je również uznać za losowe. Zastanawiam się, czy ma to związek z właściwością
Michael Hardy,
1
@AdrienRenaud - Wydaje mi się, że twoja odpowiedź jest zbyt uproszczona, biorąc pod uwagę asymetryczne aspekty, o których wspomniałem, i zostały podkreślone przez ćwiczenie ładowania, które ilustrował Roland. Jeśli nie pytam zbyt wiele, być może mógłbyś rozwinąć podejście oparte na prawdopodobieństwie, o którym wspomniałeś.
Marc w pudełku

Odpowiedzi:

9

Jak obliczyć przedział ufności przecięcia X w regresji liniowej?

Założenia

  • Użyj prostego modelu regresji yja=α+βxja+εja.
  • Błędy mają rozkład normalny zależny od regresorów ϵ|XN.(0,σ2)jan)
  • Dopasuj używając zwykłego najmniejszego kwadratu

3 procedury do obliczania przedziału ufności na przechwytywaniu x

Rozszerzenie Taylor pierwszego rzędu

Twój model to Y=zaX+b z szacowanym odchyleniem standardowym σza i σb na za i b parametry i szacowana kowariancja σzab. Ty rozwiązujesz

zaX+b=0X=-bza.

Następnie odchylenie standardowe σX na X jest dany przez:

(σXX)2)=(σbb)2)+(σzaza)2)-2)σzabzab.

MIB

Zobacz kod od Marc w polu „ Jak obliczyć przedział ufności punktu przecięcia x w regresji liniowej”? .

CAPITANI-POLLASTRI

CAPITANI-POLLASTRI zapewnia funkcję skumulowanego rozkładu i funkcję gęstości dla stosunku dwóch skorelowanych normalnych zmiennych losowych. Można go użyć do obliczenia przedziału ufności punktu przecięcia x w regresji liniowej. Ta procedura daje (prawie) identyczne wyniki jak te z MIB.

Rzeczywiście, używając zwykłego najmniejszego kwadratu i zakładając normalność błędów, β^N.(β,σ2)(XT.X)-1) (zweryfikowane) i β^są skorelowane (zweryfikowane).

Procedura jest następująca:

  • pobierz estymator OLS dla za i b.
  • uzyskaj macierz wariancji-kowariancji i wyodrębnij, σza,σb,σzab=ρσzaσb.
  • Zakładać, że za i b postępuj zgodnie z dwuwymiarowym skorelowanym rozkładem normalnym, N.(za,b,σza,σb,ρ). Następnie funkcja gęstości i funkcja rozkładu skumulowanego zxjantmirdomipt=-bza są podawane przez CAPITANI-POLLASTRI.
  • Użyj funkcji kumulatywnej dystrybucji z xjantmirdomipt=-bza obliczyć pożądane kwantyle i ustawić przedział ufności.

Porównanie 3 procedur

Procedury są porównywane przy użyciu następującej konfiguracji danych:

  • x <- 1:10
  • a <- 20
  • b <- -2
  • y <- a + b * x + rnorm (długość (x), średnia = 0, sd = 1)

10000 różnych próbek jest generowanych i analizowanych przy użyciu 3 metod. Kod (R) używany do generowania i analizy można znaleźć na stronie : https://github.com/adrienrenaud/stackExchange/blob/master/crossValidated/q221630/answer.ipynb

  • MIB i CAPITANI-POLLASTRI dają równoważne wyniki.
  • Ekspansja Taylora pierwszego rzędu różni się znacznie od dwóch pozostałych metod.
  • MIB i CAPITANI-POLLASTRI cierpią z powodu niedostatecznego zasięgu. Stwierdzono, że 68% (95%) ci zawiera prawdziwą wartość 63% (92%) czasu.
  • Rozszerzenie Taylor pierwszego rzędu cierpi z powodu nadmiernego zasięgu. Stwierdzono, że 68% (95%) ci zawiera prawdziwą wartość 87% (99%) czasu.

Wnioski

Rozkład punktów przecięcia x jest asymetryczny. Uzasadnia to asymetryczny przedział ufności. MIB i CAPITANI-POLLASTRI dają równoważne wyniki. CAPITANI-POLLASTRI mają fajne uzasadnienie teoretyczne i daje podstawy dla MIB. MIB i CAPITANI-POLLASTRI cierpią z powodu umiarkowanego niedostatecznego zasięgu i mogą być używane do ustalania przedziałów ufności.

Adrien Renaud
źródło
Dzięki za tę miłą odpowiedź. Czy ta metoda sugeruje, że standardowy błąd przechwytu x jest symetryczny? Interwały przewidywania na mojej figurze sugerują, że tak nie jest i widziałem odniesienie do tego gdzie indziej.
Marc w pudełku
Tak, oznacza to symetryczny interwał. Jeśli chcesz mieć asymetryczny, możesz użyć prawdopodobieństwa profilu, traktując parametry modelu jako parametry uciążliwe. Ale to więcej pracy :)
Adrien Renaud
Czy możesz wyjaśnić bardziej szczegółowo, w jaki sposób otrzymujesz to wyrażenie (σX/X)2)?
@fcop To rozszerzenie Taylora. Spójrz na en.wikipedia.org/wiki/Propagation_of_uncertainty
Adrien Renaud
2

Polecam ładowanie resztek:

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")

wynikowa fabuła

Na wykresie pokazano punkty, w których dolna / górna granica przedziału ufności prognoz przecina oś. Nie sądzę, że są to granice ufności przechwytywania, ale może są to przybliżone przybliżenie.

Roland
źródło
Świetnie - to już wygląda rozsądniej niż w przykładzie z twojego komentarza. Dzięki jeszcze raz.
Marc w pudełku