[R] lme4/lmer: P-Values from mcmc samples or chi2-tests?

[Manuel Morales]

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