[R-sig-ME] help with getting pvalues in lme4 - problems with mcmcpvalue, aovlmer.fnc, and HPDinterval

[cmk6 at umd.edu]

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