"ape" (Paradis et al. 2004)
"geiger" (Harmon et al. 2008)
"nlme" (Pinheiro et al. 2014)
"picante" (Kembel et al. 2010)
"caper" (Orme et al. 2012)
"gtools" (Warnes et al. 2013)
"AICcmodavg" (Mazerolle 2013)
Species-specific trait data ("primate_data.txt"), a tab separated text file, species-specific data for body mass, gestation length, home range, longevity and social group size in primates
Sample of phylogenetic trees ("primate_tree.nex"), 100 phylogenetic hypotheses in nexus format, phylogeny of primates from Arnold et al. (2010)
# activate libraries library(ape) library(geiger) library(nlme) library(picante) library(caper) library(gtools) # import the first tree from the sample of phylogenies tree <- read.nexus("primate_tree.nex")[[1]] # import and transform data xdata <- read.table("primate_data.txt", sep = "\t", header = TRUE) xdata = data.frame(xdata, log.MaxLongevity_m = log(xdata$MaxLongevity_m), log.AdultBodyMass_g = log(xdata$AdultBodyMass_g), log.SocialGroupSize = log(xdata$SocialGroupSize), log.HomeRange = log(xdata$HomeRange_km2)) # prune phylogeny to the data rownames(xdata) = xdata$Binomial tree <- drop.tip(tree, setdiff(tree$tip.label, rownames(xdata)))
We first generate models with all possible combinations of predictors
pred.vars = c("GestationLen_d", "log.AdultBodyMass_g", "log.SocialGroupSize", "log.HomeRange") m.mat = permutations(n = 2, r = 4, v = c(F, T), repeats.allowed = T) models = apply(cbind(T, m.mat), 1, function(xrow) { paste(c("1", pred.vars)[xrow], collapse = "+") }) models = paste("log.MaxLongevity_m", models, sep = "~")
Then we define objects, in which we later store model outputs
# AIC of models all.aic = rep(NA, length(models)) # Estiamted lambdas all.lambda = rep(NA, length(models)) # Which predictors are estimated in the models beside the intercept m.mat = cbind(1, m.mat) colnames(m.mat) = c("(Intercept)", pred.vars) # number of parameters estimated in the models n.par = 2 + apply(m.mat, 1, sum)
Finally, we run all models and store the parameters of interests (AIC, lambda and model coefficients). This is the part where we will use the PGLS functions available in "caper". Basically, we use the same model in a loop that goes through all defined combinations of predictors.
primate <- comparative.data(phy = tree, data = xdata, names.col = Binomial, vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE) coefs=m.mat# define an object to store the coefficients for (k in 1:length(models)) { res = try(pgls(as.formula(models[k]), data = primate, lambda = "ML")) if (class(res) != "try-error") { all.aic[k] = -2 * logLik(res)[1] + 2 * n.par[k] all.lambda[k] = summary(res)$param[2] xx = coefficients(res) coefs[k, match(names(xx), colnames(m.mat))] = xx } }
If we want to know which model is associated with the lowest AIC:
min(all.aic)
## [1] 20.32
# this corresponds to the model models[which(all.aic == min(all.aic))]
## [1] "log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize"
# with the following parameters coefs[which(all.aic == min(all.aic)), ]
## (Intercept) GestationLen_d log.AdultBodyMass_g ## 4.61488 0.00000 0.11905 ## log.SocialGroupSize log.HomeRange ## 0.06502 0.00000
all.lambda[which(all.aic == min(all.aic))]
## [1] 0.742
If we want to do model averaging over the entire candidate model set:
# get summed Akaike weights per model aic_delta = all.aic - min(all.aic) w = exp(-0.5 * aic_delta)/sum(exp(-0.5 * aic_delta)) # get model averaged parameters (with estimate set to 0 inmodels in which # the term doesn't appear): unlist(lapply(1:ncol(coefs), function(x) { weighted.mean(x = coefs[, x], w = w, na.rm = T) }))
## [1] 4.7503043 0.0001867 0.1065480 0.0310901 0.0147193
# get model averaged parameters (considering only the models in which the # term appears): unlist(lapply(1:ncol(coefs), function(x) { x = coefs[, x] x[x == 0] = NA weighted.mean(x = x, w = w, na.rm = T) }))
## [1] 4.7503043 0.0006305 0.1147252 0.0595306 0.0319888
This is very similar to the above exercise, but here we will rely on the PGLS functions available in "nlme" in combination with "ape". Moreover, we can use the "AICcmodavg" functions for multimodel inference.
library(AICcmodavg) # define candidate models as above pred.vars = c("GestationLen_d", "log.AdultBodyMass_g", "log.SocialGroupSize", "log.HomeRange") m.mat = permutations(n = 2, r = 4, v = c(F, T), repeats.allowed = T) models = apply(cbind(T, m.mat), 1, function(xrow) { paste(c("1", pred.vars)[xrow], collapse = "+") }) models = paste("log.MaxLongevity_m", models, sep = "~") # Which predictors are estimated in the models beside the intercept m.mat = cbind(1, m.mat) colnames(m.mat) = c("(Intercept)", pred.vars) # run all models Cand.models = list() for (k in 1:length(models)) { Cand.models[[k]] = gls(as.formula(models[k]), data = xdata, method = "ML", correlation = corPagel(value = 0.5, tree, fixed = F)) } # get the AIC and AICc tables aictab(cand.set = Cand.models, modnames = models, sort = F, second.ord = F)
## ## Model selection based on AIC : ## ## K ## log.MaxLongevity_m~1 3 ## log.MaxLongevity_m~1+log.HomeRange 4 ## log.MaxLongevity_m~1+log.SocialGroupSize 4 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 5 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 4 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 5 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 5 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d 4 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 5 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 5 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 5 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 7 ## AIC ## log.MaxLongevity_m~1 31.70 ## log.MaxLongevity_m~1+log.HomeRange 25.75 ## log.MaxLongevity_m~1+log.SocialGroupSize 28.51 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 26.20 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 21.19 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 21.00 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 20.32 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 21.59 ## log.MaxLongevity_m~1+GestationLen_d 29.89 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 25.28 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 27.97 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 26.18 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 23.08 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 22.85 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 22.30 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 23.53 ## Delta_AIC ## log.MaxLongevity_m~1 11.38 ## log.MaxLongevity_m~1+log.HomeRange 5.42 ## log.MaxLongevity_m~1+log.SocialGroupSize 8.19 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 5.88 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 0.86 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 0.68 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 0.00 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 1.27 ## log.MaxLongevity_m~1+GestationLen_d 9.57 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 4.96 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 7.64 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 5.86 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 2.75 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 2.53 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 1.98 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 3.21 ## AICWt ## log.MaxLongevity_m~1 0.00 ## log.MaxLongevity_m~1+log.HomeRange 0.02 ## log.MaxLongevity_m~1+log.SocialGroupSize 0.00 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 0.01 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 0.15 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 0.17 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 0.23 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.12 ## log.MaxLongevity_m~1+GestationLen_d 0.00 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 0.02 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 0.01 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 0.01 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 0.06 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 0.07 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 0.09 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.05 ## LL ## log.MaxLongevity_m~1 -12.85 ## log.MaxLongevity_m~1+log.HomeRange -8.87 ## log.MaxLongevity_m~1+log.SocialGroupSize -10.25 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange -8.10 ## log.MaxLongevity_m~1+log.AdultBodyMass_g -6.59 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange -5.50 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize -5.16 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange -4.79 ## log.MaxLongevity_m~1+GestationLen_d -10.95 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange -7.64 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize -8.98 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange -7.09 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g -6.54 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange -5.43 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize -5.15 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange -4.77
aictab(cand.set = Cand.models, modnames = models, sort = T, second.ord = T)
## ## Model selection based on AICc : ## ## K ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 5 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 4 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 5 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 5 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 7 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 5 ## log.MaxLongevity_m~1+log.HomeRange 4 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 5 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 6 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 5 ## log.MaxLongevity_m~1+log.SocialGroupSize 4 ## log.MaxLongevity_m~1+GestationLen_d 4 ## log.MaxLongevity_m~1 3 ## AICc ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 21.17 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 21.74 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 21.85 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 22.79 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 23.50 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 23.92 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 24.05 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 25.16 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 26.13 ## log.MaxLongevity_m~1+log.HomeRange 26.30 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 27.04 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 27.38 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 28.81 ## log.MaxLongevity_m~1+log.SocialGroupSize 29.07 ## log.MaxLongevity_m~1+GestationLen_d 30.45 ## log.MaxLongevity_m~1 32.03 ## Delta_AICc ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 0.00 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 0.57 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 0.68 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 1.62 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 2.33 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 2.75 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 2.89 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 3.99 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 4.96 ## log.MaxLongevity_m~1+log.HomeRange 5.14 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 5.88 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 6.22 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 7.64 ## log.MaxLongevity_m~1+log.SocialGroupSize 7.90 ## log.MaxLongevity_m~1+GestationLen_d 9.28 ## log.MaxLongevity_m~1 10.86 ## AICcWt ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 0.24 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 0.18 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 0.17 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.11 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 0.07 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 0.06 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 0.06 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.03 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 0.02 ## log.MaxLongevity_m~1+log.HomeRange 0.02 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 0.01 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 0.01 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 0.01 ## log.MaxLongevity_m~1+log.SocialGroupSize 0.00 ## log.MaxLongevity_m~1+GestationLen_d 0.00 ## log.MaxLongevity_m~1 0.00 ## Cum.Wt ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize 0.24 ## log.MaxLongevity_m~1+log.AdultBodyMass_g 0.42 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange 0.59 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.70 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize 0.77 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g 0.84 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange 0.89 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange 0.92 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange 0.94 ## log.MaxLongevity_m~1+log.HomeRange 0.96 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange 0.98 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange 0.99 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize 0.99 ## log.MaxLongevity_m~1+log.SocialGroupSize 1.00 ## log.MaxLongevity_m~1+GestationLen_d 1.00 ## log.MaxLongevity_m~1 1.00 ## LL ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize -5.16 ## log.MaxLongevity_m~1+log.AdultBodyMass_g -6.59 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.HomeRange -5.50 ## log.MaxLongevity_m~1+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange -4.79 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize -5.15 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g -6.54 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.HomeRange -5.43 ## log.MaxLongevity_m~1+GestationLen_d+log.AdultBodyMass_g+log.SocialGroupSize+log.HomeRange -4.77 ## log.MaxLongevity_m~1+GestationLen_d+log.HomeRange -7.64 ## log.MaxLongevity_m~1+log.HomeRange -8.87 ## log.MaxLongevity_m~1+log.SocialGroupSize+log.HomeRange -8.10 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize+log.HomeRange -7.09 ## log.MaxLongevity_m~1+GestationLen_d+log.SocialGroupSize -8.98 ## log.MaxLongevity_m~1+log.SocialGroupSize -10.25 ## log.MaxLongevity_m~1+GestationLen_d -10.95 ## log.MaxLongevity_m~1 -12.85
Model averaging
# get model averaged parameters (considering only the # models in which the term appears): modav_pars = matrix(NA, length(colnames(m.mat)), 2) for (i in 1:length(colnames(m.mat))) { modav_par = modavg(parm = colnames(m.mat)[i], cand.set = Cand.models, modnames = models, second.ord = F) modav_pars[i, ] = cbind(modav_par$Mod.avg.beta, modav_par$Uncond.SE) } data.frame(row.names = colnames(m.mat), Estimate = modav_pars[, 1], SE = modav_pars[, 2])
## Estimate SE ## (Intercept) 4.7503044 0.375728 ## GestationLen_d 0.0006305 0.001641 ## log.AdultBodyMass_g 0.1147252 0.040836 ## log.SocialGroupSize 0.0595306 0.041366 ## log.HomeRange 0.0319887 0.025469
# get model averaged parameters (with estimate set to 0 in models in which # the term doesn't appear): for (i in 1:length(colnames(m.mat))) { modav_par = modavg.shrink(parm = colnames(m.mat)[i], cand.set = Cand.models, modnames = models, second.ord = F) modav_pars[i, ] = cbind(modav_par$Mod.avg.beta, modav_par$Uncond.SE) } data.frame(row.names = colnames(m.mat), Estimate = modav_pars[, 1], SE = modav_pars[, 2])
## Estimate SE ## (Intercept) 4.7503044 0.3757283 ## GestationLen_d 0.0001867 0.0009382 ## log.AdultBodyMass_g 0.1065479 0.0491932 ## log.SocialGroupSize 0.0310901 0.0421649 ## log.HomeRange 0.0147193 0.0235090