### set working directory ### ## require and source packages and functions ## require(geiger); require(phytools); require(fitContinuousMCMC); require(auteur); source("fitContinuous_paleoModels.R") ## tree ## phy <- read.tree("mammal.phy"); phy <- drop.tip(phy, "Sinoconodon") phy ## extant and extinct masses ## all.data <- read.csv("mammal_mass.csv", row.names=1) ##### Fit models to all mammaliaformes ##################### all.data <- treedata(phy, all.data)$data bm.fit <- fitContinuous_paleo(phy, all.data[,1], model = "BM",meserr = as.numeric(all.data[,2])) acdc.prior <- fitContinuousMCMC:::ACDC.prior(phy) acdc.fit <- fitContinuous_paleo(phy, all.data[,1], model = "EB", bounds = list(a = c(acdc.prior[[1]], acdc.prior[[2]])),meserr = as.numeric(all.data[,2])); ou.fit <- fitContinuous_paleo(phy, all.data[,1], model = "OU",meserr = as.numeric(all.data[,2])) white.fit <- fitContinuous_paleo(phy, all.data[,1], model = "white",meserr = as.numeric(all.data[,2])); trend.fit<- fitContinuous_paleo(phy, all.data[,1], model = "trend",meserr = as.numeric(all.data[,2])); shift.fit <- fitContinuous_paleo(phy, all.data[,1], model = "timeshift",meserr = as.numeric(all.data[,2]), shift.time = 66) release.fit <- fitContinuous_paleo(phy, all.data[,1], model = "release",meserr = as.numeric(all.data[,2]), bounds = list(alpha = c(0, 10)), shift.time = 66) releaseradiate.fit <- fitContinuous_paleo(phy, all.data[,1], model = "releaseradiate",meserr = as.numeric(all.data[,2]), bounds = list(alpha = c(0, 10)), shift.time = 66) aic.c <- c(bm.fit$Trait1$aicc, trend.fit$Trait1$aicc, acdc.fit$Trait1$aicc, ou.fit$Trait1$aicc, white.fit$Trait1$aicc, shift.fit$Trait1$aicc, release.fit$Trait1$aicc, releaseradiate.fit$Trait1$aicc); names(aic.c) <- c("Brownian\nMotion", "Directional \nTrend"," AC \n/DC", "Ornstein \nUhlenbeck", "White \nNoise", "K-Pg \nshift", "K-Pg \nrelease", "K-Pg release \n& radiate"); weights <- AkaikeWeights(aic.c); quartz(height = 5, width = 7); par(mar = c(4,4,1,1)) barplot(weights$Weights, names.arg = rownames(weights), ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, ylim = c(0,1)); mtext("Model", side = 1, line = 2.5, cex = 1.4) write.csv(weights, "alldata_aic_weights.csv"); # write results to csv file ###### what happens if we only use the extant taxa ####### ephy <- drop.fossil(phy, tol = 1); mean.extant <- all.data[match(ephy$tip.label, rownames(all.data)), 1] extant.sem <- all.data[match(ephy$tip.label, rownames(all.data)), 2] ex.bm.fit <- fitContinuous_paleo(ephy, mean.extant, model = "BM",meserr = as.numeric(extant.sem)) acdc.prior <- fitContinuousMCMC:::ACDC.prior(ephy) ex.acdc.fit <- fitContinuous_paleo(ephy, mean.extant, model = "EB", bounds = list(a = c(acdc.prior[[1]], acdc.prior[[2]])),meserr = as.numeric(extant.sem)); ex.ou.fit <- fitContinuous_paleo(ephy, mean.extant, model = "OU",meserr = as.numeric(extant.sem)) ex.white.fit <- fitContinuous_paleo(ephy, mean.extant, model = "white",meserr = as.numeric(extant.sem)); ex.trend.fit<- fitContinuous_paleo(ephy, mean.extant, model = "trend",meserr = as.numeric(extant.sem)); ex.shift.fit <- fitContinuous_paleo(ephy, mean.extant, model = "timeshift",meserr = as.numeric(extant.sem), shift.time = 66) ex.release.fit <- fitContinuous_paleo(ephy, mean.extant, model = "release",meserr = as.numeric(extant.sem), bounds = list(alpha = c(0, 10)), shift.time = 66) ex.releaseradiate.fit <- fitContinuous_paleo(ephy, mean.extant, model = "releaseradiate",meserr = as.numeric(extant.sem), bounds = list(alpha = c(0, 10)), shift.time = 66) aic.c <- c(ex.bm.fit$Trait1$aicc, ex.trend.fit$Trait1$aicc, ex.acdc.fit$Trait1$aicc, ex.ou.fit$Trait1$aicc, ex.white.fit$Trait1$aicc, ex.shift.fit$Trait1$aicc, ex.release.fit$Trait1$aicc, ex.releaseradiate.fit$Trait1$aicc); names(aic.c) <- c("Brownian\nMotion", "Directional \nTrend"," AC \n/DC", "Ornstein \nUhlenbeck", "White \nNoise", "K-Pg \nshift", "K-Pg \nrelease", "K-Pg release \n& radiate"); weights_extant <- AkaikeWeights(aic.c); quartz(height = 5, width = 7); par(mar = c(4,4,1,1)) barplot(weights_extant $Weights, names.arg = rownames(weights_extant), ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, ylim = c(0,1)); mtext("Model", side = 1, line = 2.5, cex = 1.4) write.csv(weights_extant, "extant_aic_weights.csv"); # write results to csv file ###### what happens if we only use the extant mammals and their fossil relatives? ####### keep <- node.leaves(phy, mrca.of.pair(phy, "Tachyglossidae", "Canidae")) ephy2 <- drop.tip(phy, phy$tip.label[-match(keep, phy$tip.label)]) mean.extantplusdesc <- all.data[match(ephy2$tip.label, rownames(all.data)), 1] extant.extantplusdesc <- all.data[match(ephy2$tip.label, rownames(all.data)),2] ex2.bm.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "BM",meserr = as.numeric(extant.extantplusdesc)) acdc.prior <- fitContinuousMCMC:::ACDC.prior(ephy2) ex2.acdc.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "EB", bounds = list(a = c(acdc.prior[[1]], acdc.prior[[2]])),meserr = as.numeric(extant.extantplusdesc)); ex2.ou.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc,model = "OU",meserr = as.numeric(extant.extantplusdesc)) ex2.white.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "white",meserr = as.numeric(extant.extantplusdesc)); ex2.trend.fit<- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "trend",meserr = as.numeric(extant.extantplusdesc)); ex2.shift.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "timeshift",meserr = as.numeric(extant.extantplusdesc), shift.time = 66) ex2.release.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "release",meserr = as.numeric(extant.extantplusdesc), bounds = list(alpha = c(0, 10)), shift.time = 66) ex2.releaseradiate.fit <- fitContinuous_paleo(ephy2, mean.extantplusdesc, model = "releaseradiate",meserr = as.numeric(extant.extantplusdesc), bounds = list(alpha = c(0, 10)), shift.time = 66) aic.c <- c(ex2.bm.fit$Trait1$aicc, ex2.trend.fit$Trait1$aicc, ex2.acdc.fit$Trait1$aicc, ex2.ou.fit$Trait1$aicc, ex2.white.fit$Trait1$aicc, ex2.shift.fit$Trait1$aicc, ex2.release.fit$Trait1$aicc, ex2.releaseradiate.fit$Trait1$aicc); names(aic.c) <- c("Brownian\nMotion", "Directional \nTrend"," AC \n/DC", "Ornstein \nUhlenbeck", "White \nNoise", "K-Pg \nshift", "K-Pg \nrelease", "K-Pg release \n& radiate"); weights_extant_plusdesc <- AkaikeWeights(aic.c); write.csv(weights_extant_plusdesc, "extant_plusdesc_aic_weights.csv"); quartz(height = 5, width = 7); par(mar = c(4,4,1,1)) barplot(weights_extant_plusdesc $Weights, names.arg = rownames(weights_extant_plusdesc), ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, ylim = c(0,1)); mtext("Model", side = 1, line = 2.5, cex = 1.4) # now barplot all weights for all models in one panel quartz(width = 5, height = 7) par(mfrow = c(3,1)); par(mar = c(4,4,1,1)); bar.order <- c("Brownian\nMotion", "Directional \nTrend", "Ornstein \nUhlenbeck", " AC \n/DC", "White \nNoise", "K-Pg \nshift", "K-Pg \nrelease", "K-Pg release \n& radiate"); barplot(weights_extant $Weights[match(bar.order, row.names(weights_extant))], names.arg = rownames(weights_extant)[match(bar.order, row.names(weights_extant))], ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, cex.axis = 0.8, ylim = c(0,1)); text(0.2, 0.92, "(a)", cex = 1.5) barplot(weights_extant_plusdesc $Weights[match(bar.order, row.names(weights_extant_plusdesc))], names.arg = rownames(weights_extant_plusdesc)[match(bar.order, row.names(weights_extant_plusdesc))], ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, cex.axis = 0.8, ylim = c(0,1)); text(0.2, 0.92, "(b)", cex = 1.5) barplot(weights$Weights[match(bar.order, row.names(weights))], names.arg = rownames(weights)[match(bar.order, row.names(weights))], ylab ="Akaike Weights", cex.names = 0.8, las = 1, cex.lab =1.4, cex.axis = 0.8, ylim = c(0,1)); text(0.2, 0.92, "(c)", cex = 1.5) mtext("model", side = 1, line = 2.5, cex = 1.1)