Estou usando o lme4
pacote no R para fazer alguma modelagem logística de efeitos mistos.
Meu entendimento era que a soma de cada efeito aleatório deveria ser zero.
Quando eu uso modelos mistos lineares de brinquedo lmer
, os efeitos aleatórios geralmente são < confirmando minha crença de que alguns dos
modelos binomiais de brinquedos But (e nos modelos dos meus dados binomiais reais) alguns dos efeitos aleatórios somam ~ 0,9.colSums(ranef(model)$groups) ~ 0
Devo me preocupar? Como eu interpreto isso?
Aqui está um exemplo de brinquedo linear
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
}
Nesse caso, o colMeans sempre sai
Aqui está um exemplo de brinquedo binomial (eu compartilharia meus dados reais, mas eles estão sendo enviados para publicação e não tenho certeza de qual é a política do diário antes de postar):
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
}
Agora, o colMeans às vezes sai acima de 0,3 e, definitivamente, mais alto, em média do que o exemplo linear.
fonte
Respostas:
Como o código do @ Hemmo ficou um pouco distorcido na caixa "Bounty", estou adicionando esta versão reformatada como "wiki da comunidade". Se este não for um uso apropriado do wiki, peço desculpas antecipadamente. Sinta-se livre para removê-lo.
fonte