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

[Henric Nilsson (Public)]

Den Må, 2007-02-12, 13:58 skrev 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])

Try `mcmcpvalue(as.matrix(markov1[,1]))'.


HTH,
Henric



>
> Error in colMeans(samp) : 'x' must be an array of at least two dimensions
>
> ______________________________________________
> R-help at 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.
>
>