##################################################################################################################################### ### CARNAges project part III – Body mass evolution. Script written by V Fischer 2017-2026 ### ### Related paper: Fischer, Solé et al. The rise of carnivoran mammals in Europe through the lens of body mass. Biology Letters. ### ##################################################################################################################################### #===================================# #============ Libraries ============# #===================================# #Get and analyze data source("http://www.graemetlloyd.com/pubdata/functions_2.r") #Age series detrending function by Dr. Graeme T. Lloyd source("Benson&al2017 Scripts S1/surface pBIC functions.R") source("Benson&al2017 Scripts S1/Trend shift functions.R") library(readxl) library(ape) library(phytools) library(vegan) library(dispRity) library(dplyr) library(castor) library(strap) library(paleotree) library(remotes) library(diptest) library(pbapply) #Graphs and figures library(ggplot2) library(ggpubr) library(gridExtra) library(geoscale) library(wesanderson) #======================================# #============ Data, Europe ============# #======================================# #============ Get and tidy bodymass data ============# #Bodymass, temporal, and climatic data bodymass<- as.data.frame(read_excel("SI2_Bodymass_data_EU.xlsx","BM")); bodymass$Clade_1 <- as.factor(bodymass$Clade_1); bodymass$Clade_2 <- as.factor(bodymass$Clade_2) num_ages <- as.data.frame(read_excel("SI2_Bodymass_data_EU.xlsx","Numerical_ages")) #numerical ages of MPs num_ages <- num_ages[order(num_ages$Midpoint,decreasing=TRUE),] westerhold <- as.data.frame(read_excel("SI4_Westerhold&al.xlsx")) westerhold <- westerhold[westerhold$Age>min(num_ages$Midpoint) & westerhold$Age0){ temp <- do.call(rbind, replicate((range[i]), bodymass[i,], simplify=FALSE)) temp[,"FAD_MP"] <- (temp[1,"FAD_MP"]+1):(temp[1,"FAD_MP"]+range[i]) supp <- rbind(supp, temp) } } supp <- supp[-1,] unique_bodymass <- bodymass$Log10BM; names(unique_bodymass) <- bodymass$Species #stores previous bodymass df for phylogenetic uses (phenograms, surface) bodymass <- rbind(bodymass,supp) bodymass <- bodymass[,-which(colnames(bodymass)=="LAD_MP")] #removes LAD_MP column as every multibin taxon is now duplicated bodymass <- bodymass[order(bodymass[,"FAD_MP"]),] bodymass <- cbind(bodymass,num_ages[bodymass$FAD_MP,"Midpoint"]) colnames(bodymass)[c(4,6)] <- c("Bin","Age") ; bodymass <- bodymass[,c(1:4,6,5)]; rownames(bodymass) <- 1:nrow(bodymass) #Optional: removes bodymasses <100g which are probably an issue of equations being unfit for small taxa (Eoproviverra, Parvagula) bodymass <- bodymass[10^(bodymass$Log10BM)>0.1,] unique_bodymass <- unique_bodymass[10^unique_bodymass>0.1] #============ Build phylogenetic data ============# # Most of the dataset originates from Faurby et al 2024 Proc B # Unsampled species in that dataset for which the genus is present were added by simply replacing one of the species names by the species we have bodymass data for # This assumes genus monophyly # Clades and lineages not sampled in Faurby et al 2024 were added manually through grafting: we use the tree structure of Spaulding et al 2009 PLoS One (Mesonychia as sister to Carnivoramorpha + Hyaenodonta + Oxyaenodonta) Faurby <- read.nexus("SI5_Faurby&al_nexus.nex") mesonychian_tree <- read.tree(text="((Pachyaena_gigantea:1,(Pachyaena_cf_ossifraga:1,Pachyaena_cf_gracilis:1):1):1,(Dissacus_europaeus:1,(Dissacus_rougierae:1,(Dissacus_raslanloubatieri:1,(Dissacus_blayaci:1,(Dissacus_progressus:1,Dissacus_filholi:1):1):1):1):1):1):1;") trees <- Faurby #Start replacing taxa by those we have data for, within the same genus or recent generic reattribution (Dissacus, see Solé et al. 2023 J of Mammalian Evolution) #step 1, old-new name list name_map <- c("Amphictis_aginensis"="Amphictis_ambigua", "Amphictis_antiqua"="Amphictis_milloquensis", "Apterodon_altidens"="Apterodon_rauenbergensis", "Bathygale_lemanensis"="Bathygale_julieni", "Cephalogale_shareri"="Cephalogale_brevirostris", "Cephalogale_ginesticus"="Cephalogale_minor", "Dissacus_blayaci"="Hyaenodictis_blayaci", "Dissacus_filholi"="Hyaenodictis_filholi", "Dissacus_progressus"="Hyaenodictis_progressus", "Dissacus_raslanloubatieri"="Hyaenodictis_raslanloubatieri", "Dissacus_rougierae"="Hyaenodictis_rougierae", "Hyaenodon_brevirostris"="Hyaenodon_exiguus", "Hyaenodon_chunkhtensis"="Hyaenodon_heberti", "Hyaenodon_crucians"="Hyaenodon_minor", "Hyaenodon_eminus"="Hyaenodon_requieni", "Miacis_deutschi"="Miacis?_macintyri", "Mustelictis_piveteaui"="Mustelictis_cf_major", "Pachycynodon_tedfordi"="Pachycynodon_crassirostris", "Palaeoprionodon_lamandini"="Palaeoprionodon_mutabilis", "Paramiacis_teilharti"="Paramiacis_teilhardi", "Plesictis_humilidens"="Plesictis_bransatensis", "Plesictis_julieni"="Plesictis_stenogalinus", "Pseudocyonopsis_antiquus"="Pseudocyonopsis_ambiguus", "Stenogale_aurelianensis"="Stenogale_bransatensis", "Stenogale_brevidens"="Stenogale_crassa", "Stenogale_julieni"="Stenogale_gracilis", "Stenoplesictis_muhoronii"="Stenoplesictis_minor", "Viverravus_acutus"="Viverravus_lawsoni") #step 2, change in all trees of the multiphylo object trees <- lapply(trees, function(tree) { tree$tip.label <- ifelse(tree$tip.label %in% names(name_map),name_map[tree$tip.label],tree$tip.label) tree }) class(trees) <- "multiPhylo" #Now let's add new tips supertrees <- pblapply(trees, function(tree) { #step 1, create a small function to place new tip at a position that is 1/2 the length of the branch where they have to go safe_position <- function(which_node) { edge <- which(tree$edge[,2]==which_node) #if (identical(edge, integer(0)) || length(edge)==0 || tree$edge.length[edge]<=0 || is.null(tree$edge.length[edge]) || is.na(tree$edge.length[edge])) return(0) if (length(edge) == 0) return(0) tree$edge.length[edge]*0.5 } #step 2, make each of the new tip within a function to be used with lapply to affect all the trees in the multiphylo object which_node <- get_mrca_of_set(tree, c("Plesictis_genetoides","Plesictis_dieupentalensis","Plesictis_stenogalinus","Plesictis_bransatensis")) tree <- bind.tip(tree, "Bavarictis_gaimersheimensis", edge.length=1,where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Cynelos_lemanensis","Cynelos_rugosidens","Cynelos_crassidens","Ysengrinia_tolosana","Crassidia_intermedia","Goupilicitis_minor","Pseudocyonopsis_ambiguus")) tree <- bind.tip(tree, "Brachycyon_gaudryi", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Brachycyon_gaudryi") tree <- bind.tip(tree, "Brachycyon_palaeolycos", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Brachycyon_gaudryi","Cynelos_rugosidens","Brachycyon_palaeolycos")) tree <- bind.tip(tree, "Harpagocyon_inusitatus", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Harpagocyon_inusitatus","Brachycyon_gaudryi","Cynelos_rugosidens","Brachycyon_palaeolycos")) tree <- bind.tip(tree, "Symplectocyon_praecursor", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Leonhardtina_godinoti","Leonhardtina_gracilis")) tree <- bind.tip(tree, "Cartierodon_egerkingenis", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Cynodictis_lacustris") tree <- bind.tip(tree, "Cynodictis_compressidens", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Cynodictis_compressidens") tree <- bind.tip(tree, "Cynodictis_palmidens", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Haplocyon_elegans") tree <- bind.tip(tree, "Haplocyon_dombrowskii", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Leonhardtina_gracilis") tree <- bind.tip(tree, "Leonhardtina_meridianum", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Matthodon_tritens") tree <- bind.tip(tree, "Matthodon_peignei", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Oxyaenoides_lindgreni") tree <- bind.tip(tree, "Oxyaenoides_aumelasiensis", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Hoplophoneus_bidentatus","Hoplophoneus_cerebralis", "Hoplophoneus_mentalis","Hoplophoneus_occidentalis", "Hoplophoneus_primaevus","Hoplophoneus_sicarius", "Hoplophoneus_villebramarensis")) tree <- bind.tip(tree, "Eusmilus_bidentatus", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Eusmilus_bidentatus") tree <- bind.tip(tree, "Eusmilus_villebramarensis", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Pachycynodon_crassirostris") tree <- bind.tip(tree, "Pachycynodon_boriei", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Pachycynodon_boriei") tree <- bind.tip(tree, "Pachycynodon_cf_dubius", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- which(tree$tip.label=="Pseudocyonopsis_ambiguus") tree <- bind.tip(tree, "Pseudocyonopsis_landesquei", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Cynodictis_lacustris","Cynodictis_palmidens","Cynodictis_compressidens")) tree <- bind.tip(tree, "Storchictis_miacinus", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Haplocyon_crassidens","Haplocyon_elegans","Haplocyon_crucians")) tree <- bind.tip(tree, "Sarcocyon_ferox", edge.length=1, where=which_node, position=safe_position(which_node)) which_node <- get_mrca_of_set(tree, c("Stenoplesictis_cayluxi","Stenogale_crassa")) tree <- bind.tip(tree, "Viretictis_gaillardi", edge.length=1, where=which_node, position=safe_position(which_node)) tree <- drop.tip(tree, tree$tip.label[!tree$tip.label %in% bodymass$Species]) tree$root.edge <- 1 which_node <- get_mrca_of_set(tree,c("Palaeonictis_gigantea","Nimravus_intermedius")) mesonychian_tree$root.edge <- 1 tree <- bind.tree(tree, mesonychian_tree, where=which_node, position=safe_position(which_node)) tree$root.edge <- 1 tree }) class(supertrees) <- "multiPhylo" #get temporal data for a posteriori timescaling FAD_LAD <- as.data.frame(read_excel("SI2_Bodymass_data_EU.xlsx","BM"))[,c(3:5)] colnames(FAD_LAD)[2] <- "MP" FAD_LAD <- FAD_LAD %>% left_join(num_ages[,1:2], by = "MP") colnames(FAD_LAD) <- c("Species","FAD","MP","Start") FAD_LAD <- FAD_LAD %>% left_join(num_ages[,c(1,3)], by = "MP") rownames(FAD_LAD)<-FAD_LAD$Species; FAD_LAD<-FAD_LAD[,4:5];colnames(FAD_LAD) <- c("FAD","LAD") FAD_LAD <- FAD_LAD[rownames(FAD_LAD) %in% bodymass$Species,] supertrees_mbl <- pblapply(supertrees, function(tree){ tree <- multi2di(tree,random=TRUE) timePaleoPhy(tree,timeData=FAD_LAD,type= "mbl",vartime=1,dateTreatment="firstLast") }) class(supertrees_mbl) <- "multiPhylo" pdf("Figures/FigS1_Timescaled_composite_tree_mbl.pdf",width=15, height=15) geoscalePhylo(ladderize(supertrees_mbl[[1]],right=TRUE),FAD_LAD,cex.ts=0.5,cex.tip=0.6) dev.off() #For Extended Hedman, it is the timing of the First appearance that matters, therefore: FAD <- FAD_LAD$FAD; names(FAD) <- rownames(FAD_LAD) cl <- makeCluster(detectCores() - 1) clusterExport(cl, c("FAD", "supertrees")) clusterEvalQ(cl, { library(ape) library(strap) }) supertrees_hedman <- pblapply(supertrees, function(tree){ tree <- multi2di(tree,random=TRUE) DatePhyloHedman(tree=multi2di(tree), tip.ages=FAD, outgroup.ages=c(63.4,57.5), t0 = 66, resolution = 1000, conservative = TRUE) #outgroups: Eoconodon 63.4 Ma minimum, Goniacodon 57.5 Ma min (Halliday et al 2017 Biol Rev), t0= K/Pg boundary }, cl=cl) stopCluster(cl) class(supertrees_hedman) <- "multiPhylo" pdf("Figures/FigS2_Timescaled_composite_tree_hedman.pdf",width=15, height=15) geoscalePhylo(ladderize(supertrees_hedman[[1]]$tree,right=TRUE),FAD_LAD,cex.ts=0.5,cex.tip=0.6) dev.off() #============================================# #============ Phenograms, Europe ============# #============================================# #compute ancestral values for nodes using maximum likelihood anc_bodymass <- fastAnc(supertrees_mbl[[1]], unique_bodymass) #Extract data from the phenogram function to add temporal information on the plots pheno_data <- phenogram(supertrees_mbl[[1]],c(unique_bodymass,anc_bodymass),spread.labels=FALSE) #!! if spread.labels is set on TRUE (default) then the x data is that of the label position, not the end point of the branch!! root <- rowMeans(FAD_LAD)[["Stenogale_bransatensis"]]+pheno_data["Stenogale_bransatensis",1] #the general value from which temporal should be substracted to be placed in the correct position on the graph; this is a weird behaviour of the phenogram function pdf("Figures/FigS4_Phenogram_mbl.pdf", width=200/24.5, height=200/24.5) phenogram(supertrees_mbl[[1]],c(unique_bodymass,anc_bodymass),ylim=c(min(unique_bodymass)-0.18,max(unique_bodymass)),ylab="log10(bodymass)",fsize = 0.1,ftype="i",spread.labels=FALSE) rect(xleft=root-strips$xmin,ybottom=min(unique_bodymass)-0.05,xright=root-strips$xmax,ytop=max(unique_bodymass)+0.05,col="#cccccc50",border="NA") #strips rect(xleft=root-c(66,56,33.9),ybottom=min(unique_bodymass)-0.25,xright=root-c(56,33.9,23.04),ytop=min(unique_bodymass)-0.05,col=c("#fda75f","#fdb46c","#fec07a"),border=NA) text(x=c(root-c(mean(c(66,56)),mean(c(56,33.9)),mean(c(33.9,24.04)))),y=rep(min(unique_bodymass)-0.15,3),labels=c("Paleocene","Eocene","Oligocene"),cex=0.7) points(pheno_data,pch=20,col=cols[bodymass[match(rownames(pheno_data),bodymass$Species),"Clade_1"]]) #greenish "#23CC8F" = non-euelasmos – # orange "#ff952c" = euelasmos – pink "#FF0D68" = aristonectines dev.off() #compute ancestral values for nodes using maximum likelihood anc_bodymass <- fastAnc(supertrees_hedman[[1]]$tree, unique_bodymass) #Extract data from the phenogram function to add temporal information on the plots pheno_data <- phenogram(supertrees_hedman[[1]]$tree,c(unique_bodymass,anc_bodymass),spread.labels=FALSE) #!! if spread.labels is set on TRUE (default) then the x data is that of the label position, not the end point of the branch!! root <- rowMeans(FAD_LAD)[["Stenogale_bransatensis"]]+pheno_data["Stenogale_bransatensis",1] #the general value from which temporal should be substracted to be placed in the correct position on the graph; this is a weird behaviour of the phenogram function pdf("Figures/FigS5_Phenogram_hedman.pdf", width=200/24.5, height=200/24.5) phenogram(supertrees_hedman[[1]]$tree,c(unique_bodymass,anc_bodymass),ylim=c(min(unique_bodymass)-0.18,max(unique_bodymass)),ylab="log10(bodymass)",fsize = 0.1,ftype="i",spread.labels=FALSE) rect(xleft=root-strips$xmin,ybottom=min(unique_bodymass)-0.05,xright=root-strips$xmax,ytop=max(unique_bodymass)+0.05,col="#cccccc50",border="NA") #strips rect(xleft=root-c(66,56,33.9),ybottom=min(unique_bodymass)-0.25,xright=root-c(56,33.9,23.04),ytop=min(unique_bodymass)-0.05,col=c("#fda75f","#fdb46c","#fec07a"),border=NA) text(x=c(root-c(mean(c(66,56)),mean(c(56,33.9)),mean(c(33.9,24.04)))),y=rep(min(unique_bodymass)-0.15,3),labels=c("Paleocene","Eocene","Oligocene"),cex=0.7) points(pheno_data,pch=20,col=cols[bodymass[match(rownames(pheno_data),bodymass$Species),"Clade_1"]]) #greenish "#23CC8F" = non-euelasmos – # orange "#ff952c" = euelasmos – pink "#FF0D68" = aristonectines dev.off() #============================================# #============ Fit models, Europe ============# #============================================# sample_size <- 100 tree_sample <- sample(1:length(supertrees),sample_size,replace=FALSE) model_fit <- mclapply(1:sample_size, function(i){ tree_mbl <- supertrees_mbl[[tree_sample[i]]] trait_mbl <- unique_bodymass[tree_mbl$tip.label] #taxon order needs to be similar as in the tree, so it is best to recompute at each iteration tree_hedman <- supertrees_hedman[[tree_sample[i]]]$tree trait_hedman <- unique_bodymass[tree_hedman$tip.label] #taxon order needs to be similar as in the tree, so it is best to recompute at each iteration #Uniform models bm_mbl<- fitContinuous(tree_mbl, trait_mbl, SE=NA, model="BM") trend_mbl<- fitContinuous(tree_mbl, trait_mbl, SE=NA, model="drift") # uniform trend bm_hedman<- fitContinuous(tree_hedman, trait_hedman, SE=NA, model="BM") trend_hedman<- fitContinuous(tree_hedman, trait_hedman, SE=NA, model="drift") # uniform trend #TrendShift models with 1,2,3 node shifts ts1_mbl<- findBestShiftP(trait_mbl, tree_mbl, se=trend_mbl$opt$SE, vs.init=trend_mbl$opt$sigsq, ncore=1) #SE is the one estimated by geiger above ts2_mbl<- findBestShiftP(trait_mbl, tree_mbl, se=trend_mbl$opt$SE, bn=ts1_mbl$bn, vs.init=trend_mbl$opt$sigsq, ncore=1) # use previous solution for bn; SE is the one estimated by geiger above ts3_mbl<- findBestShiftP(trait_mbl, tree_mbl, se=trend_mbl$opt$SE, bn=ts2_mbl$bn, vs.init=trend_mbl$opt$sigsq, ncore=1) # use previous solution for bn; SE is the one estimated by geiger above ts1_hedman<- findBestShiftP(trait_hedman, tree_hedman, se=trend_hedman$opt$SE, vs.init=trend_hedman$opt$sigsq, ncore=1) #SE is the one estimated by geiger above ts2_hedman<- findBestShiftP(trait_hedman, tree_hedman, se=trend_hedman$opt$SE, bn=ts1_hedman$bn, vs.init=trend_hedman$opt$sigsq, ncore=1) # use previous solution for bn; SE is the one estimated by geiger above ts3_hedman<- findBestShiftP(trait_hedman, tree_hedman, se=trend_hedman$opt$SE, bn=ts2_hedman$bn, vs.init=trend_hedman$opt$sigsq, ncore=1) # use previous solution for bn; SE is the one estimated by geiger above #Surface surface_pbic_mbl <- runSurface_pBIC(nameNodes(tree_mbl),as.data.frame(trait_mbl),plotaic=FALSE) full_mbl_summary <- surfaceSummary(surface_pbic_mbl$fwd,surface_pbic_mbl$bwd) surface_pbic_hedman <- runSurface_pBIC(nameNodes(tree_hedman),as.data.frame(trait_hedman),plotaic=FALSE) full_hedman_summary <- surfaceSummary(surface_pbic_hedman$fwd,surface_pbic_hedman$bwd) #AICc vectors list(aic_mbl = c(bm_mbl$opt$aicc, trend_mbl$opt$aicc, ts1_mbl$AICc, ts2_mbl$AICc, ts3_mbl$AICc, min(full_mbl_summary$aics)), aic_hedman = c(bm_hedman$opt$aicc, trend_hedman$opt$aicc, ts1_hedman$AICc, ts2_hedman$AICc, ts3_hedman$AICc, min(full_hedman_summary$aics)), models_mbl = list(bm=bm_mbl, trend=trend_mbl, ts1=ts1_mbl, ts2=ts2_mbl, ts3=ts3_mbl, surface=surface_pbic_mbl), models_hedman = list(bm=bm_hedman, trend=trend_hedman, ts1=ts1_hedman, ts2=ts2_hedman, ts3=ts3_hedman, surface=surface_pbic_hedman)) }, mc.cores=detectCores()-1) #proportions of models with lowest AICc: aic_vector_mbl <- as.data.frame(t(sapply(model_fit, function(x) x$aic_mbl))); colnames(aic_vector_mbl) <- c("BM", "Trend", "nodeShift1", "nodeShift2", "nodeShift3", "surface") table(apply(aic_vector_mbl,1,function(row) names(which.min(row))))/nrow(aic_vector_mbl) aic_vector_hedman <- as.data.frame(t(sapply(model_fit, function(x) x$aic_hedman))); colnames(aic_vector_hedman) <- c("BM", "Trend", "nodeShift1", "nodeShift2", "nodeShift3", "surface") table(apply(aic_vector_hedman,1,function(row) names(which.min(row))))/nrow(aic_vector_hedman) #hedman write.csv(cbind(Timescaling=rep(c("mbl","hedman"),each=sample_size),(rbind(aic_vector_mbl,aic_vector_hedman))), "Stats/TableS1O_phylomodels.csv") #figures of best models for each timescaling method min_aic_mbl <- which(aic_vector_mbl==min(aic_vector_mbl), arr.ind = TRUE) #minimal AICc of all models, all trees: 3 node shift, tree sample 100 best_model_tree_mbl <- supertrees_mbl[[min_aic_mbl[1]]] best_model_trait_mbl <- unique_bodymass[best_model_tree_mbl$tip.label] best_model_mbl <- model_fit[[min_aic_mbl[1]]]$models_mbl[[min_aic_mbl[2]]] min_aic_hedman <- which(aic_vector_hedman==min(aic_vector_hedman), arr.ind = TRUE) #minimal AICc of all models, all trees: : 3 node shift, tree sample 65 best_model_tree_hedman <- supertrees_hedman[[min_aic_hedman[1]]]$tree best_model_trait_hedman <- unique_bodymass[best_model_tree_hedman$tip.label] best_model_hedman <- model_fit[[min_aic_hedman[1]]]$models_hedman[[min_aic_hedman[2]]] pdf("Figures/FigS6_best_model_mbl.pdf",height=400/24.5,width=400/24.5) showTrendShift(best_model_trait_mbl, best_model_tree_mbl, w=best_model_mbl) # best model overall mbl on our iteration: tree 100, ts3 dev.off() pdf("Figures/FigS7_best_model_hedman.pdf",height=400/24.5,width=400/24.5) showTrendShift(best_model_trait_hedman, best_model_tree_hedman, w=best_model_hedman) # best model overall hedman on our iteration: tree 65, ts3 dev.off() #====================================================================# #============ Build matrix of bodymasses per bin, Europe ============# #====================================================================# # Generate empty data frames bm_bin_EU <- bm_bin_EU_Hyaenodonta <- bm_bin_EU_early_Carnivoramorpha <- bm_bin_EU_Carnivora <- bm_bin_EU_All_Carnivoramorpha <- bm_bin_EU_Oxyaenodonta <- bm_bin_EU_Mesonychia <- data.frame(Base=num_ages$Start, End=num_ages$End,Age=num_ages$Midpoint,sampled_div=0, mean=NA, lowerQuant=NA, upperQuant=NA, min=NA, max=NA, variance=NA,row.names=num_ages$MP) for (i in 1:nrow(num_ages)){ temp_bin <- bodymass[bodymass$Bin==i,] if(nrow(temp_bin)>0){ bm_bin_EU[i,4:10] <- temp_bin%>%summarize(n(),mean(Log10BM),quantile(Log10BM,probs=0.25),quantile(Log10BM,probs=0.75),min(Log10BM),max(Log10BM),var(Log10BM)) temp_df <- temp_bin%>%group_by(Clade_1)%>%summarize(n(),mean(Log10BM),quantile(Log10BM,probs=0.25),quantile(Log10BM,probs=0.75),min(Log10BM),max(Log10BM),var(Log10BM)) bm_bin_EU_All_Carnivoramorpha[i,4:10] <- temp_bin[temp_bin$Clade_1=="Carnivoramorpha"| temp_bin$Clade_1=="Carnivora",]%>%summarize(n(),mean(Log10BM),quantile(Log10BM,probs=0.25),quantile(Log10BM,probs=0.75),min(Log10BM),max(Log10BM),var(Log10BM)) if("Hyaenodonta" %in% temp_df$Clade_1){ bm_bin_EU_Hyaenodonta[i,4:10] <- as.data.frame(temp_df %>% filter(Clade_1=="Hyaenodonta"))[,2:8]} if("Carnivoramorpha" %in% temp_df$Clade_1){ bm_bin_EU_early_Carnivoramorpha[i,4:10] <- as.data.frame(temp_df %>% filter(Clade_1=="Carnivoramorpha"))[,2:8]} if("Carnivora" %in% temp_df$Clade_1){ bm_bin_EU_Carnivora[i,4:10] <- as.data.frame(temp_df %>% filter(Clade_1=="Carnivora"))[,2:8]} if("Oxyaenodonta" %in% temp_df$Clade_1){ bm_bin_EU_Oxyaenodonta[i,4:10] <- as.data.frame(temp_df %>% filter(Clade_1=="Oxyaenodonta"))[,2:8]} if("Mesonychia" %in% temp_df$Clade_1){ bm_bin_EU_Mesonychia[i,4:10] <- as.data.frame(temp_df %>% filter(Clade_1=="Mesonychia"))[,2:8]} } } write.csv(bm_bin_EU,"Stats/TableS1A_bodymass_over_time_all_EU.csv") write.csv(bm_bin_EU_Hyaenodonta,"Stats/TableS1E_bodymass_over_time_Hyaenodonta.csv") write.csv(bm_bin_EU_early_Carnivoramorpha,"Stats/TableS1B_bodymass_over_time_early_carnivoramorpha.csv") write.csv(bm_bin_EU_Carnivora,"Stats/TableS1C_bodymass_over_time_Carnivora.csv") write.csv(bm_bin_EU_All_Carnivoramorpha,"Stats/TableS1D_bodymass_over_time_All_Carnivoramorpha.csv") write.csv(bm_bin_EU_Oxyaenodonta,"Stats/TableS1F_bodymass_over_time_Oxyaenodonta.csv") write.csv(bm_bin_EU_Mesonychia,"Stats/TableS1G_bodymass_over_time_Mesonychia.csv") #============ Body mass variance per bin and sampling intensity – correlation test ============# correl <- matrix(nrow=2,ncol=2) rownames(correl) <- c("Raw_variance_VS_sampling","Gen_diff_variance_VS_sampling") colnames(correl) <- c("Spearman_rho","p") #Using raw values a <- cor.test(bm_bin_EU[-9,"variance"],bm_bin_EU[-9,"sampled_div"],method="spearman") #line 9 is MP15 for which we have no data correl[1,1] <- a$estimate correl[1,2] <- a$p.value #After generalised differencing of the 'time' series gd.var<- gen.diff(bm_bin_EU[-9,"variance"],bm_bin_EU[-9,"Age"]) gd.sample<- gen.diff(bm_bin_EU[-9,"sampled_div"],bm_bin_EU[-9,"Age"]) gd <- as.data.frame(cbind(num_ages[c(1:8,10:(dim(num_ages)[1]-1)),"Midpoint"],gd.var,gd.sample)) colnames(gd) <- c("Age","gd.var","gd.sample") b <- cor.test(gd$gd.var,gd$gd.sample,method="spearman") correl[2,1] <- b$estimate correl[2,2] <- b$p.value write.csv(correl,"Stats/TableS1N_correlation_bias.csv") ggplot(gd,aes(x=gd.sample,y=gd.var,colour=Age))+ geom_point(size=3)+ scale_colour_gradientn(colours=colorRampPalette(c("yellow","darkorange"),space="rgb")(dim(gd)[1]))+ geom_smooth(method='lm', se=FALSE,colour="black")+ xlab("Gen. diff. sampled diversity")+ ylab("Gen. diff. variance")+ annotate("text",x=-5,y=0.025,label=paste("italic(R) ^ 2 == ",round(cor(gd$gd.var,gd$gd.sample)^2,digits=5)),parse=TRUE)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) ggsave(file="Figures/FigS3_VarianceVSdiv.pdf",device="pdf",units="mm",width=200,height=180) #============ Body mass through Time figures ============# #All taxa, with regression lines fitted for each major clade regress_clade1_plot <- ggplot(data=na.omit(bodymass),aes(x=Age,y=Log10BM,color=Clade_1,fill=Clade_1))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ annotate("text", x=c(42,32,32,55,55), y=c(-0.8,2.1,-0.1,2.3,0.75), label=c("'Early carnivoramorphans'","Carnivora","Hyaenodonta","Mesonychia","Oxyenidae"),color=c(cols[2],cols[1],cols[3:5]),hjust=0,size=3.5)+ geom_smooth(show.legend=FALSE)+ geom_point(size=1,alpha=0.5,show.legend=FALSE)+ scale_color_manual(values=cols)+ scale_fill_manual(values=cols)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text",x=num_ages$Midpoint,y=-1,label=num_ages$MP,size=2)+ #adds MP numbers annotate("text", x=41.7, y=2.5, label="MECO",color="black",size=3.5)+ annotate("text", x=31, y=2.5, label="EOT",color="black",size=3.5)+ scale_x_reverse(limits=c(23,57.4))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Bodymass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),legend.position="none",legend.direction="horizontal",legend.title=element_blank()) ggsave(filename="Figures/FigS8_regression_Age_clades.pdf",regress_clade1_plot,device="pdf",units="mm",width=180,height=100) #for info: per subfamilies ggplot(bodymass,aes(x=Age,y=Log10BM,colour=Clade_2))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_point(alpha=0.5)+ geom_smooth(method="lm",na.rm=TRUE)+ scale_x_reverse(limits=c(23,57.4))+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=1)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=1)+ annotate("text",x=c(58,num_ages$Midpoint),y=-1.7,label=c("MPs:",num_ages$MP))+ #adds MP numbers annotate("text", x=41.7, y=2.5, label="MECO",color="black")+ annotate("text", x=30, y=2.5, label="EOT",color="black")+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Bodymass (log10, kg)")+ guides(color=guide_legend(title="Clade (lower rank)"))+ theme(panel.background = element_rect(fill="white"),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),legend.position="bottom",legend.direction="horizontal",legend.title=element_blank()) ### Plot of mean bodymass, variance, and sampling #Mean bodymass, all taxa mean_bodymass_all <- ggplot(bm_bin_EU,aes(x=Age,y=mean))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data=bm_bin_EU,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data=bm_bin_EU,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ annotate("text",y=1.5,x=46,label="Max bodymass",alpha=0.5,size=1.7)+ geom_line(data=bm_bin_EU,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ annotate("text",y=-0.9,x=47,label="Min bodymass",alpha=0.5,size=1.7)+ geom_ribbon(data=bm_bin_EU,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[6],alpha=0.5)+ geom_line(data=bm_bin_EU,aes(x=Age,y=mean),color=cols[6],size=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass, carnivoraforms mean_bodymass_carni <- ggplot()+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ geom_line(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ geom_ribbon(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[1],alpha=0.5)+ geom_line(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,y=mean),color=cols[1],linewidth=1)+ geom_point(data=bodymass[bodymass$Age>40 & bodymass$Age<=45.6 & bodymass$Clade_1=="Carnivoramorpha",5:6],aes(x=Age,y=Log10BM),size=1,color=cols[2])+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass hyaenodontids mean_bodymass_hyaeno <- ggplot(bm_bin_EU_Hyaenodonta,aes(x=Age,y=mean))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data=bm_bin_EU_Hyaenodonta,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data=bm_bin_EU_Hyaenodonta,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ geom_line(data=bm_bin_EU_Hyaenodonta,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ geom_ribbon(data=bm_bin_EU_Hyaenodonta,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[3],alpha=0.5)+ geom_line(data=bm_bin_EU_Hyaenodonta,aes(x=Age,y=mean),color=cols[3],linewidth=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #============ Bodymass histograms tests & figures ============# ### Tests for multimodality using Hartigans' Dip Test Dip_Test_table <- data.frame(D=rep(NA,15),pvalue=rep(NA,15), row.names=c("All","All_preMECO","All_postMECO","All_preGC","All_postGC","Carnivoramorpha","Carnivoramorpha_preMECO","Carnivoramorpha_postMECO","Carnivoramorpha_preGC","Carnivoramorpha_postGC","Hyaenodonta","Hyaenodonta_preMECO","Hyaenodonta_postMECO","Hyaenodonta_preGC","Hyaenodonta_postGC")) #All taxa, before and after MECO then GC test <- dip.test(bodymass[,"Log10BM"]) # * Dip_Test_table[1,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[which(bodymass$Bin<10),"Log10BM"]) # * Dip_Test_table[2,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[which(bodymass$Bin>10),"Log10BM"]) # *** Dip_Test_table[3,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[which(bodymass$Bin<15),"Log10BM"]) # NS Dip_Test_table[4,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[which(bodymass$Bin>15),"Log10BM"]) # ** Dip_Test_table[5,] <- c(test$statistic,test$p.value) #Carnivoramorphans, before and after MECO then GC test <- dip.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha"),"Log10BM"]) # ** Dip_Test_table[6,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin>10,"Log10BM"]) # *** Dip_Test_table[7,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin<10,"Log10BM"]) # NS Dip_Test_table[8,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin>15,"Log10BM"]) # ** Dip_Test_table[9,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin<15,"Log10BM"]) # NS Dip_Test_table[10,] <- c(test$statistic,test$p.value) #Hyaenodonts, before and after MECO then GC test <- dip.test(bodymass[bodymass$Clade_1=="Hyaenodonta","Log10BM"]) # *** Dip_Test_table[11,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin>10,"Log10BM"]) # ** Dip_Test_table[12,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin<10,"Log10BM"]) # NS Dip_Test_table[13,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin>15,"Log10BM"]) # *** Dip_Test_table[14,] <- c(test$statistic,test$p.value) test <- dip.test(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin<15,"Log10BM"]) # NS Dip_Test_table[15,] <- c(test$statistic,test$p.value) write.csv(Dip_Test_table,"Stats/TableS1P_test_multimodality.csv") ### Figures temperatures_plot <- ggplot(data=westerhold,aes(x=Age,y=d18O_10pts_smooth,color=d18O_10pts_smooth))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_point(size=1,alpha=0.5)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ scale_x_reverse(limits=c(23,57.4))+ scale_y_reverse()+ scale_color_gradient(low="red", high="blue")+ ylab("LOESS smoothed d18O")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),legend.position="none") hist_MECO_all <- ggplot()+ geom_histogram(data=bodymass[bodymass$Bin>10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7, binwidth=0.5, fill=cols[6],linewidth=0.5,boundary=1) + geom_step(data=bodymass[bodymass$Bin<10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid", binwidth=0.5, color=cols[6], linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-MECO",color=cols[6],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-MECO",color=cols[6],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) hist_GC_all <- ggplot()+ geom_histogram(data=bodymass[bodymass$Bin>15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7, binwidth=0.5, fill=cols[6],linewidth=0.5,boundary=1) + geom_step(data=bodymass[bodymass$Bin<15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid", binwidth=0.5, color=cols[6], linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-GC",color=cols[6],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-GC",color=cols[6],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) hist_MECO_hyaeno <- ggplot()+ geom_histogram(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin>10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7,binwidth=0.5,fill=cols[3],linewidth=0.6,boundary=1) + geom_step(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin<10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid",binwidth=0.5,color=cols[3],linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-MECO",color=cols[3],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-MECO",color=cols[3],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) hist_GC_hyaeno <- ggplot()+ geom_histogram(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin>15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7,binwidth=0.5,fill=cols[3],linewidth=0.6,boundary=1) + geom_step(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Bin<15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid",binwidth=0.5,color=cols[3],linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-GC",color=cols[3],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-GC",color=cols[3],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) hist_MECO_carni <- ggplot()+ geom_histogram(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin>10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7,binwidth=0.5,fill=cols[1],linewidth=0.5,boundary=1) + geom_step(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin<10,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid",binwidth=0.5,color=cols[2],linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-MECO",color=cols[2],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-MECO",color=cols[1],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) hist_GC_carni <- ggplot()+ geom_histogram(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin>15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",alpha=0.7,binwidth=0.5,fill=cols[1],linewidth=0.5,boundary=1) + geom_step(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Bin<15,], mapping=aes(x=Log10BM, y=after_stat(ncount)),stat="bin",direction="mid",binwidth=0.5,color=cols[2],linewidth=0.7,boundary=1) + annotate("text", x=-0.8, y=0.9, label="Pre-GC",color=cols[2],size=1.7)+ annotate("text", x=2.5, y=0.9, label="Post-GC",color=cols[1],size=1.7)+ scale_x_continuous(limits=c(-1, 2.5))+ labs(x="",y="")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Assemble into Fig1 save_fig <- grid.arrange(temperatures_plot+labs(tags="A"), mean_bodymass_all+labs(tags="B"), hist_MECO_all+labs(tags="C"), hist_GC_all+labs(tags="D"), mean_bodymass_carni+labs(tags="E"), hist_MECO_carni+labs(tags="F"), hist_GC_carni+labs(tags="G"), mean_bodymass_hyaeno+labs(tags="H"), hist_MECO_hyaeno+labs(tags="I"), hist_GC_hyaeno+labs(tags="J"), heights=rep(1,8),widths=c(2,1),layout_matrix=cbind(c(1,1,2,2,5,5,8,8),c(NA,NA,3,4,6,7,9,10))) ggsave(file="Figures/Fig1_temp.pdf",save_fig,device="pdf",units="mm",width=180,height=250) #Variance & sampled div bm_df_to_plot <- rbind(bm_bin_EU,bm_bin_EU_All_Carnivoramorpha,bm_bin_EU_Hyaenodonta) bm_df_to_plot$clade <- factor(c(rep("All",nrow(bm_bin_EU)),rep("Carnivoramorpha",nrow(bm_bin_EU_All_Carnivoramorpha)),rep("Hyaenodonta",nrow(bm_bin_EU_Hyaenodonta)))) sampled_div <- ggplot(bm_df_to_plot,aes(x=Age,y=sampled_div,color=clade,width=0.35))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='white', alpha=0.7,inherit.aes=FALSE)+ geom_line(size=1,show.legend=FALSE)+ scale_color_manual(values=cols)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",size=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",size=0.5)+ scale_x_reverse(limits=c(23,57.4))+ labs(title="Sampled diversity, EU",x="Age (Myr) (data binned per MP 6 to 30)",y="Number of taxa")+ theme(axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) variance <- ggplot(bm_df_to_plot,aes(x=Age,y=variance,color=clade,width=0.35))+ geom_rect(data=strips, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='white', alpha=0.7,inherit.aes=FALSE)+ geom_line(size=1,show.legend=FALSE)+ scale_color_manual(values=cols)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",size=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",size=0.5)+ scale_x_reverse(limits=c(23.1,57.4))+ labs(title="Body mass variance, EU",x="Age (Myr) (data binned per MP 6 to 30)",y="Variance")+ theme(axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #==========================================================================# #============ Tests for non-phylogenetic regime shifts, Europe ============# #==========================================================================# ###Regime shifts using serialised statistic: two-sample Kolmogorov–Smirnov test #Create blank matrix KS <- matrix(nrow=dim(bm_bin_EU)[1],ncol=8); rownames(KS) <- rownames(bm_bin_EU); colnames(KS) <- c("All_D","All_p-value","All_Carnivoramorpha_D","All_Carnivoramorpha_p-value","Non_carnivoran_Carnivoramorpha_D","Non_carnivoran_Carnivoramorpha_p-value","Hyaenodonta_D","Hyaenodonta_p-value") #All taxa for (i in 1:(nrow(num_ages))-1){ if (sum(bodymass$Age>num_ages[i,"Midpoint"])>1){ test <- ks.test(bodymass[bodymass$Age>num_ages[i,"Midpoint"],"Log10BM"],bodymass[bodymass$Agenum_ages[i,"Midpoint"])>1){ test <- ks.test(bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Age>num_ages[i,"Midpoint"],"Log10BM"],bodymass[(bodymass$Clade_1=="Carnivora" | bodymass$Clade_1=="Carnivoramorpha") & bodymass$Agenum_ages[i,"Midpoint"])>1){ if (sum(bodymass$Clade_1=="Carnivoramorpha" & bodymass$Age1){ test <- ks.test(bodymass[bodymass$Clade_1=="Carnivoramorpha" & bodymass$Age>num_ages[i,"Midpoint"],"Log10BM"],bodymass[bodymass$Clade_1=="Carnivoramorpha" & bodymass$Agenum_ages[i,"Midpoint"])>1){ test <- ks.test(bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Age>num_ages[i,"Midpoint"],"Log10BM"],bodymass[bodymass$Clade_1=="Hyaenodonta" & bodymass$Age 23.03,c(3,4,10,11)] Juhn_df <- na.omit(Juhn_df) Juhn_df$mid_age <- rowMeans(Juhn_df[,1:2]) #Get temporal data sub-Nalma ages from Gunnell 2008 nalma <- as.data.frame(read_excel("SI3_Juhn&al2024data.xlsx",sheet="NALMA",col_names=TRUE)) rownames(nalma) <- nalma[,1]; nalma <- nalma[,-1] strips_USA <- data.frame(xmin=nalma$End[c(TRUE,FALSE)],xmax=nalma$Base[c(TRUE,FALSE)],ymin=0,ymax=35,key="xmin") strips_USA <- strips_USA[order(strips_USA$xmin),] #Build body mass over time datasets bm_bin_USA <- bm_bin_USA_Hyaenodonta <- bm_bin_USA_early_Carnivoramorpha <- bm_bin_USA_Carnivora <- bm_bin_USA_All_Carnivoramorpha <- bm_bin_USA_Oxyaenodonta <- data.frame(Base=nalma$Base, End=nalma$End,Age=rowMeans(nalma[,c("Base","End")]),sampled_div=0, mean=NA, lowerQuant=NA, upperQuant=NA, min=NA, max=NA, variance=NA,row.names=rownames(nalma)) for (i in 1:nrow(nalma)){ temp <- Juhn_df[Juhn_df$Max_Age >= nalma[i,"End"] & Juhn_df$Min_Age <= nalma[i,"Base"],c("Clade","Log10BM")] bm_bin_USA[i,4:10] <- temp%>%summarize(n(),mean(Log10BM),quantile(Log10BM,probs=0.25),quantile(Log10BM,probs=0.75),min(Log10BM),max(Log10BM),var(Log10BM)) test <- temp%>%group_by(Clade)%>%summarize(n(),mean(Log10BM),quantile(Log10BM,probs=0.25),quantile(Log10BM,probs=0.75),min(Log10BM),max(Log10BM),var(Log10BM)) if("Hyaenodonta" %in% test$Clade){ bm_bin_USA_Hyaenodonta[i,4:10] <- as.data.frame(test %>% filter(Clade=="Hyaenodonta"))[,2:8]} if("Carnivoramorpha" %in% test$Clade){ bm_bin_USA_early_Carnivoramorpha[i,4:10] <- as.data.frame(test %>% filter(Clade=="Carnivoramorpha"))[,2:8]} if("Carnivora" %in% test$Clade){ bm_bin_USA_Carnivora[i,4:10] <- as.data.frame(test %>% filter(Clade=="Carnivora"))[,2:8]} if("Oxyaenodonta" %in% test$Clade){ bm_bin_USA_Oxyaenodonta[i,4:10] <- as.data.frame(test %>% filter(Clade=="Oxyaenodonta"))[,2:8]} bm_bin_USA_All_Carnivoramorpha[i,4:10] <- as.data.frame(test %>% filter(Clade=="Carnivoramorpha"| Clade=="Carnivora"))[,2:8] } rownames(bm_bin_USA) <- rownames(bm_bin_USA_Hyaenodonta) <- rownames(bm_bin_USA_early_Carnivoramorpha) <- rownames(bm_bin_USA_Carnivora) <- rownames(bm_bin_USA_All_Carnivoramorpha) <- rownames(bm_bin_USA_Oxyaenodonta)<- rownames(nalma) write.csv(bm_bin_USA,"Stats/TableS1H_USA_bodymass_over_time_all.csv") write.csv(bm_bin_USA_Hyaenodonta,"Stats/TableS1L_USA_bodymass_over_time_Hyaenodonta.csv") write.csv(bm_bin_USA_early_Carnivoramorpha,"Stats/TableS1I_USA_bodymass_over_time_early_carnivoramorpha.csv") write.csv(bm_bin_USA_Carnivora,"Stats/TableS1J_USA_bodymass_over_time_Carnivora.csv") write.csv(bm_bin_USA_All_Carnivoramorpha,"Stats/TableS1K_USA_bodymass_over_time_All_Carnivoramorpha.csv") write.csv(bm_bin_USA_Oxyaenodonta,"Stats/TableS1M_USA_bodymass_over_time_Oxyaenodonta.csv") #Ggplots #Mean bodymass, all taxa mean_bodymass_all_USA <- ggplot(bm_bin_USA,aes(x=Age,y=mean))+ geom_rect(data=strips_USA, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data=bm_bin_USA,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data=bm_bin_USA,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ annotate("text",y=1.8,x=46,label="Max bodymass",alpha=0.5,size=1.7)+ geom_line(data=bm_bin_USA,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ annotate("text",y=-0.6,x=47,label="Min bodymass",alpha=0.5,size=1.7)+ geom_ribbon(data=bm_bin_USA,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[6],alpha=0.5)+ annotate("text",y=1.3,x=46,label="Q3",alpha=0.5,color=cols[6],size=1.7)+ annotate("text",y=-0.1,x=46,label="Q1",alpha=0.5,color=cols[6],size=1.7)+ geom_line(data=bm_bin_USA,aes(x=Age,y=mean),color=cols[6],size=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per sub-NALMA)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass, carnivoraforms mean_bodymass_carni_USA <- ggplot()+ geom_rect(data=strips_USA, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data= bm_bin_USA_All_Carnivoramorpha,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data= bm_bin_USA_All_Carnivoramorpha,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ geom_line(data= bm_bin_USA_All_Carnivoramorpha,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ geom_ribbon(data=bm_bin_USA_All_Carnivoramorpha,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[1],alpha=0.5)+ annotate("text",y=0.3,x=46,label="Q3",alpha=0.5,color=cols[1],size=1.7)+ annotate("text",y=-0.2,x=46,label="Q1",alpha=0.5,color=cols[1],size=1.7)+ geom_line(data= bm_bin_USA_All_Carnivoramorpha,aes(x=Age,y=mean),color=cols[1],linewidth=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass hyaenodontids mean_bodymass_hyaeno_USA <- ggplot(bm_bin_USA_Hyaenodonta,aes(x=Age,y=mean))+ geom_rect(data=strips_USA, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='grey95', alpha=0.7,inherit.aes=FALSE)+ geom_ribbon(data= bm_bin_USA_Hyaenodonta,aes(x=Age,ymin=min,ymax=max),fill="grey80",alpha=0.5)+ geom_line(data= bm_bin_USA_Hyaenodonta,aes(x=Age,y=max),color="black",linewidth=0.5,alpha=0.5)+ geom_line(data= bm_bin_USA_Hyaenodonta,aes(x=Age,y=min),color="black",linewidth=0.5,alpha=0.5)+ geom_ribbon(data=bm_bin_USA_Hyaenodonta,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[3],alpha=0.5)+ annotate("text",y=1.6,x=46,label="Q3",alpha=0.5,color=cols[3],size=1.7)+ annotate("text",y=0.5,x=46,label="Q1",alpha=0.5,color=cols[3],size=1.7)+ geom_line(data= bm_bin_USA_Hyaenodonta,aes(x=Age,y=mean),color=cols[3],linewidth=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=41.7, y=2.25, label="MECO",color="black")+ annotate("text", x=30, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) ggarrange(mean_bodymass_all_USA, mean_bodymass_carni_USA, mean_bodymass_hyaeno_USA,ncol=1) ggsave(file="Figures/FigS9_BM_USA.pdf",device="pdf",units="mm",width=200,height=180) #Mean bodymass, all taxa EU vs USA mean_bodymass_all_EUvsUSA <- ggplot()+ geom_ribbon(data=bm_bin_USA,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill="grey",alpha=0.3)+ geom_line(data=bm_bin_USA,aes(x=Age,y=mean),color="grey",size=1)+ geom_ribbon(data=bm_bin_EU,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[6],alpha=0.3)+ geom_line(data=bm_bin_EU,aes(x=Age,y=mean),color=cols[6],size=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=42.5, y=2.25, label="MECO",color="black")+ annotate("text", x=32, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per sub-NALMA)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass, carnivoraforms EU vs USA mean_bodymass_carni_EUvsUSA <- ggplot()+ geom_ribbon(data=bm_bin_USA_All_Carnivoramorpha,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill="grey",alpha=0.3)+ geom_line(data= bm_bin_USA_All_Carnivoramorpha,aes(x=Age,y=mean),color="grey",linewidth=1)+ geom_ribbon(data=bm_bin_EU_All_Carnivoramorpha,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[1],alpha=0.3)+ geom_line(data= bm_bin_EU_All_Carnivoramorpha,aes(x=Age,y=mean),color=cols[1],linewidth=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=42.5, y=2.25, label="MECO",color="black")+ annotate("text", x=32, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Mean bodymass, hyaenodonts EU vs USA mean_bodymass_hyaeno_EUvsUSA <- ggplot()+ geom_ribbon(data=bm_bin_USA_Hyaenodonta,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill="grey",alpha=0.3)+ geom_line(data= bm_bin_USA_Hyaenodonta,aes(x=Age,y=mean),color="grey",linewidth=1)+ geom_ribbon(data=bm_bin_EU_Hyaenodonta,aes(x=Age,ymin=lowerQuant,ymax=upperQuant),fill=cols[3],alpha=0.3)+ geom_line(data= bm_bin_EU_Hyaenodonta,aes(x=Age,y=mean),color=cols[3],linewidth=1)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",linewidth=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",linewidth=0.5)+ annotate("text", x=42.5, y=2.25, label="MECO",color="black")+ annotate("text", x=32, y=2.25, label="EOT",color="black")+ scale_x_reverse(limits=c(23,57.4))+ scale_y_continuous(limits=c(-1, 2.3))+ xlab("Age (Myr) (data binned per MP 6 to 30)")+ ylab("Body mass (log10, kg)")+ theme(panel.background = element_rect(fill="white"),axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) ggarrange(mean_bodymass_all_EUvsUSA, mean_bodymass_carni_EUvsUSA, mean_bodymass_hyaeno_EUvsUSA,ncol=1) ggsave(file="Figures/FigS10_MeanBM_EUvsUSA.pdf",device="pdf",units="mm",width=200,height=400) #Variance & sampled div bm_df_to_plot_USA <- rbind(bm_bin_USA,bm_bin_USA_All_Carnivoramorpha,bm_bin_USA_Hyaenodonta) bm_df_to_plot_USA$clade <- factor(c(rep("All",nrow(bm_bin_USA)),rep("Carnivoramorpha",nrow(bm_bin_USA_All_Carnivoramorpha)),rep("Hyaenodonta",nrow(bm_bin_USA_Hyaenodonta)))) sampled_div_USA <- ggplot(bm_df_to_plot_USA,aes(x=Age,y=sampled_div,color=clade,width=0.35))+ geom_rect(data=strips_USA, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='white', alpha=0.7,inherit.aes=FALSE)+ geom_line(size=1,show.legend=FALSE)+ scale_color_manual(values=cols)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",size=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",size=0.5)+ scale_x_reverse(limits=c(23,57.4))+ labs(title="Sampled diversity, USA",x="Age (Myr) (data binned per MP 6 to 30)",y="Number of taxa")+ theme(axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) variance_USA <- ggplot(bm_df_to_plot_USA,aes(x=Age,y=variance,color=clade,width=0.35))+ geom_rect(data=strips_USA, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf), fill='white', alpha=0.7,inherit.aes=FALSE)+ geom_line(size=1,show.legend=FALSE)+ scale_color_manual(values=cols)+ geom_vline(xintercept=33.9, linetype="longdash",color="black",size=0.5)+ geom_vline(xintercept=40, linetype="longdash",color="black",size=0.5)+ scale_x_reverse(limits=c(23,57.4))+ labs(title="Body mass variance, USA",x="Age (Myr) (data binned per MP 6 to 30)",y="Variance")+ theme(axis.text=element_text(size=6),axis.title=element_text(size=6),axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #==========================================================================# #============ Tests for non-phylogenetic regime shifts, USA ===============# #==========================================================================# ###Regime shifts using serialised statistic: two-sample Kolmogorov–Smirnov test #Create blank matrix KS_USA <- matrix(nrow=dim(bm_bin_USA)[1],ncol=8); rownames(KS_USA) <- rownames(bm_bin_USA); colnames(KS_USA) <- c("All_D","All_p-value","All_Carnivoramorpha_D","All_Carnivoramorpha_p-value","Non_carnivoran_Carnivoramorpha_D","Non_carnivoran_Carnivoramorpha_p-value","Hyaenodonta_D","Hyaenodonta_p-value") #All taxa for (i in 1:(nrow(nalma))-1){ if (sum(Juhn_df$Max_Age>nalma[i,"Base"])>1){ test <- ks.test(Juhn_df[Juhn_df$Max_Age>nalma[i,"Base"],"Log10BM"],Juhn_df[Juhn_df$Max_Agenalma[i,"Base"])>1){ test <- ks.test(Juhn_df[(Juhn_df$Clade=="Carnivora" | Juhn_df$Clade=="Carnivoramorpha") & Juhn_df$Max_Age>nalma[i,"Base"],"Log10BM"],Juhn_df[(Juhn_df$Clade=="Carnivora" | Juhn_df$Clade=="Carnivoramorpha") & Juhn_df$Max_Agenalma[i,"Base"])>1){ if (sum(Juhn_df$Clade=="Carnivoramorpha" & Juhn_df$Max_Age1){ test <- ks.test(Juhn_df[Juhn_df$Clade=="Carnivoramorpha" & Juhn_df$Max_Age>nalma[i,"Base"],"Log10BM"],Juhn_df[Juhn_df$Clade=="Carnivoramorpha" & Juhn_df$Max_Agenalma[i,"Base"])>1){ if (sum(Juhn_df$Clade=="Hyaenodonta" & Juhn_df$Max_Age1){ test <- ks.test(Juhn_df[Juhn_df$Clade=="Hyaenodonta" & Juhn_df$Max_Age>nalma[i,"Base"],"Log10BM"],Juhn_df[Juhn_df$Clade=="Hyaenodonta" & Juhn_df$Max_Age