完全分离 - 如何在 bglmer 中施加 'zero-mean Normal priors on fixed effects'

Complete separation - how to impose 'zero-mean Normal priors on fixed effects' in bglmer

我的 glmer 模型包含两个预测变量和一个交互项,存在完全分离的问题。按照 Ben Bolker 的建议 here and here,然后我用 bglmer 拟合模型,对固定效应施加零均值正态先验。我的代码如下:

bglmer(Binary_outcome ~ (1|Subject) + Factor1 + Factor2 + Factor1:Factor2, 
       mydata, 
       control=glmerControl(optimizer="bobyqa"), 
       family = binomial, 
       fixef.prior = normal(sd = c(3, 3, 3)))

Factor1和Factor2都是因子变量,各有四个水平。对于我的代码,我遵循示例 here。据我所知,我现在将 SD 为 3 的零均值正态先验放在我的固定效应结构的所有元素上。

代码似乎有效,但我完全不确定我所做的实际上是否正确。 3 SD 是帮助完全分离的一般建议吗?我将如何指定 fixef.priors 仅在交互项上进行? (完全分离涉及 Factor1Factor2 的特定组合,而不是一般的 Factor1Factor2。或者,如果涉及交互,我是否必须对所有三个元素都设置固定效应先验?

如有疑问,请进行实验。 tl;dr 惩罚一切可能没问题,但如果你喜欢惩罚相互作用而不惩罚主要影响(或几乎;你必须给出一个有限的 sd 但它可以是很大)。 sd=3 合理; sd=2.5 是截距以外参数的默认值。除非你正在拟合一个预测模型并且想不厌其烦地做 cross-validation 来选择最佳的惩罚强度,否则并没有真正自动的决定方式;您只需要选择一个 sd,使您的所有参数都“合理”,而不会将 well-determined 参数过度压缩到零。

我模拟的数据有点像你的(2 个因素,每个因素有 4 个水平,二元响应,一个随机截距项)并尝试了不同的惩罚场景。我设置了真实参数,以便大多数参数都是合理的 (beta=1),但有一个大的主效应和一个大的交互参数 (beta=12)。

  • 无惩罚(glmer
  • 惩罚sd=3
  • 所有条款的惩罚 sd=1
  • 主效应项的惩罚很弱(sd=1e4),交互作用的惩罚中等(sd=2)[“混合”]

来自 help("bmerDist-class"):

When specifying standard deviations, a vector of length less than the number of fixed effects will have its tail repeated, while the first element is assumed to apply only to the intercept term. So in the default of ‘c(10, 2.5)’, the intercept receives a standard deviation of 10 and the various slopes are all given a standard deviation of 2.5

library(lme4)
library(blme)
library(broom.mixed)
library(dotwhisker)
library(colorspace)

dd <- expand.grid(Subject=factor(1:10),
                  Factor1=letters[1:4],
                  Factor2=LETTERS[1:4],
                  rep=1:20
                  )
bvec <- rep(1,16)
bvec[c(2,10)] <- 12
form <- response ~ (1|Subject) + Factor1 + Factor2 + Factor1:Factor2
set.seed(101)
dd$response <- simulate(form[-2],
         newdata=dd,
         newparams=list(beta=bvec, theta=1),
         family=binomial,
         weights=rep(1,nrow(dd)))[[1]]

b0 <- glmer(form, 
            dd,
            family = binomial)
bfun <- function(sd) {
   bglmer(form, dd, family = binomial, fixef.prior = normal(sd = sd))
}
b1 <- bfun(3)
b2 <- bfun(1)
## eight intercept + main effects first, then eight interaction parameters
b3 <- bfun(rep(c(1e4,2),c(8,8))

theme_set(theme_bw())
dwplot(list(unpenalized=b0,sd3=b1,sd1=b2,mixed=b3),effect="fixed") +
    coord_cartesian(xlim=c(-5,5))+
    geom_vline(xintercept=c(0,1),lty=2,colour="darkgray") +
    scale_colour_discrete_qualitative(guide=guide_legend(reverse=TRUE))