[R-sig-ME] lme4, glmmPQL and calculating QAICc

[Antonio.Gasparrini at lshtm.ac.uk]

Dear Ben and R users,
 
I put the pdSymm constructor in glmmPQL just because I tried also different correlation structures.
If I change it to the simpler 'random=~time|region' I got slightly different results and the notification of 'General positive-definite, Log-Cholesky parametrization' instead of the simple 'General positive-definite'.
 
Anyway, nothing affecting the results substantially.
I can't include the original data, so I simulated them (trends of Poisson-distributed counts with overdispersion for 20 regions) with parameters similar to the estimates I found. Here is the code:
 
##################################################################

set.seed(13041975)
library(MASS);library(lattice)

# BUILD THE DESIGN MATRIX FOR FIXED AND RANDOM EFFECTS
region <- rep(1:20,each=60)
time <- rep(1:60,20)
X <- cbind(rep(1,20*60), time)
Z <- diag(rep(1,20)) %x% cbind(rep(1,60),1:60)

# SET FIXED AND RANDOM EFFECTS (FROM MV-NORMAL)
fixef <- c(-6.5,0.001)
ranef <- mvrnorm(20,c(0,0),matrix(c(0.14^2,0,0,0.003^2),2,2))

# POPULATION (DIFFERENT SIZE)
pop <- rep(runif(20,10^4,10^6),each=60)

# PREDICTED
pred <- exp(X%*%fixef + Z%*%as.vector(t(ranef))) * pop

# OBSERVED (FROM NEGATIVE BINOMIAL WITH OVERDISPERSION PHI)
phi <- 1.4^2
obs <- rnbinom(length(pop), mu=pred, size=pred/(phi-1))

# PLOT (RATES)
xyplot(obs/pop*10^5+pred/pop*10^5~time|region,type=c("p","l"),
 distribute.type=T)

# MODELS

library(nlme)
pql <- glmmPQL(obs ~ offset(log(pop)) + time, random=~time|region,
 family=poisson)
pql$sigma
# THIS IS OK

library(lme4)
glmer <- glmer(obs ~ offset(log(pop)) + time + (time|region),
 family=quasipoisson)
glmer
lme4:::sigma(glmer)
# STRANGE RESULTS

##################################################################

After this results, I think there's something odd with glmer and quasipoisson family.
I'm currently running an analysis on real data, I hope someone can give me some hints.

Thanks

Antonio Gasparrini
Public and Environmental Health Research Unit (PEHRU)
London School of Hygiene & Tropical Medicine
Keppel Street, London WC1E 7HT, UK
Office: 0044 (0)20 79272406 - Mobile: 0044 (0)79 64925523

>>> Ben Bolker <bolker at ufl.edu> 12/02/2010 03:59 >>>
Antonio.Gasparrini at lshtm.ac.uk wrote:
> Dear R users,
>  
> following the recent discussion on calculating QAICc in lme4, I report the weird results I got comparing glmer and glmmPQL.
>  
> I ran these models: 
>  
> pql.model <- glmmPQL(outcome ~ offset(log(pop)) + time + 
>  harmonic(month,3,12), random=list(region=pdSymm(~time)), family=poisson, data)
>  
> glmer.model <- glmer(outcome ~ offset(log(pop)) + time + 
>  harmonic(mm,3,12) + (time|region), family=poisson, data)
> The first model with glmmPQL estimates a sigma (within-group error) anyway, both with poisson or quasipoisson family.
> Its value is 1.40
> The value of sigma^2 is equal to the overdispersion parameter of simpler glm-gam models (~1.96), which makes sense.
>  
> The second model with glmer doesn't estimate a sigma (correctly), but when the family is set to quasipoisson the estimated sigma [lme4:::sigma(glmer.model)] is 15.8, which is simply unbelievable. The standard errors are therefore incredibly huge.
>  
> I couldn't find a reason for that.
> Any comment/suggestion is more than welcome.
> Thanks

  This is interesting.  Can you provide a reproducible example?  Is
random=list(region=pdSymm(~time)) really equivalent to (time|region)?