Sources

R packages

"lme4" (Bates et al 2011)

"rptR" (Nakagawa and Schielzeth 2010)

"MCMCglmm" (Hadfield 2010)

Data

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

Codes

Data input

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 ...

rptR

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 []

lme4

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

MCMCglmm

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

References