Functions that aid in exploration of models from the packages MCMCglmm and lme4.
# Install this development version from GitHub:
# install.packages("devtools")
devtools::install_github("Euphrasiologist/VCVglmm")Please install them as follows:
install.packages("MCMCglmm"); install.packages("lme4")
library(MCMCglmm); library(lme4)
install.packages("data.table"); install.packages("tidyverse"); install.packages("ape"); install.packages("coda"); install.packages("lattice")
library(data.table); library(ggplot2); library(plyr); library(dplyr); library(reshape2); library(ape); library("coda"); library("lattice")Our toy data set will contain six parasites and six hosts, with 15 replicates of each parasite on each host. The response variable, Pterm, is simulated through normal deviations, with each parasite having different means, regardless of host (as we shall see…).
We can fit models to a very simple dataset about growth of parasites on host species.
head(dat.1); tail(dat.1)
#> Parasite Hosts Pterm
#> 1: P1 H1 3.641533
#> 2: P2 H1 12.280153
#> 3: P3 H1 11.834130
#> 4: P4 H1 19.699902
#> 5: P5 H1 23.683628
#> 6: P6 H1 27.577708
#> Parasite Hosts Pterm
#> 1: P1 H6 7.947823
#> 2: P2 H6 8.176867
#> 3: P3 H6 15.271369
#> 4: P4 H6 20.314061
#> 5: P5 H6 20.519166
#> 6: P6 H6 28.016953Firstly we will do a model in lme4.
library(lme4)
#> Loading required package: Matrix
mod.1 <- lmer(Pterm ~ 1 + (1 | Hosts) + (1 | Parasite), data = dat.1)
#> boundary (singular) fit: see ?isSingular
summary(mod.1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Pterm ~ 1 + (1 | Hosts) + (1 | Parasite)
#> Data: dat.1
#>
#> REML criterion at convergence: 2380.9
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.78200 -0.70653 0.02668 0.66659 2.71914
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> Hosts (Intercept) 0.000 0.000
#> Parasite (Intercept) 88.197 9.391
#> Residual 4.474 2.115
#> Number of obs: 540, groups: Hosts, 6; Parasite, 6
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 17.520 3.835 4.568
#> convergence code: 0
#> boundary (singular) fit: see ?isSingularSo this shows that almost all of the variance in our Pterm is being explained by Parasite, and the rest by the residual. But let’s check that further…
library(VCVglmm)
VarExpl.mer(mod.1)
#> Var % Var Explained
#> Hosts 0.00000 0
#> Parasite 88.19652 100And if we want that pesky p-value that lme4 doesn’t give you in the
summary(), don’t worry.
calc_pvals(mod.1)
#> Levels Estimate Std. Error t value pval_upperdf pval_lowerdf
#> 1: (Intercept) 17.51988 3.835063 4.568343 3.048217e-06 0.0003226639And if we quickly translate that to MCMCglmm() language…
library(MCMCglmm)
#> Loading required package: coda
#> Loading required package: ape
mod.2 <- MCMCglmm(Pterm ~ 1,
random = ~ Hosts + Parasite,
data = dat.1,
pr = TRUE,
verbose = FALSE)
summary(mod.2)
#>
#> Iterations = 3001:12991
#> Thinning interval = 10
#> Sample size = 1000
#>
#> DIC: 2349.483
#>
#> G-structure: ~Hosts
#>
#> post.mean l-95% CI u-95% CI eff.samp
#> Hosts 0.001058 3.097e-17 0.0006547 151.2
#>
#> ~Parasite
#>
#> post.mean l-95% CI u-95% CI eff.samp
#> Parasite 144.7 16.05 388.7 1000
#>
#> R-structure: ~units
#>
#> post.mean l-95% CI u-95% CI eff.samp
#> units 4.509 4.062 5.087 1000
#>
#> Location effects: Pterm ~ 1
#>
#> post.mean l-95% CI u-95% CI eff.samp pMCMC
#> (Intercept) 17.360 6.128 26.497 1000 0.012 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looks sensible to me, we might want to add a prior for real data
however. Below, using Solapply() we can extract the posterior modes of
each of the levels of the random effects to see which levels are having
the greatest impact.
Solapply(mod.2)
#> Variable Grouped_Value Group
#> 1: Hosts.H1 -6.041321e-04 Hosts
#> 2: Hosts.H2 -5.608347e-07 Hosts
#> 3: Hosts.H3 5.569833e-04 Hosts
#> 4: Hosts.H4 -3.068158e-04 Hosts
#> 5: Hosts.H5 3.354911e-04 Hosts
#> 6: Hosts.H6 -7.292986e-04 Hosts
#> 7: Parasite.P1 -1.172435e+01 Parasite
#> 8: Parasite.P2 -8.089398e+00 Parasite
#> 9: Parasite.P3 -1.141637e+00 Parasite
#> 10: Parasite.P4 3.126249e+00 Parasite
#> 11: Parasite.P5 8.828054e+00 Parasite
#> 12: Parasite.P6 1.366192e+01 ParasiteAgain we see that ~95%+ of the variability in our data is explained by Parasite, which was pretty much 100% in the lmer model.
MCMCRepnorm2(mod = mod.2)
#> % Var explained
#> Hosts -5.659259e-10
#> Parasite 9.798530e+01
#> units 2.014334e+00
#> $Hosts
#> lower upper
#> var1 2.80567e-18 0.0004912091
#> attr(,"Probability")
#> [1] 0.95
#>
#> $Parasite
#> lower upper
#> var1 89.11348 99.43389
#> attr(,"Probability")
#> [1] 0.95
#>
#> $units
#> lower upper
#> var1 0.5661093 10.88652
#> attr(,"Probability")
#> [1] 0.95glmm_explain() works on MCMCglmm, lmer and glmer models. It reports:
- convergence and fit checks;
- each fixed effect on the latent (link) scale and on the data scale. Rate ratios and changes in expected counts for Poisson models, odds ratios and changes in probability for binary models, all averaged over the random effects;
- the variance explained by each random term on the latent, link and data scales;
- marginal and conditional R².
glmm_explain(mod.2)glmm_varpart() generalises MCMCRepnorm2() and MCMCReppois() to Gaussian, Poisson, binary (logit/probit), ordinal and threshold models. It also handles random slopes and idh()/us() structures, with each term's variance averaged over the observations. Proportions are given on three scales:
- latent;
- link, which adds the distribution-specific variance;
- data, the observed scale, computed exactly or by Gauss-Hermite quadrature.
glmm_R2() gives Nakagawa & Schielzeth's marginal and conditional R².
These work on multi-response models (cbind(y1, y2) ~ trait - 1 + ...), giving results per trait with that trait's family. Ordinal and threshold models with more than two categories are also supported: their data scale uses the cutpoints and refers to the category number. Categorical responses with more than two categories get the latent scale per log-odds contrast, and the data scale per category: the variance in each category indicator explained by each term, computed jointly over the contrasts by simulation. glmm_explain() also reports how each covariate changes every category's probability. glmm_cor() extracts every correlation in the covariance structures, for example among-group vs residual correlations between traits, or intercept-slope correlations.
vp <- glmm_varpart(mod.2, fixed = TRUE)
vp
plot(vp)
glmm_R2(mod.2)
glmm_cor(mv.mod) # correlations, e.g. for a multi-response modelMCMCdiag(mod.2) # ESS, autocorrelation, Geweke, Heidelberger-Welch, boundary variances, correlated fixed effects...
MCMCtrace(mod.2) # ggplot2 traces with running means; pass a list of models to overlay chains
MCMCacfplot(mod.2) # autocorrelation by lag
MCMCcorplot(mod.2) # posterior correlations between fixed effects
MCMCpriorplot(mod.2, prior) # is the posterior any different from the prior?
MCMCgelman(chain1, chain2, chain3) # Gelman-Rubin R-hat across chainsmerdiag(mod.1) # convergence, singular fits, few levels, overdispersion, collinearity, separation
merdiagplot(mod.1) # residuals vs fitted, QQ plots of residuals and random effects, scale-locationPlease tell me where I’ve gone wrong and contribute if you can/want to!
- Common et al., 2020 Diversity in CRISPR‐based immunity protects susceptible genotypes by restricting phage spread and evolution, doi: https://doi.org/10.1111/jeb.13638