GLMM(使用 lme4)中的随机效应总和应接近于零的程度

机器算法验证 r lme4-nlme 随机效应模型 咕噜咕噜
2022-03-24 12:51:51

我正在使用lme4R 中的包进行一些逻辑混合效果建模。
我的理解是每个随机效应的总和应该为零。

当我使用 制作玩具线性混合模型lmer时,随机效应通常 <证实了我的信念, 但在玩具二项式模型(以及我的真实二项式数据模型)中,一些随机效应总和为 ~ 0.9。1010colSums(ranef(model)$groups) ~ 0

我应该担心吗?我如何解释这个?

这是一个线性玩具示例

toylin<-function(n=30,gn=10,doplot=FALSE){
 require(lme4)
 x=runif(n,0,1000)
 y1=matrix(0,gn,n)
 y2=y1
 for (gx in 1:gn)
 {
   y1[gx,]=2*x*(1+(gx-5.5)/10) + gx-5.5  + rnorm(n,sd=10)
   y2[gx,]=3*x*(1+(gx-5.5)/10) * runif(1,1,10)  + rnorm(n,sd=20)
 }
 c1=y1*0;
 c2=y2*0+1;
 y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),]))
 g=rep(1:gn,each=n,times=2)
 x=rep(x,times=gn*2)
 c=c(c1,c2)
 df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g)))
 (m=lmer(y~x*c + (x*c|g),data=df))
 if (doplot==TRUE)
  {require(lattice)
   df$fit=fitted(m)
   plot1=xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1)
   plot2=xyplot(y ~ x|g,data=df,group=c)
   print(plot1+plot2)
  }
 print(colMeans(ranef(m)$g))
 m
}

在这种情况下,colMeans 总是出现<106

这是一个二项式玩具示例(我会分享我的实际数据,但它正在提交出版,我不确定期刊的事先发布政策是什么):


toybin<-function(n=100,gn=4,doplot=FALSE){
  require(lme4)
x=runif(n,-16,16) y1=matrix(0,gn,n) y2=y1 for (gx in 1:gn) { com=runif(1,1,5) ucom=runif(1,1,5) y1[gx,]=tanh(x/(com+ucom) + rnorm(1)) > runif(x,-1,1) y2[gx,]=tanh(2*(x+2)/com + rnorm(1)) > runif(x,-1,1) } c1=y1*0; c2=y2*0+1; y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),])) g=rep(1:gn,each=n,times=2) x=rep(x,times=gn*2) c=c(c1,c2) df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g))) (m=lmer(y~x*c + (x*c|g),data=df,family=binomial)) if (doplot==TRUE) {require(lattice) df$fit=fitted(m) print(xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1)) } print(colMeans(ranef(m)$g)) m }

现在 colMeans 有时会高于 0.3,而且绝对高于线性示例的平均值。

1个回答

由于@Hemmo 的代码在“Bounty”框中略有损坏,因此我将此重新格式化的版本添加为“community wiki”。如果这不是对 wiki 的适当使用,我提前道歉。随意删除它。

library(mvabund)
library(lme4) 
data(spider) 
Y <- as.matrix(spider$abund)
X <- spider$x 
X <- X[ ,c(1, 4, 5, 6)] 
X <- rbind(X, X, X, X, X, X, X, X, X, X, X, X) 
site <- rep(seq(1, 28), 12) 
dataspider <- data.frame(c(Y), X, site) 
names(dataspider) <- c("Y","soil.dry", "moss", "herb.layer", "reflection", "site") 
fit <- glmer(
  Y ~ soil.dry + moss + herb.layer + reflection + (1|site), 
  family = poisson(link = log), 
  data = dataspider,
  control = glmerControl(optimizer = "bobyqa")
)