[R] lme4/lmer: P-Values from mcmc samples or chi2-tests?
Manuel Morales
Manuel.A.Morales at williams.edu
Wed Feb 14 21:54:41 CET 2007
On Tue, 2007-02-13 at 14:45 +0100, Christoph Scherber wrote:
> 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 *
<snip>
> 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:
<snip>
[MCMC]
> logpatch 0.004
> loghab 0.0598
> landscape_diversity 0.0074
<snip>
[single-term deletions]
> logpatch 0.007106
> loghab 0.1704
> landscape_diversity 0.01276
<snip>
> 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".
The P-values from the MCMC chain look suspiciously like 1/2 the others.
Are you effectively doing a 1-tailed test in your MCMC comparison?
--
Manuel A. Morales
http://mutualism.williams.edu
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20070214/8ea87faa/attachment.bin
More information about the R-help
mailing list