rtanarm für die Bayessche hierarchische Modellierung von Binomialexperimenten

Question

Angenommen, es werden drei Binomialexperimente chronologisch durchgeführt. Für jedes Experiment kenne ich die # of trials sowie die #von Erfolgen. Um die ersten beiden älteren Experimente als vorherige für das dritte Experiment zu verwenden, möchte ich "ein Bayes'sches hierarchisches Modell an die beiden älteren Experimente anpassen und die posteriore Form verwenden, die als Prior für das dritte Experiment verwendet wurde".

Angesichts meiner verfügbaren Daten (unten) ist meine Frage: ist mein rstanarm Code unter Erfassung, was ich oben beschrieben?

Study1_trial = 70 
Study1_succs = 27 
#================== 
Study2_trial = 84 
Study2_succs = 31 
#================== 
Study3_trial = 100 
Study3_succs = 55 

Was ich in Paket versucht habe rstanarm:

library("rstanarm") 

data <- data.frame(n = c(70, 84, 100), y = c(27, 31, 55)); 
mod <- stan_glm(cbind(y, n - y) ~ 1, prior = NULL, data = data, family = binomial(link = 'logit')) 

## can I use a beta(1.2, 1.2) as prior for the first experiment? 

Answer 1

TL; DR: Wenn Sie direkt die Wahrscheinlichkeit des Erfolgs wurden die Vorhersage, würde das Modell eine Bernoulli-Wahrscheinlichkeit mit dem Parameter Theta sein (die Erfolgswahrscheinlichkeit), die Werte zwischen null und eins annehmen können. Sie könnten in diesem Fall eine Beta für Theta verwenden. Aber mit einem logistischen Regressionsmodell modellieren Sie tatsächlich die Log-Erfolgsquoten, die einen beliebigen Wert von "-Inf" bis "Inf" annehmen können, also ein Prior mit einer normalen Verteilung (oder ein anderer Prior, der einen echten Wert annehmen kann) ein Bereich, der durch die verfügbaren vorherigen Informationen bestimmt wird, wäre angemessener.

---

Für ein Modell, bei dem der einzige Parameter, der Schnittpunkt ist, ist die vor der Wahrscheinlichkeitsverteilung für die Protokolle Chancen auf Erfolg. Mathematisch ist das Modell:

log(p/(1-p)) =  a 

Wo p die Wahrscheinlichkeit des Erfolgs ist und a, der Parameter Sie schätzen, ist der Abschnitt, der jede reelle Zahl sein kann. Wenn die Erfolgsaussichten 1: 1 (dh p = 0,5) sind, dann a = 0. Wenn die Chancen größer als 1: 1 sind, ist a positiv. Wenn die Chancen weniger als 1: 1 sind, ist a negativ.

Da wir einen Prior für a wollen, brauchen wir eine Wahrscheinlichkeitsverteilung, die jeden realen Wert annehmen kann. Wenn wir nichts über die Erfolgsaussichten wissen würden, könnten wir einen sehr schwach informativen Prior wie eine Normalverteilung mit, sagen wir, Mittelwert = 0 und sd = 10 verwenden (das ist der rstanarm Standard), was bedeutet, dass eine Standardabweichung würde umfassen Erfolgsquoten von ca. 22000: 1 bis 1: 22000! Dieser Prior ist also im Wesentlichen flach.

Wenn wir Ihre ersten beiden Studien nehmen die vor zu konstruieren, können wir die Wahrscheinlichkeitsdichte verwenden diesen Studien basierend auf und es dann in die Log-Odds-Transformation skalieren:

# Possible outcomes (that is, the possible number of successes) 
s = 0:(70+84) 

# Probability density over all possible outcomes 
dens = dbinom(s, 70+84, (27+31)/(70+84)) 

Angenommen, wir ein normales verwenden werden Bei der Verteilung für den Prior wollen wir die wahrscheinlichste Erfolgswahrscheinlichkeit (die für den Prior der Mittelwert ist) und die Standardabweichung des Mittelwerts.

# Prior parameters 
pp = s[which.max(dens)]/(70+84) # most likely probability 
psd = sum(dens * (s/max(s) - pp)^2)^0.5 # standard deviation 

# Convert prior to log odds scale 
pp_logodds = log(pp/(1-pp)) 
psd_logodds = log(pp/(1-pp)) - log((pp-psd)/(1 - (pp-psd))) 

c(pp_logodds, psd_logodds) 

[1] -0.5039052 0.1702006 

Sie könnten erzeugen im Wesentlichen den gleichen Stand von stan_glm auf den ersten beiden Studien mit dem Standardlauf (flach) vor:

prior = stan_glm(cbind(y, n-y) ~ 1, 
       data = data[1:2,], 
       family = binomial(link = 'logit')) 

c(coef(prior), se(prior)) 

[1] -0.5090579 0.1664091 

Jetzt lass uns f Es verwendet das Modell aus den Daten von Studie 3 unter Verwendung der Standardvorgabe und der Vorabversion, die wir gerade erstellt haben. Ich habe zu einem Standard-Datenrahmen gewechselt, da stan_glm scheinbar fehlschlägt, wenn der Datenrahmen nur eine Zeile hat (wie in data = data[3, ]).

Zum Vergleich, lassen Sie uns auch ein Modell mit allen drei Studien und der Standard-Wohnung vorher generieren. Wir würden dieses Modell praktisch die gleichen Ergebnisse wie mod2 zu geben, erwarten:

mod3 <- stan_glm(cbind(y, n - y) ~ 1, 
       data = data, 
       family = binomial(link = 'logit')) 

Jetzt ist die drei Modelle lassen vergleichen:

library(tidyverse) 

list(`Study 3, Flat Prior`=mod1, 
    `Study 3, Prior from Studies 1 & 2`=mod2, 
    `All Studies, Flat Prior`=mod3) %>% 
    map_df(~data.frame(log_odds=coef(.x), 
        p_success=predict(.x, type="response")[1]), 
     .id="Model") 

       Model log_odds p_success 
1    Study 3, Flat Prior 0.2008133 0.5500353 
2 Study 3, Prior from Studies 1 & 2 -0.2115362 0.4473123 
3   All Studies, Flat Prior -0.2206890 0.4450506 

Für Studie 3 mit der flachen vor (Zeile 1), ist die vorhergesagte Erfolgswahrscheinlichkeit wie erwartet 0,55, da dies die Daten aussagen und der Prior keine zusätzliche Information liefert.

Für Studie 3 mit einem früheren auf Studien 1 und 2 basiert, ist die Wahrscheinlichkeit des Erfolgs 0,45. Die geringere Erfolgswahrscheinlichkeit beruht auf der geringeren Erfolgswahrscheinlichkeit in den Studien 1 und 2, die zusätzliche Informationen hinzufügen. In der Tat ist die Wahrscheinlichkeit des Erfolgs von mod2 genau das, was Sie direkt aus den Daten berechnen würden: with(data, sum(y)/sum(n)). mod3 setzt alle Informationen in die Wahrscheinlichkeit ein, anstatt sie zwischen dem Vorherigen und der Wahrscheinlichkeit zu teilen, ist aber ansonsten im Wesentlichen dasselbe wie mod2.