Dlaczego przedział wiarygodności Bayesa w tej regresji wielomianowej jest tendencyjny, a przedział ufności jest prawidłowy?


9

Rozważ poniższy wykres, w którym symulowałem dane, w następujący sposób. Patrzymy na wynik binarny dla którego prawdziwe prawdopodobieństwo bycia 1 wskazuje czarna linia. Zależność funkcjonalna między współzmienną i jest wielomianem trzeciego rzędu z łączem logistycznym (więc jest nieliniowa w podwójnym kierunku).yobsxp(yobs=1|x)

Zielona linia to dopasowanie regresji logistycznej GLM, gdzie jest wprowadzany jako wielomian trzeciego rzędu. Przerywane zielone linie to 95% przedziały ufności wokół prognozy , gdzie dopasowane współczynniki regresji. Użyłem i do tego.xp(yobs=1|x,β^)β^R glmpredict.glm

Podobnie linia pruple jest średnią z odcinka tylnego z 95% wiarygodnym przedziałem dla bayesowskiego modelu regresji logistycznej z użyciem jednolitego wcześniejszego. Użyłem do tego pakietu z funkcją (ustawienie daje jednolity nieinformacyjny wcześniej).p(yobs=1|x,β)MCMCpackMCMClogitB0=0

Czerwone kropki oznaczają obserwacje w zbiorze danych, dla których , czarne kropki to obserwacje z . Zauważ, że jak często w klasyfikacji / analizie dyskretnej obserwuje się ale nie .yobs=1yobs=0yp(yobs=1|x)

wprowadź opis zdjęcia tutaj

Można zobaczyć kilka rzeczy:

  1. Celowo zasymulowałem, że jest rzadki na lewej ręce. Chcę, aby pewność siebie i wiarygodny przedział stały się tutaj szerokie z powodu braku informacji (obserwacji).x
  2. Oba przewidywania są tendencyjne w górę po lewej stronie. Ten błąd jest spowodowany czterema czerwonymi punktami oznaczającymi obserwacje, co błędnie sugeruje, że trafiłaby tutaj prawdziwa funkcjonalna forma. Algorytm nie ma wystarczających informacji, aby stwierdzić, że prawdziwa funkcjonalna forma jest wygięta w dół.yobs=1
  3. Przedział ufności zwiększa się zgodnie z oczekiwaniami, natomiast przedział wiarygodności nie . W rzeczywistości przedział ufności obejmuje całą przestrzeń parametrów, tak jak powinno to wynikać z braku informacji.

Wygląda na to, że wiarygodny przedział jest niewłaściwy / zbyt optymistyczny dla części . To naprawdę niepożądane zachowanie, aby wiarygodny przedział był wąski, gdy informacje stają się rzadkie lub całkowicie nieobecne. Zwykle nie tak reaguje wiarygodny interwał. Czy ktoś może wyjaśnić:x

  1. Jakie są tego powody?
  2. Jakie kroki mogę podjąć, aby uzyskać bardziej wiarygodny odstęp? (tzn. taki, który zawiera co najmniej prawdziwą formę funkcjonalną, lub lepiej staje się tak szeroki jak przedział ufności)

Kod umożliwiający uzyskanie przedziałów predykcji w grafice jest wydrukowany tutaj:

fit <- glm(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
x_pred <- seq(0, 1, by=0.01)
pred <- predict(fit, newdata = data.frame(x=x_pred), se.fit = T)
plot(plogis(pred$fit), type='l')
matlines(plogis(pred$fit + pred$se.fit %o% c(-1.96,1.96)), type='l', col='black', lty=2)


library(MCMCpack)
mcmcfit <- MCMClogit(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
gibbs_samps <- as.mcmc(mcmcfit)
x_pred_dm <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=x_pred))
gibbs_preds <- apply(gibbs_samps, 1, `%*%`, t(x_pred_dm))
gibbs_pis <- plogis(apply(gibbs_preds, 1, quantile, c(0.025, 0.975)))
matlines(t(gibbs_pis), col='red', lty=2)

Dostęp do danych : https://pastebin.com/1H2iXiew dzięki @DeltaIV i @AdamO


Jeśli ktoś mógłby mi wyjaśnić, jak udostępnić tabelę z danymi, mogę to zrobić.
tomka

Możesz użyć dputw ramce danych zawierającej dane, a następnie dołączyć dputdane wyjściowe jako kod w swoim poście.
DeltaIV

1
@ tomka oh widzę. Nie jestem ślepy na kolory, ale bardzo trudno mi dostrzec różnicę zielono-niebieską!
AdamO

