Dear R users,
I have now tried out several options of obtaining p-values for
(quasi)poisson lmer models, including Markov-chain Monte Carlo sampling
and single-term deletions with subsequent chi-square tests (although I
am aware that the latter may be problematic).
However, I encountered several problems that can be classified as
(1) the quasipoisson lmer model does not give p-values when summary() is
called (see below)
(2) dependence on the size of the mcmc sample
(3) lack of correspondence between different p-value estimates.
How can I proceed, left with these uncertainties in the estimations of
the p-values?
Below is the corresponding R code with some output so that you can see
it all for your own:
##
m1<-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),primula,poisson,method="ML")
m2<-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),primula,quasipoisson,method="ML")
summary(m1);summary(m2)
#m1: [...]
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.40302 0.57403 -0.7021 0.48262
logpatch 0.10915 0.04111 2.6552 0.00793 **
loghab 0.08750 0.06128 1.4279 0.15331
landscape_diversity 1.02338 0.40604 2.5204 0.01172 *
#m2: [...] #for the quasipoisson model, no p values are shown?!
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.4030 0.6857 -0.5877
logpatch 0.1091 0.0491 2.2228
loghab 0.0875 0.0732 1.1954
landscape_diversity 1.0234 0.4850 2.1099
##
anova(m1)
#here, no p-values appear (only in the current version of lme4)
Analysis of Variance Table
Df Sum Sq Mean Sq
logpatch 1 11.6246 11.6246
loghab 1 6.0585 6.0585
landscape_diversity 1 6.3024 6.3024
anova(m2)
Analysis of Variance Table
Df Sum Sq Mean Sq
logpatch 1 11.6244 11.6244
loghab 1 6.0581 6.0581
landscape_diversity 1 6.3023 6.3023
So I am left with the p-values estimated from the poisson model;
single-term deletion tests for each of the individual terms lead to
different p-values; I restrict this here to m2:
##
m2a=update(m2,~.-loghab)
anova(m2,m2a,test="Chi")
Data: primula
Models:
m2a: number_pollinators ~ logpatch + landscape_diversity + (1 | site)
m2: number_pollinators ~ logpatch + loghab + landscape_diversity + (1 |
site)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
m2a 4 84.713 91.850 -38.357
m2 5 84.834 93.755 -37.417 1.8793 1 0.1704
##
m2b=update(m2,~.-logpatch)
anova(m2,m2b,test="Chi")
Data: primula
Models:
m2b: number_pollinators ~ loghab + landscape_diversity + (1 | site)
m2: number_pollinators ~ logpatch + loghab + landscape_diversity +
m2b: (1 | site)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
m2b 4 90.080 97.217 -41.040
m2 5 84.834 93.755 -37.417 7.246 1 0.007106 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
m2c=update(m2,~.-landscape_diversity)
anova(m2,m2c,test="Chi")
Data: primula
Models:
m2c: number_pollinators ~ logpatch + loghab + (1 | site)
m2: number_pollinators ~ logpatch + loghab + landscape_diversity +
m2c: (1 | site)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
m2c 4 89.036 96.173 -40.518
m2 5 84.834 93.755 -37.417 6.2023 1 0.01276 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The p-values from mcmc are:
##
markov1=mcmcsamp(m2,5000)
HPDinterval(markov1)
lower upper
(Intercept) -1.394287660 0.6023229
logpatch 0.031154910 0.1906861
loghab 0.002961281 0.2165285
landscape_diversity 0.245623183 1.6442544
log(site.(In)) -41.156007604 -1.6993996
attr(,"Probability")
[1] 0.95
##
mcmcpvalue(as.matrix(markov1[,1])) #i.e. the p value for the intercept
[1] 0.3668
> mcmcpvalue(as.matrix(markov1[,2])) #i.e. the p-value for logpatch
[1] 0.004
> mcmcpvalue(as.matrix(markov1[,3])) #i.e. the p-value for loghab
[1] 0.0598
> mcmcpvalue(as.matrix(markov1[,4])) #i.e. the p-value for landscape.div
[1] 0.0074
If one runs the mcmcsamp function for, say, 50,000 runs, the p-values
are slightly different (not shown here).
##here are the p-values summarized in tabular form:
(Intercept) 0.3668
logpatch 0.004
loghab 0.0598
landscape_diversity 0.0074
##and from the single-term deletions:
(Intercept) N.A.
logpatch 0.007106
loghab 0.1704
landscape_diversity 0.01276
## Compare this with the p-values from the poisson model:
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.40302 0.57403 -0.7021 0.48262
logpatch 0.10915 0.04111 2.6552 0.00793 **
loghab 0.08750 0.06128 1.4279 0.15331
landscape_diversity 1.02338 0.40604 2.5204 0.01172 *
##
To summarize, at least for quasipoisson models, the p-values obtained
from mcmcpvalue() are quite different from those obtained using
single-term deletions followed by a chisquare test.
Especially in the case of "loghab", the difference is so huge that one
could tend to interpret one of the p-values as "marginally significant".
However, the mcmc chains look allright.
What would your suggestions be on how to proceed?
Thanks a lot for your help!
Best wishes,
Christoph.
##
I am using R 2.4.1 and the lme4 version I use is 0.9975-11 (Date:
2007-01-25)
______________________________________________
R-help@stat.math.ethz.ch 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.
