[R-sig-ME] Variance estimation in gamma GLMMs

[Björn Johansson]

In the program below I am making comparisons between lme4::glmer and glmmTMB for some simple randomly generated data with just an intercept and a single random effect, for the distributions Poisson, binomial and gamma with log link. For the Poisson and binomial data the results are almost identical. For the gamma data they are very close, except for the estimates of the variance of the random effects (the actual variance is 0.15 in all cases). I would like to understand why.

Kind regards

-----------------------------------------------------------------

library(lme4)
library(glmmTMB)

# Poisson

beta <- 3  # Intercept
varv <- 0.15 # Variance of random effect
mlf <- sample(1:100,10000,replace=TRUE) # Levels of the random effect
simdata <- data.frame(mlf)
v <- rnorm(100,sd=sqrt(varv))
simdata$v <- v[simdata$mlf]
simdata$expx <- exp(beta + simdata$v)
simdata$x <- rpois(nrow(simdata),simdata$expx)

est_glmer <- glmer(x ~ 1 + (1|mlf), family = poisson, data = simdata, nAGQ=1L)
est_tmb <- glmmTMB(x ~ 1 + (1|mlf), data=simdata, family=poisson)
summary(est_glmer)$varcor$mlf[[1]]
unname(unlist(summary(est_tmb)$varcor))

# Binomial

beta <- 0.5 # Intercept
varv <- 0.15  # Variance of random effect
mlf <- sample(1:100,10000,replace=TRUE) # Levels for the random effect
simdata <- data.frame(mlf)
v <- rnorm(100,sd=sqrt(varv))
simdata$v <- v[simdata$mlf]
simdata$p <- exp(beta+simdata$v)/(1+exp(beta+simdata$v))
simdata$x <- rbinom(nrow(simdata),1,simdata$p)

est_glmer <- glmer(x ~ 1 + (1|mlf), family = binomial, data = simdata, nAGQ=1L)
est_tmb <- glmmTMB(x ~ 1 + (1|mlf), data=simdata, family=binomial)
summary(est_glmer)$varcor$mlf[[1]]
unname(unlist(summary(est_tmb)$varcor))

# Gamma

mlf <- sample(1:100,10000,replace=TRUE)
beta <- 7  # Intercept
phi <- 0.5 # 1/(shape parameter) in the gamma distribution
varv=0.15    # Variance of the random effect
simdata <- data.frame(mlf)
v <- rnorm(100,sd=sqrt(varv))
simdata$v <- v[simdata$mlf]
simdata$expx <- exp(beta + simdata$v)
simdata$x <- rgamma(nrow(simdata),1/phi,1/(simdata$expx*phi))

est_glmer <- glmer(x ~ 1 + (1|mlf), family = Gamma(link="log"), data = simdata, nAGQ=1L)
est_tmb <- glmmTMB(x ~ 1 + (1|mlf), data=simdata, family=Gamma(link="log"))
summary(est_glmer)$varcor$mlf[[1]]
unname(unlist(summary(est_tmb)$varcor))


Skickades från E-post för Windows


	[[alternative HTML version deleted]]