7
$\begingroup$

in GLMM faq page http://glmm.wikidot.com/faq there is a statement about overfitting:

"One alternative (suggested by Robert LaBudde) is to "fit the model with the random factor as a fixed effect, get the level coefficients in the sum to zero form, and then compute the standard deviation of the coefficients." This is appropriate for users who are (a) primarily interested in measuring variation (i.e. the random effects are not just nuisance parameters, and the variability [rather than the estimated values for each level] is of scientific interest)"

What exactly does this mean in practice, doe sit mean using offset() for the random effects that are now fixed effects?

Also is there a way to diagnose overfitting other than very small random effect variance components? For instance how do you get the df for a GLMM that are reported in papers?

$\endgroup$
6
$\begingroup$

Computing the standard deviation (or variance of the coefficients) essentially means getting the fixed-effect estimates for each level (analogous to the BLUPs/conditional modes in a mixed model) and computing their variance. You can do this by appropriately setting contrasts to contr.sum (sum-to-zero contrasts) (in this case you'll still have to reconstruct the value of one level, since the model will only fit n-1 coefficients in a model with an intercept), and/or appropriate use of -1 or +0 in the model to fit a no-intercept model where the coefficients are computed for every level. Or, as shown below, you can just use brute force via predict (or e.g. via the lsmeans package) to compute values for each level ...

Make up data with only two levels of the RE grouping variable:

dd <- expand.grid(f1=factor(1:3),f2=factor(1:2),rep=1:10)
library(lme4)
simList <- 
   suppressMessages(simulate(~f1+(1|f2),
                 newdata=dd,
                 family="gaussian",
                 newparams=list(theta=1,beta=c(0,1,2),sigma=1),
                 seed=101,n=500))

Fit f2 as a random effect and retrieve estimated variance:

sumfun1 <- function(y0) {
    m <- lmer(y~f1+(1|f2),data=transform(dd,y=y0))
    unlist(VarCorr(m))
}

library(plyr)
r1 <- laply(simList,sumfun1,.progress="text")

This actually works surprisingly well given the small number of levels:

mean(r1)  ## 0.98
confint(lm(r1~1)) 
##                 2.5 %   97.5 %
## (Intercept) 0.9248779 1.189029

But we often get zero estimates of the variance:

sum(r1==0)  ## 60

(and a handful of very small values)

sum(log10(r1)<(-6))  ## 69

Now try it via fixed effects:

sumfun2 <- function(y0) {
    lm1 <- lm(y~f1+f2,data=transform(dd,y=y0))
    pframe <- data.frame(f1="1",f2=levels(dd$f2))
    var(predict(lm1,newdata=pframe))
}

r2 <- laply(simList,sumfun2,.progress="text")
mean(r2)  ## 1.01294
confint(lm(r2~1))
##             2.5 %   97.5 %
## (Intercept) 0.89081 1.135071

r1[log10(r1)< (-6)] <- 1e-6
p0 <- rbind(data.frame(m="f1=random",r=r1),
            data.frame(m="f1=fixed",r=r2))
library(ggplot2); theme_set(theme_bw())
ggplot(p0,aes(x=log10(r),fill=m))+
    geom_histogram(alpha=0.5,position="identity")+
        geom_vline(xintercept=0,lty=2)

enter image description here

The fixed-effect approach actually works better than I expected ...

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.