# [R] lmer and estimation of p-values: error with mcmcpvalue()

Christoph Scherber Christoph.Scherber at agr.uni-goettingen.de
Mon Feb 12 13:58:23 CET 2007

```Dear all,

I am currently analyzing count data from a hierarchical design, and I?ve
tried to follow the suggestions for a correct estimation of p-values as
discusssed at R-Wiki
(http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests&s=lme%20and%20aov).

However, I have the problem that my model only consists of parameters
with just 1 d.f. (intercepts, slopes), so that the "mcmcpvalue" function
defined below obviously produces error messages.

How can I proceed in estimating the p-values, then?

I very much acknowledge any suggestions.

Best regards
Christoph.

##
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) }

m1<-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),quasipoisson)

Generalized linear mixed model fit using Laplace
Formula: number_pollinators ~ logpatch + loghab + landscape_diversity +
(1 | site)
Data: primula
AIC   BIC logLik deviance
84.83 93.75 -37.42    74.83
Random effects:
Groups   Name        Variance Std.Dev.
site     (Intercept) 0.036592 0.19129
Residual             1.426886 1.19452
number of obs: 44, groups: site, 15

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

Correlation of Fixed Effects:
(Intr) lgptch loghab
logpatch     0.091
loghab      -0.637 -0.121
lndscp_dvrs -0.483 -0.098 -0.348

markov1=mcmcsamp(m1,5000)
HPDinterval(markov1)

mcmcpvalue(as.matrix(markov1)[,1])

Error in colMeans(samp) : 'x' must be an array of at least two dimensions

```