0
$\begingroup$

I want to estimate the fixed effects and the covariance matrix (or standard deviation of the random effects terms) for a GLMM with a probit link using glmer. Most of the documentation that I can find is regarding lmer and I don't understand the differences. Se two minimal example below where the corresponding LMM is unbiasedly estimated by lmer but GLMM is not by glmer.

The first example is a random intercept mixed model with fixed effets c(2.5, -0.1, 0.2), random intercept sd=2, residual sd=1. The output of glmer is on average c(3.98, -0.026, 0.054) for fixed effects and 4.29 for sdcor from VarCorr.

#Minimal example

nrep=10;N<-1000;n<-4  
Sglmer<-as.data.frame(matrix(nrow=nrep,ncol=4))   
Slmer<-as.data.frame(matrix(nrow=nrep,ncol=5)) 

for(i in 1:nrep){   
D<-data.frame(X1<-rep(rbinom(N,1,0.3),n),X2<-rep(rnorm(N),n),id<-rep(1:N,n))  
err<-rep(rnorm(N,0,2),n)+rnorm(n*N)   
D$zstar<-2.5-0.1*X1+0.2*X2+err    
D$z<-D$zstar>0
Tlmer<-lmer(zstar~X1+X2+(1|id),data=D)  
Tglmer<-glmer(z~X1+X2+(1|id),data=D,family=binomial(probit))
Slmer[i,4:5]<-as.data.frame(VarCorr(Tlmer))[1:2,'sdcor']  
Slmer[i,1:3]<-fixef(Tlmer)  
Sglmer[i,4]<-as.data.frame(VarCorr(Tglmer))["sdcor"]  
Sglmer[i,1:3]<-fixef(Tglmer) }

summary(Slmer) 
summary(Sglmer)

The second example is a simplification of what I really want to estimate, a GLMM with fixed intercept and random slope, with fixed effets c(2.5, -0.1, 0.2), random slope sd=2, residual sd=1. The output of glmer is on average c(3.15, -0.05, 0.32) fixed effects and 2.65 sdcor from VarCorr.

nrep=10;N<-1000;n<-4  
Sglmer2<-as.data.frame(matrix(nrow=nrep,ncol=4))  
Slmer2<-as.data.frame(matrix(nrow=nrep,ncol=5))  
for(i in 1:nrep){  
D<-data.frame(X1<-rep(rbinom(N,1,0.3),n),time<-c(rep(0,N),rep(2,N),rep(4,N),rep(6,N)),id<-rep(1:N,n))  
err<-rnorm(N,0,2)*D$time+rnorm(n*N,0,1)  
zstar<-2.5-0.1*X1+0.2*time+err  
z<-zstar>0  
Tlmer<-lmer(zstar~X1+time+(0+time|id),data=D)  
Tglmer<-glmer(z~X1+time+(0+time|id),data=D,family=binomial(probit))  
Slmer2[i,4:5]<-as.data.frame(VarCorr(Tlmer))[1:2,'sdcor']  
Slmer2[i,1:3]<-fixef(Tlmer)  
Sglmer2[i,4]<-as.data.frame(VarCorr(Tglmer))["sdcor"]  
Sglmer2[i,1:3]<-fixef(Tglmer)    }  

summary(Slmer2)
summary(Sglmer2)
New contributor
Minna is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$

Your Answer

Minna is a new contributor. Be nice, and check out our Code of Conduct.

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.