Interwał prognoz dla przyszłej proporcji sukcesów w ustawieniach dwumianowych

Załóżmy, że dopasowuję regresję dwumianową i uzyskuję oszacowania punktowe i macierz wariancji-kowariancji współczynników regresji. To pozwoli mi uzyskać CI dla oczekiwanego odsetka sukcesów w przyszłym eksperymencie, $p$ , ale potrzebuję CI dla obserwowanej proporcji. Opublikowano kilka powiązanych odpowiedzi, w tym symulację (załóżmy, że nie chcę tego robić) oraz link do Krishnamoorthya i in. (Który nie do końca odpowiada na moje pytanie).

Moje rozumowanie jest następujące: jeśli używamy tylko modelu dwumianowego, jesteśmy zmuszeni to założyć $p$ jest próbkowany z rozkładu normalnego (z odpowiednim Wald CI) i dlatego niemożliwe jest uzyskanie CI dla obserwowanej proporcji w postaci zamkniętej. Jeśli to założymy $p$ jest próbkowany z rozkładu beta, wtedy rzeczy są znacznie łatwiejsze, ponieważ liczba sukcesów będzie podążać za rozkładem dwumianowym. Będziemy musieli założyć, że nie ma niepewności w szacowanych parametrach beta, $\alpha$ i $\beta$ .

Istnieją trzy pytania:

1) Teoretyczny: czy można stosować tylko oszacowania punktowe parametrów beta? Wiem, że aby zbudować CI do przyszłej obserwacji w wielokrotnej regresji liniowej

$Y = x'\beta + \epsilon, \epsilon \sim N(0, \sigma^2)$

robią tę niepoprawność wariancji terminu, $\sigma^2$ . Rozumiem (popraw mnie, jeśli się mylę), że w praktyce jest to uzasadnienie $\sigma^2$ jest szacowany z dużo większą precyzją niż współczynniki regresji i nie zyskamy wiele, próbując uwzględnić niepewność $\sigma^2$ . Czy podobne uzasadnienie ma zastosowanie do szacowanych parametrów beta, $\alpha$ i $\beta$ ?

2) Jaki pakiet jest lepszy (R: gamlss-bb, betareg, aod ?; Mam również dostęp do SAS).

3) Biorąc pod uwagę oszacowane parametry beta, czy istnieje (przybliżony) skrót, aby uzyskać kwantyle (2,5%, 97,5%) w celu zliczenia przyszłych sukcesów lub, jeszcze lepiej, w odniesieniu do odsetka przyszłych sukcesów w rozkładzie dwumianowym.

— James
źródło

Na pytanie pierwsze, tak, jest to ważna rzecz, którą ludzie robią, nazywa się to Empirical Bayes: en.wikipedia.org/wiki/Empirical_Bayes_method

— Paul

Nie sądzę, że użycie metody XYZ do oszacowania parametru modelu może automatycznie sugerować, że można zignorować niepewność oszacowania podczas tworzenia CI dla przyszłych obserwacji. Np. W wielokrotnej regresji liniowej używają OLS zamiast EB, a niepewność w

σ

$\sigma$ jest również ignorowane. Dlaczego? Ponadto, ten artykuł Wiki nigdy nie sugeruje, że w EB precyzja szacowania hiperparametrów najwyższego poziomu jest zwykle o wiele wyższa, że można uznać je za ustalone ze względów praktycznych.

— James

„Kiedy prawdziwa dystrybucja

p (η ∣ y)

$p(\eta \mid y)$ jest ostro zakończony, co stanowi całkę

p (θ ∣ y)

$p(\theta \mid y)$ może być niewiele zmieniony poprzez zastąpienie rozkładu prawdopodobieństwa

η

$\eta$ z oszacowaniem punktowym

η^{*}

$\eta^*$ reprezentujący szczyt rozkładu ”. To, czy tak jest w twoim przypadku, zależy od specyfiki Twojej problematycznej domeny.

— Paul

Dobre pytanie! Nie można uzyskać osi przestawnej, ale co z wykorzystaniem prawdopodobieństwa profilu? Zobacz Jakie nie-bayesowskie metody wnioskowania predykcyjnego? .

— Scortchi - Przywróć Monikę

Odpowiem na wszystkie 3 części.

Istnieją dwa skonfliktowane problemy, po pierwsze metoda stosowana w tym przypadku w celu dopasowania modelu regresji. Po drugie, jak podzielić szacunki na podstawie szacunków, aby przewidzieć nowe oszacowanie.

jeśli twoje zmienne odpowiedzi są rozkładane dwumianowo, zwykle używałbyś regresji logistycznej lub regresji probit (glm z normalnym cdf jako funkcją link).

W przypadku regresji logistycznej odpowiedzią będzie stosunek obserwowanych zliczeń podzielony przez znaną górną granicę, tj. $y_i/n_i$ . Następnie weź swoje predyktory / zmienne towarzyszące i umieść je w wywołaniu R funkcji glm. Zwrócony obiekt zawiera wszystko, czego potrzebujesz do wykonania pozostałych obliczeń.

x<- rnorm(100, sd=2)
prob_true <- 1/(1+exp(-(1+5*x)))
counts <- rbinom(100, 50,prob_true)
print(d.AD <- data.frame(counts,x))
glm.D93 <- glm(counts/50 ~ x, family = binomial() )

W przypadku modelu regresji liniowej wzór na przedział predykcji jest następujący:

$\hat{y}_i \pm t_{n-p}s_y\sqrt{1+\frac{1}{n}+\frac{(x_i-\bar{x})^2}{(n-1)s^2_x}}$

Możesz użyć modelu regresji liniowej jako przybliżenia glm. Aby to zrobić, należy wykonać formułę regresji liniowej dla liniowej kombinacji predyktorów przed wykonaniem odwrotnej transformacji łącza w celu przywrócenia prawdopodobieństwa z powrotem w skali 0-1. Kod służący do tego jest zapisywany w funkcji przewidywanej R.glm () R. Oto przykładowy kod, który również stworzy niezłą fabułę. ( EDYCJA : Ten kod dotyczy przedziału ufności, a nie przedziału prognozowania)

y_hat <- predict(glm.D93, type="link", se.fit=TRUE)
t_np<- qt(.975, 100-2, ncp=0)

ub <- y_hat$fit + t_np * y_hat$se.fit
lb <- y_hat$fit - t_np * y_hat$se.fit

point <- y_hat$fit

p_hat <- glm.D93$family$linkinv(point)
p_hat_lb <- glm.D93$family$linkinv(lb)
p_hat_ub <- glm.D93$family$linkinv(ub)

plot(x,p_hat)
points(x, p_hat_ub, col='red')
points(x, p_hat_lb, col='blue')

Możesz zrobić to samo dla dowolnego glm, np. Poissona, odwrotnego gaussa, gamma itp. W każdym przypadku wykonaj przedział predykcji na skali liniowej kombinacji predyktorów. Po uzyskaniu dwóch punktów końcowych przedziału prognozy przekształcasz te punkty końcowe za pomocą łącza odwrotnego. Dla każdego z wymienionych przeze mnie glms link odwrotny może być inny niż logit, który tu napisałem. Mam nadzieję że to pomoże.

— Lucas Roberts
źródło