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