dear R-users,
I am trying to understand what is the
sigma parameter returner by glmer
I thought it was (an estimate of) the sigma parameter
defined by Mc Cullagh & Nelder (e.g. p 126 of 2nd edition)
but I ran some simulations and it seems that this is
something else.
I simulated data corresponding to a binomial model,
intended to be fitted by this command:
glmer(cbind(success,failure)~X+(1|group),family=binomial)
but I instead fitted the following model:
glmer(cbind(success,failure)~X+(1|group),family=quasibinomial)
(and repeated this process 500 times)
I expected sigma to be close to 1 but I found
that the mean sigma was about 0.05 (sd = 0.003)
If I do the analogous simulation study with glm,
that is, I simulate binomial data and fit them
with family=quasibinomial instead of binomial,
I find a mean dispersion parameter = 0.9999
(sd=0.09).
changing parameters values does not alter this
pattern.
In both cases, the fixed effects parameters are
correctly estimated.
here is the function I used to simulate data (taking
0 as the standard deviation <sigmag> of random effect
provides data suitable to glm)
function(x,theta,sigmag,nb.groups=10,size=50)
#-----------------------------
# sim.data.mixed
#-----------------------------
# simulates data for glmer
# Y is Binom(p,size)
# with logit(p) = theta1 + theta2*X + B
# where B is Norm(0,sigmag)
#-----------------------------
# x : the x values (same for each group)
# length(x) is the number of observations per group
# theta: the fixed effects parameters (intercept & slope)
# sigmag : the random effects standard deviation
# size : the binomial parameter (same for everybody)
{
group<-rep(1:nb.groups,rep(length(x),nb.groups))
random.effect<-rnorm(nb.groups,mean=0,sd=sigmag)
xmat<-expand.grid(x,random.effect)
eta<-theta[1]+theta[2]*xmat[,1]+xmat[,2]
y<-rbinom(length(eta),size=size,prob=invlogit(eta))
return(data.frame(success=y,failure=size-y,x=xmat[,1],group=group,b=xmat[,2]))
#------------------------------
# b (random effects) is returned here but not used by glmer
}