完全分离-如何在bglmer中对固定效果施加“零均值先验先验”

问题描述

我的glmer模型包含两个预测变量和一个交互项，因此完全分离。遵循Ben Bolker的建议here和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)))

因子1和因子2都是因子变量，每个变量具有四个级别。对于我的代码，我遵循示例here。据我了解，我现在将固定效应结构的所有元素的SD均为3的零均值普通先验。

该代码似乎有效，但是我完全不确定我所做的是否正确。一般建议使用3 SD来帮助完全分离吗？以及我该如何指定fixef.priors仅在交互项上使用？ （完全分隔涉及Factor1和Factor2的特定组合，通常不涉及Factor1或Factor2）。还是要考虑相互作用，我是否必须对这三个要素都施加固定的先验先验？

解决方法

如有疑问，请尝试。 tl; dr 最好对所有内容进行惩罚，但是如果您喜欢对交互进行惩罚而不对主要效果进行惩罚，则可以这样做（或者几乎可以；您必须提供有限的sd，但是它可能非常大） 。 sd = 3是合理的； sd = 2.5是截距以外的其他参数的默认值。除非您要拟合预测模型并想进行交叉验证以选择最佳惩罚强度，否则就没有真正的自动决策方法；除非您使用预测模型，否则不做任何决定。您只需要选择一个使所有参数“合理”的sd，而不会将确定好的参数过度压缩为零。

我模拟了像您一样的数据（2个因子，每个因子有4个等级，二进制响应，一个随机拦截项），并尝试了不同的惩罚方案。我设置了真实的参数，以便大多数参数都是合理的（beta = 1），但是有一个很大的主效应和一个很大的交互参数（beta = 12）。

• 不受到处罚（glmer）
• 用sd=3进行惩罚
• 使用sd=1惩罚所有术语
• 主效应项（sd=1e4）的惩罚很弱，而交互作用（sd=2的惩罚是中等的[[mixed]]

来自help("bmerDist-class")：

指定标准时 偏差，长度小于固定数量的向量 效果将重复其尾巴，而第一个元素 假定仅适用于拦截项。所以在 默认值为“ c（10，2.5）”，则拦截器会收到一个标准 10的偏差和各种斜率都给出了 标准偏差为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,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))