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

[Christoph Scherber]

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
  Family: quasipoisson(log link)
    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