"lme4"
(Bates et al 2011)
"rptR"
(Nakagawa and Schielzeth 2010)
"MCMCglmm"
(Hadfield 2010)
brain size and body size on primates ("primate_ind.txt"
, a tab separated text file), with more than one observation for the majority of species, individual-specific data
Import, transform and check individual specific data:
data.ind = read.table("primate_ind.txt", sep = "\t", header = TRUE) data.ind$lg_body = log(data.ind$body) data.ind$lg_brain = log(data.ind$brain) str(data.ind)
## 'data.frame': 894 obs. of 5 variables: ## $ species : Factor w/ 86 levels "Allenopithecus_nigroviridis",..: 1 2 3 3 3 3 3 3 3 3 ... ## $ body : num 2920 826 777 2310 4309 ... ## $ brain : num 57 45 42.5 52 54.2 ... ## $ lg_body : num 7.98 6.72 6.66 7.75 8.37 ... ## $ lg_brain: num 4.04 3.81 3.75 3.95 3.99 ...
May not compatible with the most recent version of R (need earlier releases to run the scripts below)
library(rptR, quietly = T) rpt(data.ind$lg_brain, data.ind$species, datatype = "Gaussian", method = "REML")
## ## Repeatability calculation using the LMM.REML method ## ## R = 0.985 ## SE = 0.002 ## CI = [0.979, 0.989] ## P = 0 [] ## 0.001 []
rpt(data.ind$lg_body, data.ind$species, datatype = "Gaussian", method = "ANOVA")
## ## Repeatability calculation using the ANOVA method ## ## R = 0.884 ## SE = 0.017 ## CI = [0.851, 0.917] ## P = 0 [] ## 0.001 []
To estimate repeatability directly from a mixed model
library(lme4, quietly = T) mod = lmer(lg_brain ~ 1 + (1 | species), data.ind) attr(lme4::VarCorr(mod)$species, "stddev")^2/(attr(lme4::VarCorr(mod)$species, "stddev")^2 + attr(lme4::VarCorr(mod), "sc")^2)
## (Intercept) ## 0.9854
To obtain P value based on randomization
null = c() for (i in 1:1000) { species_random = sample(data.ind$species, length(data.ind$species)) null[i] = anova(lmer(lg_brain ~ 1 + (1 | species_random), data.ind), mod, test = "Chisq")$Pr[2] } sum(null > 0.05)/length(null)
## [1] 0
To get 95% confidence intervals based on parametric bootstrap
rvalues <- numeric() for (i in 1:1000) { y = simulate(mod) mboot <- lmer(y[, 1] ~ 1 + (1 | species), data.ind) rvalues[i] = attr(lme4::VarCorr(mboot)$species, "stddev")^2/(attr(lme4::VarCorr(mboot)$species, "stddev")^2 + attr(lme4::VarCorr(mboot), "sc")^2) } quantile(rvalues, c(0.025, 0.975))
## 2.5% 97.5% ## 0.9798 0.9891
library(MCMCglmm, quietly = T) prior <- list(R = list(V = 1, nu = 0.002), G = list(G1 = list(V = 1, nu = 0.002))) mod.mcmc = MCMCglmm(lg_body ~ 1, random = ~species, data = data.ind, prior = prior, verbose = F) Rvalue = mod.mcmc$VCV[, 1]/(mod.mcmc$VCV[, 1] + mod.mcmc$VCV[, 2]) mean(Rvalue)
## [1] 0.8513
quantile(Rvalue, c(0.025, 0.975))
## 2.5% 97.5% ## 0.8086 0.8896