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]]
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
