estrain at postoffice.utas.edu.au
2007-Mar-12 23:30 UTC
[R] Lmer Mcmc Summary and p values
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 aretotally 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
estrain at postoffice.utas.edu.au
2007-Mar-13 00:40 UTC
[R] Lmer Mcmc Summary and p values
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
estrain at postoffice.utas.edu.au
2007-Mar-13 01:00 UTC
[R] Lmer Mcmc Summary and p values
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 aretotally 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