
Dear R users,
Please ignore my first email titled Lme and p values it was an unfinished 
email. I am trying to calculate p values for lmer. I have read the posts on 

http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests&s=lme%20and%20aov.My question is why are the pvalue for my model so different between
summary and the mcmcpvalue commands. Please find my code below. Any help
would be much appreciated.

Thanks in advance.
Beth

My model

library(lme4)
fit<-lmer(End~Treatment+offset(log(Area))+(1|Site/Treatment),family=poisson)

> summary(fit)
Generalized linear mixed model fit using Laplace
Formula: End ~ T + offset(log(Area)) + (1 | S/T)
 Family: poisson(log link)
   AIC   BIC logLik deviance
 28.92 35.99  -8.46    16.92
Random effects:
 Groups Name        Variance Std.Dev.
 T:S    (Intercept) 5e-10    2.2361e-05
 S      (Intercept) 5e-10    2.2361e-05
number of obs: 24, groups: T:S, 8; S, 2

Estimated scale (compare to  1 )  0.7792541

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -3.8290     0.4995  -7.666 1.77e-14 ***
T2            3.7970     0.5050   7.519 5.51e-14 ***
T3            0.2409     0.6704   0.359    0.719
T4           -0.2483     0.8661  -0.287    0.774
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
   (Intr) T2     T3
T2 -0.989
T3 -0.745  0.737
T4 -0.577  0.570  0.430


 ##To find pvalues
> library(coda)
> 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)
> HPDinterval(fitSI)
                    lower       upper
(Intercept)    -4.0778905  -3.1866836
T2              3.4455972   4.3196598
T3              0.3993082   1.2877477
T4             -1.7898933  -0.2980325
log(T:S.(In)) -22.2198233 -19.7342530
log(S.(In))   -28.7857601 -23.0952939
attr(,"Probability")
[1] 0.95
> mcmcpvalue(as.matrix(fitSI[,1]))
[1] 0
> mcmcpvalue(as.matrix(fitSI[,2]))
[1] 0
> mcmcpvalue(as.matrix(fitSI[,3]))
[1] 0.0075
> mcmcpvalue(as.matrix(fitSI[,4]))
[1] 0.00758
> mcmcpvalue(as.matrix(fitSI[,5]))
[1] 0
> mcmcpvalue(as.matrix(fitSI[,6]))
[1] 0

