Christoph Scherber
2007-Feb-12 12:58 UTC
[R] lmer and estimation of p-values: error with mcmcpvalue()
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
Henric Nilsson (Public)
2007-Feb-12 13:26 UTC
[R] lmer and estimation of p-values: error with mcmcpvalue()
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. > >