[R-sig-ME] variance estimates differ glmer vs lme
Vincenzo Ellis
vincenzoaellis at gmail.com
Wed Feb 15 17:47:54 CET 2017
Dear mixed effects modelers,
I'm interested in doing variance partitioning of parasite prevalence (a
proportion calculated as the number of infected host individuals divided by
the total number of host individuals sampled) among different host
taxonomic levels.
I've run the analysis in lme4 using the glmer function and I've also
arcsine square-root transformed the response variable (prevalence) and run
that in the nlme package using the lme function (with weights set to the
square root of sample size).
These two approaches give different answers. I assumed they would differ a
bit, but I expected them to be more similar than they are. The glmer
approach puts 66% of the variance in prevalence at the species within
genera level, while the lme approach puts 92% of the variance at the
species within genera level.
*Question:* Am I wrong to think that these two approaches should have given
more similar results? And is one of these approaches clearly correct or
incorrect?
Thanks for any insights!
Vincenzo
Reproducible example:
#### Compare hierarchical variance partitioning results between glmer() and
lme()
## install packages if necessary
install.packages(c("gsheet", "lme4", "nlme"))
## load data
library(gsheet)
dat <- gsheet2tbl('
https://docs.google.com/spreadsheets/d/1EhkrHG19pESBCTQ1mMa8od4lbyF-N484RchskCrPUqo/edit?usp=sharing
')
## glmer variance partitioning (takes about 30s to run on my computer)
library(lme4)
mod.glmm <- glmer(cbind(Infected, Sample_Size-Infected) ~ 1 +
(1|Family/Genus/Species),
data = dat, family = binomial)
## variance estimates
print(VarCorr(mod.glmm), comp=c("Variance"))
## proportional variance
mod.glmm.var <- c(VarCorr(mod.glmm)[[1]][1,1], VarCorr(mod.glmm)[[2]][1,1],
VarCorr(mod.glmm)[[3]][1,1])
round(mod.glmm.var/sum(mod.glmm.var), 3)
## lme variance partitioning
library(nlme)
prev <- with(dat, Infected/Sample_Size)
mod.lmm <- lme(asin(sqrt(prev)) ~ 1, random = ~ 1|Family/Genus,
data = dat, weights = ~ 1/sqrt(Sample_Size))
## variance estimates
VarCorr(mod.lmm)
## proportional variance
mod.lmm.var <- as.numeric(VarCorr(mod.lmm)[c(5,4,2),1])
round(mod.lmm.var/sum(mod.lmm.var), 3)
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list