1

背景： McEearth（2016 年）在他的重新思考书籍第 158-159 页中，使用索引变量而不是虚拟编码来对称为“clade”的 3 类变量进行预测“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 )
4

1 回答 1

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 ）中的分层模型被参数化为与公共参数的偏差，因此a我没有使用 0.6 的期望值，而是对全局截距使用 0.6 的期望值，而先验a值是正常的，期望值为 0。这意味着您需要调用coef(fit1)而不是ranef(fit1)查看“截距”，因为它们在原始模型中已参数化。
于 2018-03-21T21:09:44.763 回答