

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

