I can?t seem to get sensible pvalues testing the effect of one factor
?landtype? in the following model using function mcmcpvalue; and it
seems from previous emails that it is no longer possible to use
aovlmer.fnc ( ) in lme4. And as of last night, I get new error
(?Error in UseMethod("HPDinterval") : no applicable method for
"HPDinterval")when trying to perform HPDinterval ().
I would appreciate any guidance on how to get pvalues to test one
factor using lmer function! In addition, I'd like to verify that the
best way to determine significant pairwise differences among factor
levels is to examine overlap of 95% CIs generated based on MCMC
(i.e., using HPDinterval or pvals.fnc).
I?m running the following model:
modlmer <- lmer(Chao1.Res ~ landtype + (1|patch), data = S.final,
REML=FALSE)
Where: landtype has 4 levels: Agriculture (as Intercept), Bauxite,
Forest, Urban.
Model output: Linear mixed model fit by maximum likelihood Formula:
Chao1.Res ~ landtype + (1 | patch) Data: S.final AIC BIC logLik
deviance REMLdev 1732 1753 -859.8 1720 1711 Random effects:
Groups Name Variance Std.Dev. patch (Intercept) 22.084
4.6993 Residual 36.615 6.0510 Number of obs: 253,
groups: patch, 99
Fixed effects: Estimate Std. Error t value (Intercept) 27.300
1.322 20.652 landtypeBauxite -6.357 1.772 -3.588
landtypeForest -1.049 1.727 -0.607 landtypeUrban -5.504
1.885 -2.920
Correlation of Fixed Effects: (Intr) lndtyB lndtyF landtypeBxt -0.746
landtypFrst -0.765 0.571 landtypUrbn -0.701 0.523 0.537
markov <- mcmcsamp(modlmer, 50000)
I get the below error (as of last night) when trying to run
HPDinterval
Error in UseMethod("HPDinterval") : no applicable method for
"HPDinterval"
So instead I ran:
Estimate MCMCmean HPD95lower HPD95upper pMCMC Pr(>|t|) (Intercept)
27.300 27.403 25.279 29.610 0.0001 0.0000 landtypeBauxite
-6.357 -6.519 -9.384 -3.603 0.0001 0.0012 landtypeForest
-1.049 -1.124 -3.864 1.790 0.4296 0.5829 landtypeUrban
-5.504 -5.610 -8.705 -2.677 0.0002 0.0077
Question: So should I look at HPD95lower and HPD95upper to determine
whether factor levels (Ag, Bauxite, Forest, or Urban) differ
pair-wise?
Trying to get pvalues using:
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) + }
I need pvalue to test landtype ? levels Agriculture (Intercept),
Bauxite, Forest, and Urban:
mcmcpvalue(as.matrix(markov)[, 1:4])
[1] 0
Pvalue of 0 does not make sense.
Verified that I wanted first 4 columns:
(Intercept) landtypeBauxite landtypeForest landtypeUrban ST1
sigma [1,] 27.30047 -6.357429 -1.048611 -5.503726
0.7766178 6.051015 [2,] 28.04026 -6.345041 -0.329192
-7.144239 0.5959559 5.379194 [3,] 28.83598 -8.727255
-4.117192 -5.349411 0.6346374 6.175085 [4,] 27.33180
-5.973496 -1.231745 -5.913463 0.3999811 6.015337 [5,]
27.68958 -5.833905 -2.323891 -5.954748 0.3312879
6.738105 [6,] 28.62220 -8.177245 -1.173439
-5.937315 0.3879350 7.291594
Just to check I tried difference columns and do not understand what
these pvalues represent:
mcmcpvalue(as.matrix(markov)[, 2:3]) #Bauxite & Forest?
mcmcpvalue(as.matrix(markov)[, 1:3]) #Ag & Forest?
mcmcpvalue(as.matrix(markov)[, 3:4]) #Forest & Urban?
mcmcpvalue(as.matrix(markov)[, 1:2]) #Ag & Bauxite?
[1] 0
Determined anti-conservative pvalue based on LRT:
mod0 <- lmer(Chao1.Res ~ 1 + (1|patch), data = S.final,
REML=FALSE) mod1 <- lmer(Chao1.Res ~ landtype + (1|patch), data =
S.final, REML=FALSE) anova(mod0,mod1)
Data: S.final Models: mod0: Chao1.Res ~ 1 + (1 | patch) mod1:
Chao1.Res ~ landtype + (1 | patch) Df AIC BIC logLik Chisq
Chi Df Pr(>Chisq) mod0 3 1743.43 1754.03 -868.71
mod1 6 1731.60 1752.80 -859.80 17.832 3 0.0004764 ***
Then tried to get pvalues based on aovlmer.fnc
mod1.aov = aovlmer.fnc (modlmer, mcmc = x$mcmc, which =
c("(Intercept)","landtypeBauxite", "landtypeForest",
"landtypeUrban" ))
Error in anova(object) : Calculated PWRSS for a LMM is negative
I have the following packages installed and loaded: require(lme4)
require(Matrix) require(lattice) require(SASmixed) require(coda)
require(languageR)
Christina Kennedy Behavior, Ecology, Evolution & Systematics
University of Maryland 3221 Biology-Psychology Building College Park,
Maryland 20742 Cell Phone: (202) 288-7483 Lab Phone: (301) 405-4512
Email: cmk6 at umd.edu _______________________________________________
R-sig-mixed-models at r-project.org mailing list