Dear R users
I am trying to obtain p-values for (quasi)poisson lmer models, including
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 found the same
results as given.
Cheers
Beth

______________________________________________
[email protected] 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.

Reply via email to