이항 실험의 베이지안 계층 적 모델링을위한 rstanarm

Question

시간 경과에 따라 수행 된 3 가지 이항 실험이 있다고 가정합니다. 각 실험에 대해 번과뿐 아니라 # 번을 알고 있습니다. 세 번째 실험에서 처음 두 이전 실험을 사용하려면 "두 이전 실험에 베이지안 계층 모델을 적용하고 세 번째 실험의 경우 이전에 후방 양식을 사용하십시오".

아래에서 제 데이터를 볼 때 제 질문은 : 위의 내용을 캡처하는 아래 코드는 rstanarm입니까?

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

내가 패키지 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 : 직접 성공의 확률을 예측한다면은 모델 (매개 변수 세타와 베르누이 가능성이 될 것입니다 성공 확률)로 0과 1 사이의 값을 취할 수 있습니다. 이 경우에는 theta에 대한 베타 사전을 사용할 수 있습니다. 그러나 로지스틱 회귀 모델을 사용하면 성공의 로그 확률을 실제로 모델링 할 수 있습니다.이 확률은 -Inf에서 Inf로 어떤 값을 취할 수 있으므로 정규 분포를 사용하여 이전 값을 구할 수 있습니다. 이용 가능한 선행 정보에 의해 결정된 일부 범위)이 더 적절할 것이다.

---

유일한 매개 변수가 절편 인 모델의 경우, 우선 순위는 로그 확률 성공 확률 분포입니다. 수학적 모델은 다음과 같습니다

p 성공과 a의 확률

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

, 당신이 추정하고있는 매개 변수는 실제 숫자가 될 수 있습니다 절편이다. 성공 확률이 1 : 1 (즉, p = 0.5)이면 a = 0입니다. 확률이 1 : 1보다 큰 경우 a은 양수입니다. 확률이 1 : 1보다 작 으면 a이 음수입니다.

a에 대한 이전 값을 원하기 때문에 실제 값을 취할 수있는 확률 분포가 필요합니다. 우리가 성공 확률에 대해 몰랐다면 mean = 0과 sd = 10 (보통은 rstanarm 기본값)과 같이 정규 분포와 같이 매우 약한 유익한 사전을 사용할 수 있습니다. 이는 하나의 표준 편차가 성공의 확률은 약 22,000 : 1에서 1 : 22000까지입니다! 따라서이 사전은 본질적으로 평평합니다. 우리가 이전을 구성하는 첫 번째 두 연구를 취할 경우

, 우리가 그 연구를 바탕으로 확률 밀도를 사용하고 로그 확률로 변환 할 수는 확장 :

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

우리는 정상을 사용합니다 가정 사전에 대한 분포, 우리는 성공 확률 (이전에 대한 평균)과 평균의 표준 편차를 원합니다.

# 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 

같은 이전 기본값 (평면) 전에 처음 두 연구에 stan_glm을 실행하여 본질적으로 생성 할 수있는 :

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 

이제하자. 그것은 우리가 방금 생성 한 사전과 사전을 사용하여 스터디 3의 데이터를 사용하는 모델입니다. 나는 stan_glm이 데이터 프레임에 단 하나의 행을 가지고있을 때 실패한 것으로 보았 기 때문에 (data = data[3, ]에서와 같이) 표준 데이터 프레임으로 전환했다.

# Default weakly informative prior 
mod1 <- stan_glm(y ~ 1, 
       data = data.frame(y=rep(0:1, c(45,55))), 
       family = binomial(link = 'logit')) 

# Prior based on studies 1 & 2 
mod2 <- stan_glm(y ~ 1, 
       data = data.frame(y=rep(0:1, c(45,55))), 
       prior_intercept = normal(location=pp_logodds, scale=psd_logodds), 
       family = binomial(link = 'logit')) 

비교를 들어, 모든 세 가지 연구와 모델과 이전 평면의 기본을 생성 할 수 있습니다. (평면 이전에 연구 3

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 

: 이제

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

의이 세 가지 모델을 비교하자 : 우리는이 모델이 거의 mod2과 같은 결과를 제공 기대 행 1) 예측 된 성공 확률은 예상대로 0.55이며 이는 데이터가 말하고 이전은 추가 정보를 제공하지 않기 때문입니다.

연구 1과 2를 바탕으로 한 연구 3의 경우, 성공 확률은 0.45입니다. 성공 확률이 낮 으면 연구 1과 2에서 성공 가능성이 낮기 때문에 추가 정보가 추가되기 때문입니다. 사실 mod2의 성공 확률은 데이터에서 직접 계산 한 것입니다 : with(data, sum(y)/sum(n)). mod3은 모든 정보를 이전 및 가능성 사이에서 분리하는 대신 가능성에 넣지만, 그렇지 않으면 본질적으로 mod2과 동일합니다.