R help - Apr 2009 - Sean / Re: question related to fitting overdispersion count data using lmer quasipoisson

Hey Buddy,

Hope you have been doing well since last contact.

If you have the answer to the following question, please let me know.

If you have chance to travel up north. let me know.

best,

-Sean


---------- Forwarded message ----------
From: Sean Zhang <seanecon@gmail.com>
Date: Sat, Apr 11, 2009 at 12:12 PM
Subject: question related to fitting overdispersion count data using lmer
quasipoisson
To: r-help@r-project.org
Cc: seanecon@gmail.com


Dear R-helpers:
I have a question related to fitting overdispersed count data using lmer.
Basically, I simulate an overdispsed data set by adding an observation-level
normal random shock
into exp(....+rnorm()).
Then I fit a lmer quasipoisson model.
The estimation results are very off (see model output of fit.lmer.over.quasi
below).
Can someone kindly explain to me what went wrong?

Many thanks in advance.

-Sean

#data simulation (modified from code at
http://markmail.org/message/j3zmgrklihe73p4p)
set.seed(100)
m <- 5
n <- 100
N <- n*m
#X <- cbind(1,runif(N))
X <- cbind(1,rnorm(N))
X <- cbind(runif(N),rnorm(N))
id <- rep(1:n,each=m)
#
Z <- kronecker(diag(n),rep(1,m))
#Possion with group level heterogeneity
z <- rpois(N, exp(X%*%matrix(c(1,2)) + Z%*%matrix(rnorm(n))))
#2*rnorm(n*m) is added to each observation to create overdispersion
z.overdis <- rpois(N, exp(X%*%matrix(c(1,2)) + Z%*%matrix(rnorm(n)) +
2*rnorm(n*m)))

#without observation-level random shock i.e., 2*rnorm(n*m), estimate results
are very accurate
(fit.lmer <- lmer(z ~ X + (1|id), family="poisson",verbose=F))
 #Generalized linear mixed model fit by the Laplace approximation
#Formula: z ~ X + (1 | id)
# AIC BIC logLik deviance
# 851 868   -422      843
#Random effects:
# Groups Name        Variance Std.Dev.
# id     (Intercept) 0.977    0.988
#Number of obs: 500, groups: id, 100
#
#Fixed effects:
#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)  -0.0128     0.1116    -0.1      0.9
#X1            1.0615     0.0601    17.7   <2e-16 ***
#X2            2.0236     0.0214    94.7   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Correlation of Fixed Effects:
#   (Intr) X1
#X1 -0.349
#X2 -0.270  0.258




#Now you can see the results are very off
(fit.lmer.over.quasi <- lmer(z.overdis ~ X + (1|id),
family=quasipoisson(link="log"),verbose=F))
#Generalized linear mixed model fit by the Laplace approximation
#Formula: z.overdis ~ X + (1 | id)
#   AIC   BIC logLik deviance
# 41867 41888 -20929    41857
#Random effects:
# Groups   Name        Variance Std.Dev.
# id       (Intercept) 175.8    13.26
# Residual              72.9     8.54
#Number of obs: 500, groups: id, 100
#
#Fixed effects:
#            Estimate Std. Error t value
#(Intercept)   1.3530     1.3492    1.00
#X1            1.0834     0.2273    4.77
#X2            1.3501     0.0783   17.25
#
#Correlation of Fixed Effects:
#   (Intr) X1
#X1 -0.099
#X2 -0.055  0.070

	[[alternative HTML version deleted]]