Подход уникальных перехватов для категориальных переменных в пакете "rstanarm" в R

Question

Справочная информация: McElearth (2016) в своих переосмысленных страницах книги 158-159 использует переменную индекса вместо фиктивного кодирования для переменной 3 категории, называемой "clade", для предсказания "kcal.per.g" (линейная регрессия).

Вопрос: мне было интересно, могли бы мы применить тот же подход в "rstanarm"? Я предоставил данные и код R для возможной демонстрации ниже.

library("rethinking") # A github package not on CRAN
data(milk)
d <- milk
d$clade_id <- coerce_index(d$clade) # Index variable maker
#[1] 4 4 4 4 4 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 # index variable

# Model Specification:
fit1 <- map(
alist(
kcal.per.g ~ dnorm( mu , sigma ) ,
mu <- a[clade_id] ,
a[clade_id] ~ dnorm( 0.6 , 10 ) ,
sigma ~ dunif( 0 , 10 )
 ) ,
data = d )

Answer 1

Наиболее аналогичным способом сделать это с помощью пакета rstanarm является

library(rstanarm)
fit1 <- stan_glmer(kcal.per.g ~ 1 + (1 | clade_id), data = milk,
                   prior_intercept = normal(0.6, 1, autoscale = FALSE), 
                   prior_aux = exponential(rate = 1/5, autoscale = FALSE),
                   prior_covariance = decov(shape = 10, scale = 1))

Однако это не совсем то же самое по следующим причинам:

1. Ограниченные одинаковые приоры на sigma не реализованы, потому что они не очень хорошая идея, поэтому я использовал экспоненциальное распределение с ожиданием 5 вместо
2. Исправление стандартного отклонения на a также не реализован, поэтому я использовал гамма-распределение с ожиданием 10
3. Иерархические модели в rstanarm (и lme4) параметризуются с отклонениями от общих параметров, вместо того, чтобы использовать ожидание 0,6 для a Я использовал ожидание 0,6 для глобального перехвата и предыдущего a нормально с ожиданием ноль. Это означает, что вам нужно позвонить coef(fit1) скорее, чем ranef(fit1) чтобы увидеть "перехватывает", как они параметризованы в исходной модели.