1
@AdamO mam nadzieję, że tak jest lepiej
tomka

Odpowiedzi:


6

Dla częstościowym modelu wariancji powiększa predykcji proporcjonalnie do kwadratu odległości od środka ciężkości . Twoja metoda obliczania przedziałów predykcyjnych dla Bayesowskiego GLM wykorzystuje kwantyle empiryczne oparte na dopasowanej krzywej prawdopodobieństwa, ale nie uwzględnia dźwigniXX

Dwumianowy częsty GLM nie różni się od GLM z łączem tożsamości, z tym wyjątkiem, że wariancja jest proporcjonalna do średniej.

Zauważ, że dowolna wielomianowa reprezentacja prawdopodobieństw logitowych prowadzi do prognoz ryzyka, które są zbieżne z 0 jako i 1 jako lub odwrotnie, w zależności od znaku najwyższego terminu zamówienia wielomianowego .XX

W przypadku częstych predykcji dominuje ta tendencja do kwadratu odchylenia (dźwigni) proporcjonalnego wzrostu wariancji prognoz. To dlatego szybkość konwergencji do przedziałów prognozowania w przybliżeniu równych [0, 1] jest większa niż konwergencja logitów wielomianowych trzeciego rzędu do prawdopodobieństwa 0 lub 1 osobno.

Nie jest tak w przypadku Bayesowskich dopasowanych tylnych kwantyli. Nie ma wyraźnego zastosowania odchylenia kwadratowego, więc polegamy na proporcji dominujących tendencji 0 lub 1 do konstruowania długoterminowych przedziałów prognozowania.

Jest to oczywiste na podstawie ekstrapolacji bardzo daleko w skrajnych .X

Korzystając z kodu, który podałem powyżej, otrzymujemy:

> x_pred_dom <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=c(1000)))
> gibbs_preds <- plogis(apply(gibbs_samps[1000:10000, ], 1, `%*%`, t(x_pred_dom))) # a bunch of 0/1s basically past machine precision
> prop.table(table(gibbs_preds))
gibbs_preds
         0          1 
0.97733585 0.02266415 
> 

A więc 97,75% czasu, trzeci termin wielomianowy był ujemny. Można to zweryfikować na podstawie próbek Gibbs:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

Stąd przewidywane prawdopodobieństwo jest zbieżne z 0, gdy przechodzi w nieskończoność. Jeśli sprawdzimy SE modelu Bayesa, okaże się, że trzeci termin wielomianowy szacuje się na -185.25 z se 108,81, co oznacza, że ​​jest to 1,70 SD od 0, więc używając normalnych praw prawdopodobieństwa, powinien on spaść poniżej 0 95,5% czasu ( nie jest to zupełnie inna prognoza oparta na 10 000 iteracjach). Po prostu inny sposób zrozumienia tego zjawiska.X

Z drugiej strony, częste dopasowanie wieje do 0,1 zgodnie z oczekiwaniami:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

daje:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1

Niemniej jednak: czy model Bayesa nie jest zbyt pewny w obszarach danych , z których nie widział żadnych przykładów? Wiem, że bayesowskie posteriori lub rozkłady predykcyjne często zachowują się bardzo różnie (tj. Bardziej jak przedział konf.). Podejrzewam, że jest pewien wpływ przeora. Jeśli manipulować w określić precyzję normalne przed i można zaobserwować dość wpływ na wiarygodny przedział. xB0MCMClogit
tomka

@tomka Nie wiem, jak dokładnie na to odpowiedzieć, ponieważ wydaje się to styczne do pytania. Najważniejsze jest wskazanie, że te metody obliczania wskaźników PI nie są tak naprawdę porównywalne, zwłaszcza, że ​​dotyczą ekstrapolacji. Oczywiście, z wnioskowaniem bayesowskim, jeśli użyjesz pouczającego przeora, zyskasz skuteczność, gdy przeor ma rację, i tracisz, gdy przeor jest niewłaściwy.
AdamO

Chcę tylko powiedzieć, że wciąż myślę o twojej odpowiedzi. Nadal uważam, że to dziwne, że tylny nie reaguje na rzadkość poprzez poszerzenie. Uważam, że w przypadku innych priorów lepsze zachowanie w rzadkim regionie. Nie mogę tego dokładnie określić; Być może wzmocnię to pytanie przykładem, w którym wiarygodny przedział działa w sposób, jakiego bym się spodziewał, nawet w przypadku ekstrapolacji (szczególnie o normalnej liniowej regresji bayesowskiej). Kiedy to zrobię, dam ci znać.
tomka
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.