二項式模型中隨機效應的估計 (lme4)

我正在隨機模擬伯努利試驗組之間，然後我將相應的模型與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 估計值接近總體比例 (在我的例子中）。增加後估計和比例之間的差異可以忽略不計，但人們總能得到這張照片。這是一個已知的事實，為什麼它成立？我希望得到以經驗比例為中心的估計。

您所看到的是一種稱為收縮的現象，這是混合模型的基本屬性；作為每個估計值的相對方差的函數，個體組估計值向整體平均值“縮小”。（雖然 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)

引用自：https://stats.stackexchange.com/questions/142810