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.
>
>