1

我想拟合 GLMM Poisson 模型计数。我有 121 个受试者 ( subject)，我观察到每个受试者 8 个泊松计数 ( count)：它们对应于 2 种类型的事件 ( event) x 4 个句点 ( period)。

'data.frame':   968 obs. of  4 variables:
 $ count  : num  4 0 2 3 3 0 1 0 8 14 ...
 $ subject: num  1 1 1 1 1 1 1 1 2 2 ...
 $ event  : Factor w/ 2 levels "call","visit": 2 2 2 2 2 2 2 2 1 1 ...
 $ period : Factor w/ 4 levels "period_1","period_2",..: 1 2 3 4 1 2 3 4 1 2 ...

我想要的是：

  • 我想用贝叶斯方法拟合 GLMM，假设 (1) subject， (2) subject:visit， (3)的随机效应subject:event。
  • 我想为所有固定效果先设置 N(0, 10^6)，
  • 我想分别为 (1) 、 (2) 、 (3)的随机效应设置 N(0, sigma2_a ), N(0, sigma2_b ), N(0, sigma2_c ) 先验，subjectsubject:visitsubject:event
  • 我想分别为sigma2_a、sigma2_b、sigma2_c设置统一的先验。

我设法得到的：

  • 我相信我正在正确设置模型公式，并且正在为固定效应参数设置所需的先验：

    ## define some priors          
prior <- c(prior_string("normal(0,10^3)", class = "b"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_2"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_3"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_4"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_2"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_3"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_4"))

## fit model
fit1 <- brm(count ~ event + period + event:period + (1|subject) + (0 + event|subject) + (0 + period|subject),
    data = data.long, family = poisson(),
    prior = prior,
    warmup = 1000, iter = 4000, chains = 4,
    cores = 7)

我挣扎的是：

  • 如何分别为 (1) 、 (2) 、 (3)的随机效应设置 N(0, sigma2_a ), N(0, sigma2_b ), N(0, sigma2_c ) 先验，subjectsubject:visitsubject:event
  • 如何分别为sigma2_a、sigma2_b、sigma2_c设置统一的先验。
4

1 回答 1

1

关于你想要什么：你能否更详细地解释你想要subject:visit什么subject:event？目前，我无法判断您的模型是否为您的目的正确指定。

关于先验：

  1. 如果设置prior_string("normal(0,10^3)", class = "b")，则无需为每个系数再次指定该先验。为整个先前的课程全局设置它就足够了。请注意，鉴于您的响应变量和预测变量的规模， normal(0, 1000) 是如此之宽，以至于它实际上没有影响。你也可以把它排除在外。
  2. 这些normal(0, sigma_x)先验是为所有随机效应项自动设置的。无需担心它们。
  3. 您可以使用 设置随机效应 SD 的先验class = "sd"。例如prior_string("cauchy(0, 5)", class = "sd"). 我明确建议您不要使用统一先验，特别是因为它们对参数（标准差）有一个硬上限，它没有理论上限。
于 2018-05-08T22:04:57.537 回答