rstanarm для байесовского иерархического моделирования биномиальных экспериментов
Предположим, есть три биномиальных эксперимента, проведенных в хронологическом порядке. Для каждого эксперимента я знаю количество испытаний и успехов. Чтобы использовать первые два более старых эксперимента, как и в предыдущем третьем эксперименте, я хочу "подогнать байесовскую иерархическую модель к двум более старым экспериментам и использовать заднюю форму, которая предшествует третьему эксперименту".
Учитывая мои доступные данные (ниже), мой вопрос: мой 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?
1 ответ
TL;DR: если бы вы непосредственно предсказывали вероятность успеха, модель имела бы вероятность Бернулли с параметром тета (вероятность успеха), которая могла бы принимать значения от нуля до единицы. В этом случае вы можете использовать бета-версию для тета. Но с помощью модели логистической регрессии вы на самом деле моделируете шансы на успех, которые могут принимать любое значение от -Inf до Inf, таким образом, априор с нормальным распределением (или какой-то другой априор, который может принимать любое реальное значение в пределах некоторый диапазон, определенный доступной предшествующей информацией), был бы более соответствующим.
Для модели, где единственным параметром является перехват, предшествующим является распределение вероятностей для логарифмических шансов на успех. Математически модель представляет собой:
log(p/(1-p)) = a
куда p вероятность успеха и aоцениваемый параметр - это перехват, который может быть любым действительным числом. Если шансы на успех составляют 1:1 (то есть р = 0,5), то a = 0, Если шансы больше 1: 1, то a положительно. Если шансы меньше 1: 1, то a отрицательно.
Так как мы хотим, чтобы до aНам нужно распределение вероятностей, которое может принимать любое реальное значение. Если бы мы ничего не знали о шансах на успех, мы могли бы использовать очень слабо информативный априор, такой как нормальное распределение с, скажем, mean=0 и sd=10 (это rstanarm по умолчанию), что означает, что одно стандартное отклонение будет включать в себя шансы на успех в диапазоне от 22000: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'))
Для сравнения давайте также сгенерируем модель со всеми тремя исследованиями и по умолчанию с фиксированным значением. Мы ожидаем, что эта модель даст практически те же результаты, что и mod2:
mod3 <- stan_glm(cbind(y, n - y) ~ 1,
data = data,
family = binomial(link = 'logit'))
Теперь давайте сравним три модели:
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
Для исследования 3 с фиксированной априорной (строка 1) прогнозируемая вероятность успеха составляет 0,55, как и ожидалось, поскольку именно об этом говорят данные, а предыдущая не дает никакой дополнительной информации.
Для исследования 3 с предварительным, основанным на исследованиях 1 и 2, вероятность успеха составляет 0,45. Меньшая вероятность успеха обусловлена меньшей вероятностью успеха в исследованиях 1 и 2 с добавлением дополнительной информации. На самом деле вероятность успеха от mod2 это именно то, что вы рассчитываете непосредственно из данных: with(data, sum(y)/sum(n)), mod3 помещает всю информацию в правдоподобие вместо того, чтобы разделять ее между предыдущим и правдоподобием, но в остальном по сути такой же, как mod2,