J'utilise le lme4
package dans R pour faire de la modélisation logistique à effets mixtes.
J'ai cru comprendre que la somme de chaque effet aléatoire devait être nulle.
Lorsque je crée des modèles mixtes linéaires de jouets à l'aide de lmer
, les effets aléatoires sont généralement < confirmant ma conviction que le
Mais dans les modèles binomiaux jouets (et dans les modèles de mes données binomiales réelles) une partie de la somme des effets aléatoires à ~ 0.9.colSums(ranef(model)$groups) ~ 0
Dois-je m'inquiéter? Comment interpréter cela?
Voici un exemple de jouet linéaire
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
}
Dans ce cas, les colMeans sortent toujours
Voici un exemple de jouet binomial (je partagerais mes données réelles, mais elles sont soumises pour publication et je ne suis pas sûr de ce que la politique de la revue est en train de publier au préalable):
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
}
Maintenant, les colMeans dépassent parfois 0,3, et nettement plus, en moyenne que l'exemple linéaire.
la source
Réponses:
Depuis que le code de @ Hemmo a été légèrement altéré dans la case "Bounty", j'ajoute cette version reformatée en "wiki communautaire". Si ce n'est pas une utilisation appropriée du wiki, je m'excuse à l'avance. N'hésitez pas à le retirer.
la source