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

[Ben Bolker]

  Since your model is relatively simple (not a GLMM, not crossed
random effects) you are likely to have much better luck getting
p values with the nlme package at this point.

  lme(fixed=Chao1Res ~ landtype, random = ~1|patch, data=S.final,
   method="ML")

 should be about right -- then I believe anova() on the result
will give you an ANOVA table.

cmk6 at umd.edu wrote:
> 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 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc