[R] lme4 and mcmcamp

[Ivar Herfindal]

Hello.

I ran into a similar situation as you did, and according to D. Bates, 
this can be due to the inclusion of the offset variable. It appears that 
the mcmcsamp does not work very well with offset variables. For more 
details, have a look at this thread:

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2007q1/000061.html

I experimented a bit with xyplot on the mcmcsamp-object, and clearly, 
the procedure gets stuck at some points. This can explain why you get 
totally different results from mcmcsamp compared to the summary-command. 
Accordingly, it seems that the summary gives the most reliable results 
when a offset variable is included.

Ivar

estrain at postoffice.utas.edu.au wrote:

>Dear R users
>I am trying to obtain p-values for (quasi)poisson lmer models, using
>Markov-chain Monte Carlo sampling and the command summary.
>  
>
>>My problems is that p values derived from both these methods are
>>    
>>
>totally different. My question is
>(1) there a bug in my code and
>  
>
>
>(2) How can I proceed, left with these uncertainties in the estimations of
>  
>
>>the p-values?
>>
>>Below is the corresponding R code with some output:
>>##
>>fit<-lmer(End~Treatment+offset(log(Area))+(1|Site/Treatment), family=poisson
>>    
>>
>
>  
>
>>## Results
>>summary(fit)
>>    
>>
>
>AIC BIC loglik deviance
>28.92 35.99 -8.46 16.92
>Random effects
>
>Groups Name                    Variance    Std dev.
>Treatment * Site Intercept  5e-10        2.2361e-05
>Site                    Intercept  5e-10       2.2361e-05
>
>
>number of obs 24 groups Treatment*Site 8 and Site 2
>
>Fixed effects
>
>             Estimate   Std error  z value   P
>Intercept -3.8290   0.4995      -7.666   1.77e-14****
>Treatment 2 3.7970 0.505      7.516     5.51e-14****
>
>Treatment 3 0.2409 9.6704 0.359 0.719 Treatment 4 -0.2483 0.8661 -0.287 0.774
>
>Correlation of fixed effects
>
>    Intra T2 T3
>T2 -0.989
>T3 -0.745 0.737
>T4 -0.577 0.570 0.430  
>
>  
>
>>The p-values from mcmc are:
>>
>>    
>>
>mcmcpvalue<-function(samp)
>{
>std<-backsolve(chol(var(samp)),
>cbind(0,t(samp))-colMeans(samp),
>transpose=TRUE)
>sqdist<-colSums(std*std)
>sum(sqdist[-1]>sqdist[1]/nrow(samp)
>}
>
>fitSI<-mcmcsamp(fit,50000)
>library(coda)
>HPDinterval(fitSI)
>
>                   lower        upper
>
>
>Intercept -4.0778905 -3.1366836
>Treatment2 3.4455972 4.3196598
>Treatment 3 0.399302 1.287747
>Treatment 4 -1.7898933 -0.2980325
>log(Treatment*Site.(in)) -22.2198233 -19.7342530 
>log(Site.(In)) -28.7857601 -23.0952939
>
>mcmcpvalue(as.Matrix(fitSI[,1]))
>etc
>
>Intecept 0
>Treatment 2 0
>Treatment 3 0.0075
>Treatment 4 0.00758
>log(Treatment*Site.(in)) 0
>log(Site.(In)) 0
>
>Any help in explaining why these two results are completely different
>would be much appreciated. I have tried the example and read all of the posts.
>
>Cheers
>Beth
>
>______________________________________________
>R-help at stat.math.ethz.ch 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.
>  
>