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 [EMAIL PROTECTED] 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 > >______________________________________________ >[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. > > ______________________________________________ [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.
