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

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).$y_{obs}$$x$$p(y_{obs}=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.$x$$p(y_{obs}=1 | x, \hat{\beta})$$\hat{\beta}$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(y_{obs}=1 | x, \beta)$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 .$y_{obs}=1$$y_{obs}=0$$y$$p(y_{obs}=1 | x)$

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ół.$y_{obs}=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

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źwigni$X$$X$

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 .$X\rightarrow -\infty$$X\rightarrow \infty$

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