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
______________________________________________
R-help@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.
