Dear r-sig-mixed-models,



After discovering the glmm wikidot page on random effects testing (http://glmm.wikidot.com/random-effects-testing) I've been trying out Fabian Scheipl's RLRsim package to explore the importance of three random effects using lmer in the lme4 package. Before this I had been using maximum likelihood and anova to compare models, but have read that p value estimates may be overly conservative.



I'm interested in how the variance in isotopic signature (a proxy for wintering location and diet) varies as a function of BirdID (due to individual differences in wintering strategies), FeatherID (due to differences in moult patterns) and Sample (potential metabolic differences within feathers). There are 2 samples per feather (Rachis and Vane), and I have sampled between 5 and 12 feathers from each of 15 birds. So the design is nested but not fully balanced:

> xtabs(~Sample + Feather, data = data)

        Feather

Sample  C01 C02 P01 P02 P03 P04 P05 P06 P07 P08 P09 P10

Rachis  15  15   3  15   3   3  15   3   3   3  15   3

Vane    15  15   3  15   3   3  15   3   3   3  15   3



I'm testing the importance of each random effect by removing the term and using exactRLRT to calculate a LR and p value as follows:



fm0 = lmer(d13C ~ 1 + (1|Sample), data = data, REML = TRUE)

fm1 = lmer(d13C ~ 1 + (1|BirdID) + (1|FeatherID) + (1|Sample), data = data, REML = TRUE)

fm2 = lmer(d13C ~ 1 + (1|BirdID) + (1|FeatherID), data = data, REML = TRUE)

exactRLRT(fm0, fm1, fm2)



However, while exactRLRT will provide p values for BirdID and FeatherID, when testing for Sample I get the following error:



Error in chol.default(cov2cor(Vr)) :

  the leading minor of order 2 is not positive definite In addition: Warning message:

In cov2cor(Vr) : diag(.) had 0 or NA entries; non-finite result is doubtful



> traceback()

4: chol.default(cov2cor(Vr))

3: chol(cov2cor(Vr))

2: RLRTSim(X, Z, qrX = qrX, sqrt.Sigma = chol(cov2cor(Vr)), lambda0 = 0,

       seed = seed, nsim = nsim, log.grid.hi = log.grid.hi, log.grid.lo = log.grid.lo,

       gridlength = gridlength)

1: exactRLRT(fm0, fm1, fm2, nsim = 1e+05)



This suggests an error in the Vr matrix:

> extract.lmerDesign(fm0)$Vr

     [,1] [,2]

[1,]    0    0

[2,]    0    0



Has anyone faced a similar issue? Any help and advice would be greatly appreciated. Is it an issue with the way exactRLRT is running, or is it inappropriate for me to extract a p value using this method as Vr is 0 in this case? I have until now been examining differences between models using anova (and ML), which suggests there is little variation explained by Sample (X^2 = 2.8, p = 0.09).



> sessionInfo()

R version 2.15.2 (2012-10-26)

Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)



locale:

[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8



attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base



other attached packages:

[1] RLRsim_2.0-12   mgcv_1.7-22     lmtest_0.9-30   zoo_1.7-9       lme4_0.999999-0

[6] Matrix_1.0-9    lattice_0.20-10



loaded via a namespace (and not attached):

[1] grid_2.15.2   nlme_3.1-105  stats4_2.15.2 tools_2.15.2



Many thanks,



James



---------------------------------------------------

Dr James Grecian

Post-doctoral Research Associate

Institute of Biodiversity, Animal Health & Comparative Medicine

University of Glasgow

	[[alternative HTML version deleted]]

