help with getting pvalues in lme4 - problems with mcmcpvalue, aovlmer.fnc, and HPDinterval

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

HPDinterval(markov)

Error in UseMethod("HPDinterval") : no applicable method for
"HPDinterval"

So instead I ran:

pvals.fnc(modlmer)$fixed

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:

head(as.matrix(markov))

(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?

[1] 0

mcmcpvalue(as.matrix(markov)[, 1:3]) #Ag & Forest?

[1] 0

mcmcpvalue(as.matrix(markov)[, 3:4]) #Forest & Urban?

[1] 0.00062

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