if (FALSE) {
## Reproducing the MCMCglmm model
require(MCMCglmm)
data(charadriiformes)
## Setting up the model parameters:
## 1 - The formula (the first three PC axes)
model_formula <- cbind(PC1, PC2, PC3) ~ trait:clade-1
## 2 - The residual term
model_residuals <- ~us(trait):units
## 3 - The random terms
## (one per clade and one for the whole phylogeny)
model_randoms <- ~ us(at.level(clade,1):trait):animal
+ us(at.level(clade,2):trait):animal
+ us(at.level(clade,3):trait):animal
+ us(trait):animal
## Flat priors for the residuals and random terms
flat_priors <- list(
## The residuals priors
R = list(
R1 = list(V = diag(3), nu = 0.002)),
## The random priors (the phylogenetic terms)
G = list(
G1 = list(V = diag(3), nu = 0.002),
G2 = list(V = diag(3), nu = 0.002),
G3 = list(V = diag(3), nu = 0.002),
G4 = list(V = diag(3), nu = 0.002)))
## Run the model for 110000 iterations
## sampled every 100 with a burnin (discard)
## of the first 10000 iterations)
model <- MCMCglmm(formula = model_formula,
rcov = model_residual,
random = model_randoms,
family = rep("gaussian", 3),
prior = flat_priors,
nitt = 110000,
burnin = 10000,
thin = 100,
pedigree = charadriiformes$tree,
data = charadriiformes$data)
}
Run the code above in your browser using DataLab