机器算法验证 - 二项式模型中随机效应的估计 (lme4) - 吾爱随笔录

二项式模型中随机效应的估计 (lme4)

机器算法验证混合模式广义线性模型二项分布随机效应模型 lme4-nlme

2022-03-05 15:45:24

我在组之间用随机模拟伯努利试验，然后我用相应的模型拟合包装： $\text{logit}\, \theta \sim {\cal N}(\text{logit}\, \theta_0, 1^2)$ lme4

library(lme4)
library(data.table)
I <- 30 # number of groups
J <- 10 # number of Bernoulli trials within each group
logit <- function(p) log(p)-log(1-p)
expit <- function(x) exp(x)/(1+exp(x))
theta0 <- 0.7
ddd <- data.table(subject=factor(1:I),logittheta=rnorm(I, logit(theta0)))[, list(result=rbinom(J, 1, expit(logittheta))), by=subject]
fit <- glmer(result~(1|subject), data=ddd, family="binomial")
props <- ddd[, list(p=mean(result)), by=subject]$p
estims <- expit(coef(fit)$subject[,1])
par(pty="s")
plot(props, estims, asp=1, xlim=c(0,1), ylim=c(0,1),
     xlab="proportion", ylab="glmer estimate")
abline(0,1)

然后我将各组的成功比例与他们的估计进行比较，我总是得到这样的结果：

在此处输入图像描述

通过“总是”，我的意思是 glmer 估计总是高于小比例的经验比例，总是低于高比例。glmer 估计值接近总体比例（在我的示例中在增加之后，估计值和比例之间的差异变得可以忽略不计，但人们总能得到这张照片。这是一个已知的事实，为什么它成立？我希望得到以经验比例为中心的估计。 $0.7$ $J$

1个回答

您看到的是一种称为收缩的现象，这是混合模型的基本属性；作为每个估计值的相对方差的函数，单个组估计值向整体平均值“缩小”。（虽然 CrossValidated 上的各种答案都讨论了收缩，但大多数都指的是套索或岭回归等技术；这个问题的答案提供了混合模型和其他收缩观点之间的联系。）

可以说收缩是可取的；它有时被称为借贷强度。特别是当我们每组的样本很少时，每组的单独估计将不如利用每个群体的一些汇集的估计精确。在贝叶斯或经验贝叶斯框架中，我们可以将人口水平分布视为群体水平估计的先验。当每组的信息量（样本大小/精度）变化很大时（例如，在人口非常少和非常大的区域的空间流行病学模型中），收缩估计特别有用/强大.

收缩属性应该适用于贝叶斯和常客拟合方法——方法之间的真正区别在于顶层（常客的“惩罚加权残差平方和”是贝叶斯在组级别的对数后验偏差...... ) 下图中显示lme4和MCMCglmm结果的主要区别在于，由于 MCMCglmm 使用随机算法，具有相同观察比例的不同组的估计值略有不同。

通过更多的工作，我认为我们可以通过比较组和整个数据集的二项式方差来计算出预期的精确收缩程度，但同时这里有一个演示（事实上 J=10 的情况看起来更少我认为比 J=20 缩小只是抽样变化）。（我不小心将模拟参数更改为均值 = 0.5，RE 标准偏差 = 0.7（在 logit 标度上）...）

在此处输入图像描述

library("lme4")
library("MCMCglmm")
##' @param I number of groups
##' @param J number of Bernoulli trials within each group
##' @param theta random effects standard deviation (logit scale)
##' @param beta intercept (logit scale)
simfun <- function(I=30,J=10,theta=0.7,beta=0,seed=NULL) {
    if (!is.null(seed)) set.seed(seed)
    ddd <- expand.grid(subject=factor(1:I),rep=1:J)
    ddd <- transform(ddd,
                     result=suppressMessages(simulate(~1+(1|subject),
                     family=binomial,
                     newdata=ddd,
                     newparams=list(theta=theta,beta=beta))[[1]]))
}
sumfun <- function(ddd) {
    fit <- glmer(result~(1|subject), data=ddd, family="binomial")
    fit2 <- MCMCglmm(result~1,random=~subject, data=ddd,
                    family="categorical",verbose=FALSE,
                    pr=TRUE)
    res <- data.frame(
        props=with(ddd,tapply(result,list(subject),mean)),
        lme4=plogis(coef(fit)$subject[,1]),
        MCMCglmm=plogis(colMeans(fit2$Sol[,-1])))
    return(res)
}
set.seed(101)
res <- do.call(rbind,
        lapply(c(10,20,50,100,500),
               function(J) {
                   data.frame(J=J,sumfun(simfun(J=J)))
               }))
library("reshape2")
m <- melt(res,id.vars=c("J","props"))
library("ggplot2"); theme_set(theme_bw())
ggplot(m,aes(props,value))+
    geom_point(aes(colour=factor(J),shape=variable))+
    geom_abline(intercept=0,slope=1,colour="gray")+
      labs(x="observed proportion",y="estimate")
ggsave("shrinkage.png",width=5,height=5)

其它你可能感兴趣的问题

上一篇在高斯混合中找到最大概率点下一篇如何解释 Harvey-Collier 检验和 Rainbow 检验的线性方向？