用混合效应模型预测R的概率

Question

我使用R函数glmer定义了一个二元响应混合效应模型，如下所示:

fit <-glmer(binary_r ~ cat1 + (1 | SUBJECTIDf) + (1 | cat2) + (1 | cat1:cat2), 
                                  family = binomial("logit"), data = mydata))

式中，cat1和cat2为分类变量，SUBJECTIDf为标记研究对象个体的因子变量。SUBJECTIDf和cat2是交叉分类因子。

我想对上述模型做以下操作:

1. 创建一个表，该表表示属于cat1和cat2的每个类别组合的积极响应的概率;

2. 创建一个图(可能是毛毛虫图)，显示a)下定义的概率;

和

为第1条定义的概率创建汇总统计。，包括cat1和cat2所有类别组合的最小和最大概率。

我假设为了实现上述结果，假设在1下定义的单个概率是合适的。，为cat1和cat2的给定组合中所有科目的平均值或中位数。

Answer 1

既然你没有给出一个可复制的例子，我将模拟一个…这部分答案只是设置了一个示例数据集。

## crossed cat2 and SUBJECTIDf
dd <- expand.grid(cat2=factor(letters[1:10]),
                  SUBJECTIDf=factor(1:10))
## each subject gets one value of cat1:
## for example, half get A while half get B
cat_tab <- data.frame(SUBJECTIDf=factor(1:10),cat1=rep(c("A","B"),5))
dd <- merge(dd,cat_tab)

我在适应cat1和cat2不是交叉分类的设计考虑方面遇到了一点麻烦。我将删掉一些类别:

dd <- with(dd,dd[!(cat1=="A" & cat2 %in% c("a","b","c")),])

现在我们已经设置了设计，我们模拟响应值:

library(lme4)
form <- binary_r ~ cat1 + (1 | SUBJECTIDf) + (1 | cat2) + (1 | cat1:cat2)
dd$binary_r <- simulate(form[-2],  ## RHS only
         family=binomial,
         newdata=dd,
         newparams=list(beta=0:1,
                        theta=c(2,4,1)),
         seed=101)[[1]]

在这一点上，我们得到了你上面建议的模型拟合。

fit <- glmer(form, family = binomial, data=dd)

创建一个表，其中表示对属于cat1和cat2的每个类别组合的正面响应的概率;

通过设置下面的re.form来排除主题ID，我们隐式地计算假设的中位数个体的值(即随机效应设为零;平均值和中位数预测在logit尺度上是一致的，但一旦我们反变换到概率尺度就不一致了)。

获取数据中cat1和cat2的唯一组合:

newdd <- unique(dd[,c("cat1","cat2")])
newdd$SUBJECTIDf <- NA  ## need to have SUBJECTIDf in the data frame ...
t1 <- predict(fit,newdata=newdd,type="response",
              re.form=~(1|cat2)+(1|cat1:cat2))
newdd <- data.frame(newdd[,c("cat1","cat2")],pred=t1)
head(newdd)
##   cat1 cat2        pred
## 4    A    d 0.215336024
## 5    A    e 0.944897414
## 6    A    f 0.036751551
## 7    A    g 0.003819873
## 8    A    h 0.970115614
## 9    A    i 0.003819873

我们还可以计算所有个体的预测如下:

## we happen to have a factorial design, but expand.grid() would
## e.g. fill in missing values
newdd2 <- unique(dd[,c("cat1","cat2","SUBJECTIDf")])
t2 <- predict(fit,newdata=newdd2,type="response",
              re.form=NULL)
newdd2$pred <- t2
head(newdd2)

为了创建图表(见下文)，我们必须对类别组合中的个体进行汇总。

创建一个图(可能是毛毛虫图)，显示a)下定义的概率;

library(ggplot2); theme_set(theme_bw())
ggplot(newdd,aes(cat2,pred,colour=cat1))+
   geom_point()+scale_colour_brewer(palette="Set1")

或者，汇总特定主题的预测:

ggplot(newdd2,aes(cat2,pred,colour=cat1))+
   stat_summary(fun.y=mean,geom="point")+
   scale_colour_brewer(palette="Set1")

我们可以在cat2类别上使用reorder()来尝试获得更合理的顺序，但由于存在cat1:cat2交互，这可能不太有效。卡特彼勒图(即获得预测的不确定性)有点棘手，因为很难获得结合条件模式(单个随机效应的值)和固定效应的不确定性的预测的不确定性。可以通过(1)假设条件模式和固定效果是独立的或(2)参数引导(bootMer)来完成，但两者都比我现在愿意采取的麻烦多一点……

为1中定义的概率创建汇总统计。，包括cat1和cat2所有类别组合的最小和最大概率。

这对我来说没有意义，除非我们走分解路线。如果我们预测了cat1和cat2的每个组合，那么我们只有每个组合的单个值(即没有"最小/最大"概率)。在基数R中聚合很容易，例如

 aggregate(pred~cat1:cat2,data=newdd2,
            FUN=function(x) c(min=min(x),max=max(x)))

或

 library(dplyr)
 newdd2 %>% group_by(cat1,cat2) %>%
     summarise(min=min(pred),max=max(pred))