--- title: "3. Major Axis - Ontogeny vs phylogeny" author: "Nathalie Feiner" date: "16/01/2021" output: html_document: toc: true toc_depth: 3 toc_float: true --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(message = FALSE) ``` The input dataset contains 374 individual embryos belonging to 15 Anolis species and 693 individual adults belonging to 267 species. Note that relative pathes have to be adjusted for reproducing the results. Note that analyses including PCA on different species were performed using two alternative methods: -using phylogenetic PCA (pPCA; results presented in the main text) -using standard PCA (results presented in the electronic supplementary material) See extended methods section in the electronic supplementary material for more explanations. **This output is based on bootstrapping with 100 iterations; Results presented in the manuscript are based on bootstrapping with 1000 iterations!** ```{r load libraries, include=F} library(ggplot2) library(openxlsx) library(pcaMethods) library(tidyr) library(dplyr) library(naniar) library(VIM) library(missMDA) library(car) library(tinytex) library(phangorn) library(geiger) library(phylobase) library(ggtree) library(phytools) library(lmerTest) library(NISTunits) library(multigroup) library(gtools) library(grid) #load workspace with previous results #load(file = "C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/OntoPhylo3.RData") n.sim <- 100 #shorthands for individual limb bone of fore- and hindlimbs Forel <- c("humerus", "ulna", "digit_1", "digit_2", "digit_3", "digit_4") n.Forel <- length(Forel) Hindl <- c("femur", "tibia", "toe_1", "toe_2", "toe_3", "toe_4", "toe_5") n.Hindl <- length(Hindl) PCS6 <- c("PC1","PC2","PC3","PC4","PC5","PC6") PCS7 <- c("PC1","PC2","PC3","PC4","PC5","PC6","PC7") #load embryonic data Data <- read.csv("C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/Data_Embryos_Final.csv", row.names = 1) #log all data Data[,c(Forel,Hindl)] <- log(Data[,c(Forel,Hindl)]) #load adult data Adult_data_ind <- read.csv("C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/Data_Adult_Individuals.csv", header=T, row.names = 1, na.strings=c("","NA")) Adult_data_spec <- aggregate(Adult_data_ind[,c(2:15)], list(Adult_data_ind$species), mean) colnames(Adult_data_spec)[1] <- "species" Adult_data <- merge(x = Adult_data_spec, y = distinct(Adult_data_ind[ , c("species", "Region","Region_Det")]), by.x = "species", by.y = "species") rownames(Adult_data) <- Adult_data[,1] Adult_data[,c(Forel,Hindl)] <- log(Adult_data[,c(Forel,Hindl)]) Adult_data_selec <- subset(Adult_data, species %in% Data$species) Adult_data_Ant <- subset(Adult_data, Region_Det %in% c("Greater Antilles","Lesser Antilles Northern")) Adult_data_MLpri <- subset(Adult_data, Region_Det == "MLpri") Adult_data_MLsec <- subset(Adult_data, Region_Det == "MLsec") Adult_data_selec$species <- droplevels(Adult_data_selec$species) Adult_data_Ant$species <- droplevels(Adult_data_Ant$species) Adult_data_MLpri$species <- droplevels(Adult_data_MLpri$species) Adult_data_MLsec$species <- droplevels(Adult_data_MLsec$species) #load Phylogenetic tree tree <- read.nexus("C:/Users/Nathalie/Dropbox/ANOLIS/Analyses_EvoRates/Anolis.nex") tree_selec <- drop.tip(tree, name.check(tree, Adult_data_selec)$tree_not_data) tree_Ant <- drop.tip(tree, name.check(tree, Adult_data_Ant)$tree_not_data) tree_MLpri <- drop.tip(tree, name.check(tree, Adult_data_MLpri)$tree_not_data) tree_MLsec <- drop.tip(tree, name.check(tree, Adult_data_MLsec)$tree_not_data) #load functions source("C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/FinalR_scripts/functions.R") ``` # Part I: Evolutionary allometry between hatchlings and adults of the 15 species (Table S9) ### Observed angles ```{r, echo=F} Data_hatch <- subset(Data, stage == "hatch") Data_hatch <- aggregate(Data_hatch[,c(Forel,Hindl)], list(Data_hatch$species), mean) colnames(Data_hatch)[1] <- "species" rownames(Data_hatch) <- Data_hatch[,1] tree_hatch <- drop.tip(tree, name.check(tree, Data_hatch)$tree_not_data) res.pPCA_selec_FL <- phyl.pca(tree_selec, Adult_data_selec[,c(Forel)], method="lambda", mode="cov") res.pPCA_selec_HL <- phyl.pca(tree_selec, Adult_data_selec[,c(Hindl)], method="lambda", mode="cov") res.pPCA_hatch_FL <- phyl.pca(tree_selec, Data_hatch[,c(Forel)], method="lambda", mode="cov") res.pPCA_hatch_HL <- phyl.pca(tree_selec, Data_hatch[,c(Hindl)], method="lambda", mode="cov") res.PCA_selec_FL <- prcomp(Adult_data_selec[,c(Forel)]) res.PCA_selec_HL <- prcomp(Adult_data_selec[,c(Hindl)]) res.PCA_hatch_FL <- prcomp(Data_hatch[,c(Forel)]) res.PCA_hatch_HL <- prcomp(Data_hatch[,c(Hindl)]) paste("observed angle FL (phyloPCA)") angle(res.pPCA_selec_FL$Evec[,1],res.pPCA_hatch_FL$Evec[,1]) paste("observed angle HL (phyloPCA)") angle(res.pPCA_selec_HL$Evec[,1],res.pPCA_hatch_HL$Evec[,1]) paste("observed angle FL (standard PCA)") angle(res.PCA_selec_FL$rotation[,1],res.PCA_hatch_FL$rotation[,1]) paste("observed angle HL (standard PCA)") angle(res.PCA_selec_HL$rotation[,1],res.PCA_hatch_HL$rotation[,1]) ``` ### Create datasets with shared allometry (rotate original datasets) ```{r, echo=F} Data_hatch_new <- Data_hatch_new2 <- Data_hatch[-c(1:nrow(Data_hatch)),] Adult_data_selec_new <- Adult_data_selec_new2 <- Adult_data_selec[-c(1:nrow(Adult_data_selec)),] Data_hatch$stage <- "hatch" Adult_data_selec$stage <- "adult" Data_onto_phylo <- smartbind(Data_hatch,Adult_data_selec) res.FCPCA_ontphy_FL <- FCPCA(Data_onto_phylo[,c(Forel)], Data_onto_phylo$stage, graph=F) res.FCPCA_ontphy_HL <- FCPCA(Data_onto_phylo[,c(Hindl)], Data_onto_phylo$stage, graph=F) #Adults X <- Adult_data_selec FL <- X[,c(Forel)] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_selec_FL$a,nrow(Adult_data_selec)),nrow(Adult_data_selec),6,byrow=TRUE)) %*% res.pPCA_selec_FL$Evec %*% t(res.FCPCA_ontphy_FL$loadings.common) HL <- X[,c(Hindl)] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_selec_HL$a,nrow(Adult_data_selec)),nrow(Adult_data_selec),7,byrow=TRUE)) %*% res.pPCA_selec_HL$Evec %*% t(res.FCPCA_ontphy_HL$loadings.common) X[,c(Forel)] <- FL_new X[,c(Hindl)] <- HL_new Adult_data_selec_new <- rbind(Adult_data_selec_new,X) #Hatchlings X <- Data_hatch FL <- X[,c(Forel)] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_hatch_FL$a,nrow(Data_hatch)),nrow(Data_hatch),6,byrow=TRUE)) %*% res.pPCA_hatch_FL$Evec %*% t(res.FCPCA_ontphy_FL$loadings.common) HL <- X[,c(Hindl)] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_hatch_HL$a,nrow(Data_hatch)),nrow(Data_hatch),7,byrow=TRUE)) %*% res.pPCA_hatch_HL$Evec %*% t(res.FCPCA_ontphy_HL$loadings.common) X[,c(Forel)] <- FL_new X[,c(Hindl)] <- HL_new Data_hatch_new <- rbind(Data_hatch_new,X) #non-phylogenetic PCA: #Adults X <- Adult_data_selec FL <- X[,c(Forel)] FL_new <- as.matrix(scale(FL, center=T, scale=F)) %*% res.PCA_selec_FL$rotation %*% t(res.FCPCA_ontphy_FL$loadings.common) HL <- X[,c(Hindl)] HL_new <- as.matrix(scale(HL, center=T, scale=F)) %*% res.PCA_selec_HL$rotation %*% t(res.FCPCA_ontphy_HL$loadings.common) X[,c(Forel)] <- FL_new X[,c(Hindl)] <- HL_new Adult_data_selec_new2 <- rbind(Adult_data_selec_new2,X) #Hatchlings X <- Data_hatch FL <- X[,c(Forel)] FL_new <- as.matrix(scale(FL, center=T, scale=F)) %*% res.PCA_hatch_FL$rotation %*% t(res.FCPCA_ontphy_FL$loadings.common) HL <- X[,c(Hindl)] HL_new <- as.matrix(scale(HL, center=T, scale=F)) %*% res.PCA_hatch_HL$rotation %*% t(res.FCPCA_ontphy_HL$loadings.common) X[,c(Forel)] <- FL_new X[,c(Hindl)] <- HL_new Data_hatch_new2 <- rbind(Data_hatch_new2,X) ``` ### Bootstrapping from datasets with shared multivariate allometry ```{r, echo=F, Warning=F} #First, parametric bootstrapping Onto_coef_sim_FL <- Phylo_coef_sim_FL <- array(NA, dim=c(2,n.Forel,n.sim)) Onto_coef_sim_HL <- Phylo_coef_sim_HL <- array(NA, dim=c(2,n.Hindl,n.sim)) sim_rSDE.selec_FL <- sim.ecov(as.matrix(Adult_data_selec_new[tree_selec$tip.label, Forel]), tree_selec, lambda.method = "ML", b = n.sim) sim_rSDE.selec_HL <- sim.ecov(as.matrix(Adult_data_selec_new[tree_selec$tip.label, Hindl]), tree_selec, lambda.method = "ML", b = n.sim) sim_rSDE.hatch_FL <- sim.ecov(as.matrix(Data_hatch_new[tree_selec$tip.label, Forel]), tree_selec, lambda.method = "ML", b = n.sim) sim_rSDE.hatch_HL <- sim.ecov(as.matrix(Data_hatch_new[tree_selec$tip.label, Hindl]), tree_selec, lambda.method = "ML", b = n.sim) Phylo_coef_sim_FL[1,1:6,] <- sim_rSDE.selec_FL$PC1 Phylo_coef_sim_HL[1,1:7,] <- sim_rSDE.selec_HL$PC1 Onto_coef_sim_FL[1,1:6,] <- sim_rSDE.hatch_FL$PC1 Onto_coef_sim_HL[1,1:7,] <- sim_rSDE.hatch_HL$PC1 #Second, non-parametric bootstrapping for(k in 1:n.sim){ Data_adul <- sample_frac(Adult_data_selec_new2, 1, replace=T) Data_onto <- sample_frac(Data_hatch_new2, 1, replace=T) Phylo_coef_sim_FL[2,,k] <- prcomp(Data_adul[,Forel])$rotation[,1] Onto_coef_sim_FL[2,,k] <- prcomp(Data_onto[,Forel])$rotation[,1] Phylo_coef_sim_HL[2,,k] <- prcomp(Data_adul[,Hindl])$rotation[,1] Onto_coef_sim_HL[2,,k] <- prcomp(Data_onto[,Hindl])$rotation[,1] } Angles_FL_sim <- Angles_HL_sim <- array(NA, dim=c(2,1,n.sim)) rownames(Angles_FL_sim) <- rownames(Angles_HL_sim) <- c("pPCA", "standard PCA") for(k in 1:n.sim){ Angles_FL_sim[1,,k] <- angle(abs(Onto_coef_sim_FL[1,,k]),abs(Phylo_coef_sim_FL[1,,k])) Angles_HL_sim[1,,k] <- angle(abs(Onto_coef_sim_HL[1,,k]),abs(Phylo_coef_sim_HL[1,,k])) Angles_FL_sim[2,,k] <- angle(abs(Onto_coef_sim_FL[2,,k]),abs(Phylo_coef_sim_FL[2,,k])) Angles_HL_sim[2,,k] <- angle(abs(Onto_coef_sim_HL[2,,k]),abs(Phylo_coef_sim_HL[2,,k])) } Angles_FL_sim_975CI <- apply(Angles_FL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) Angles_HL_sim_975CI <- apply(Angles_HL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) paste("95% CI of simulated (parallel) angles between onto and phylo axis (15 species) - forelimb") Angles_FL_sim_975CI paste("95% CI of simulated (parallel) angles between onto and phylo axis (15 species) - hindlimb") Angles_HL_sim_975CI ``` # Part II: Are ontogenetic & evolutionary coefficients different from isometric? (Tables S10 and S11) Note that coefficients are given as the raw values divided by p^-0.5^ with p being the number of elements. Thus, values below 1 are negative allometric, above 1 positive allometric (relative to other elements). ```{r echo=F, fig.show="hold", warning=FALSE, out.width="50%"} Data <- subset(Data, species != "roquet") #drop A. roquet Data$species <- droplevels(Data$species) #bootstrap datasets and extract coefficients Ont_Phy_coef_boot_FL <- array(NA, dim=c(7,n.Forel,n.sim)) Ont_Phy_coef_boot_HL <- array(NA, dim=c(7,n.Hindl,n.sim)) rownames(Ont_Phy_coef_boot_FL) <- rownames(Ont_Phy_coef_boot_HL) <- c("Ontogeny","Phylogeny - Ant (pPCA)","Phylogeny - MLpri (pPCA)","Phylogeny - MLsec (pPCA)","Phylogeny - Ant (standard PCA)","Phylogeny - MLpri (standard PCA)","Phylogeny - MLsec (standard PCA)") #observed coefficients: Ont_Phy_coef_boot_FL[1,1:6,1] <- FCPCA(Data[,Forel], Data$species, graph=F)$loadings.common[,1] Ont_Phy_coef_boot_HL[1,1:7,1] <- FCPCA(Data[,Hindl], Data$species, graph=F)$loadings.common[,1] Ont_Phy_coef_boot_FL[2,1:6,1] <- phyl.pca(tree_Ant, Adult_data_Ant[,Forel], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_FL[3,1:6,1] <- phyl.pca(tree_MLpri, Adult_data_MLpri[,Forel], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_FL[4,1:6,1] <- phyl.pca(tree_MLsec, Adult_data_MLsec[,Forel], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_HL[2,1:7,1] <- phyl.pca(tree_Ant, Adult_data_Ant[,Hindl], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_HL[3,1:7,1] <- phyl.pca(tree_MLpri, Adult_data_MLpri[,Hindl], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_HL[4,1:7,1] <- phyl.pca(tree_MLsec, Adult_data_MLsec[,Hindl], method="lambda", mode="cov")$Evec[,1] Ont_Phy_coef_boot_FL[5,1:6,1] <- prcomp(Adult_data_Ant[,Forel])$rotation[,1] Ont_Phy_coef_boot_FL[6,1:6,1] <- prcomp(Adult_data_MLpri[,Forel])$rotation[,1] Ont_Phy_coef_boot_FL[7,1:6,1] <- prcomp(Adult_data_MLsec[,Forel])$rotation[,1] Ont_Phy_coef_boot_HL[5,1:7,1] <- prcomp(Adult_data_Ant[,Hindl])$rotation[,1] Ont_Phy_coef_boot_HL[6,1:7,1] <- prcomp(Adult_data_MLpri[,Hindl])$rotation[,1] Ont_Phy_coef_boot_HL[7,1:7,1] <- prcomp(Adult_data_MLsec[,Hindl])$rotation[,1] #non-parametric bootstrapping for(k in 2:n.sim){ Data_sub <- sample_frac(Data, 1, replace=T) Adult_data_Ant_sub <- sample_frac(Adult_data_Ant, 1, replace=T) Adult_data_MLpri_sub <- sample_frac(Adult_data_MLpri, 1, replace=T) Adult_data_MLsec_sub <- sample_frac(Adult_data_MLsec, 1, replace=T) temp <- Data_sub %>% group_by(species) %>% summarise(Unique_Elements = n_distinct(Individ_ID)) if ( min(temp[,2]) <= 2){ k-1 }else{ Ont_Phy_coef_boot_FL[1,,k] <- FCPCA(Data_sub[,Forel], Data_sub$species, graph=F)$loadings.common[,1] Ont_Phy_coef_boot_HL[1,,k] <- FCPCA(Data_sub[,Hindl], Data_sub$species, graph=F)$loadings.common[,1] Ont_Phy_coef_boot_FL[5,,k] <- prcomp(Adult_data_Ant_sub[,Forel])$rotation[,1] Ont_Phy_coef_boot_FL[6,,k] <- prcomp(Adult_data_MLpri_sub[,Forel])$rotation[,1] Ont_Phy_coef_boot_FL[7,,k] <- prcomp(Adult_data_MLsec_sub[,Forel])$rotation[,1] Ont_Phy_coef_boot_HL[5,,k] <- prcomp(Adult_data_Ant_sub[,Hindl])$rotation[,1] Ont_Phy_coef_boot_HL[6,,k] <- prcomp(Adult_data_MLpri_sub[,Hindl])$rotation[,1] Ont_Phy_coef_boot_HL[7,,k] <- prcomp(Adult_data_MLsec_sub[,Hindl])$rotation[,1] } } #parametric bootstrapping of phylogenetic pPC1 Ont_Phy_coef_boot_FL[2,1:6,2:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant[tree_Ant$tip.label, Forel]), tree_Ant, lambda.method = "ML", b = n.sim-1)$PC1 Ont_Phy_coef_boot_FL[3,1:6,2:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri[tree_MLpri$tip.label, Forel]), tree_MLpri, lambda.method = "ML", b = n.sim-1)$PC1 Ont_Phy_coef_boot_FL[4,1:6,2:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec[tree_MLsec$tip.label, Forel]), tree_MLsec, lambda.method = "ML", b = n.sim-1)$PC1 Ont_Phy_coef_boot_HL[2,1:7,2:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant[tree_Ant$tip.label, Hindl]), tree_Ant, lambda.method = "ML", b = n.sim-1)$PC1 Ont_Phy_coef_boot_HL[3,1:7,2:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri[tree_MLpri$tip.label, Hindl]), tree_MLpri, lambda.method = "ML", b = n.sim-1)$PC1 Ont_Phy_coef_boot_HL[4,1:7,2:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec[tree_MLsec$tip.label, Hindl]), tree_MLsec, lambda.method = "ML", b = n.sim-1)$PC1 paste("Forelimb observed coefficients") round(t(abs(Ont_Phy_coef_boot_FL[,1:6,1]))/6^-0.5,3) paste("Forelimb lower CI") round(t(apply(abs(Ont_Phy_coef_boot_FL), c(1,2), quantile, na.rm=TRUE, probs=0.025))/6^-0.5,3) paste("Forelimb upper CI") round(t(apply(abs(Ont_Phy_coef_boot_FL), c(1,2), quantile, na.rm=TRUE, probs=0.975))/6^-0.5,3) paste("Hindlimb observed coefficients") round(t(abs(Ont_Phy_coef_boot_HL[,1:7,1]))/7^-0.5,3) paste("Hindlimb lower CI") round(t(apply(abs(Ont_Phy_coef_boot_HL), c(1,2), quantile, na.rm=TRUE, probs=0.025))/7^-0.5,3) paste("Hindlimb upper CI") round(t(apply(abs(Ont_Phy_coef_boot_HL), c(1,2), quantile, na.rm=TRUE, probs=0.975))/7^-0.5,3) ``` # Part III: Alignment between ontogenetic and evolutionary allometries ### Observed angles ```{r, echo=F, Warning=F} paste("Angle between ontogenetic common axis and phylo axis - forelimb") Angles_FL_obs <- Angles_HL_obs <- array(NA, dim=c(6,1)) rownames(Angles_FL_obs) <- rownames(Angles_HL_obs) <- c("Onto - Ant (pPCA)","Onto - MLpri (pPCA)","Onto - MLsec (pPCA)","Onto - Ant (standard PCA)","Onto - MLpri (standard PCA)","Onto - MLsec (standard PCA)") paste("Observed angles between ontogenetic common axis and evolutionary axis - forelimb") Angles_FL_obs[1,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[2,1:6,1])) #7.35 Angles_FL_obs[4,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[5,1:6,1])) #6.85 Angles_FL_obs[2,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[3,1:6,1])) #8.62 Angles_FL_obs[5,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[6,1:6,1])) #8.71 Angles_FL_obs[3,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[4,1:6,1])) #8.11 Angles_FL_obs[6,1] <- angle(abs(Ont_Phy_coef_boot_FL[1,1:6,1]),abs(Ont_Phy_coef_boot_FL[7,1:6,1])) #8.16 Angles_FL_obs paste("Observed angles between ontogenetic common axis and evolutionary axis - hindlimb") Angles_HL_obs[1,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[2,1:7,1])) #7.32 Angles_HL_obs[4,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[5,1:7,1])) #7.64 Angles_HL_obs[2,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[3,1:7,1])) #7.46 Angles_HL_obs[5,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[6,1:7,1])) #8.22 Angles_HL_obs[3,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[4,1:7,1])) #5.28 Angles_HL_obs[6,1] <- angle(abs(Ont_Phy_coef_boot_HL[1,1:7,1]),abs(Ont_Phy_coef_boot_HL[7,1:7,1])) #5.39 Angles_HL_obs ``` ### Create dataset with shared allometry (rotate original datasets) ```{r, echo=F, Warning=F} #run common PCA on onto & phylo together #fuse datasets Adult_data_Ant_noSpec <- Adult_data_Ant Adult_data_Ant_noSpec$species <- "adults" Adult_data_MLpri_noSpec <- Adult_data_MLpri Adult_data_MLpri_noSpec$species <- "adults" Adult_data_MLsec_noSpec <- Adult_data_MLsec Adult_data_MLsec_noSpec$species <- "adults" Data_noSpec <- Data Data_noSpec$species <- "embryos" Data_onto_phylo_Ant <- smartbind(Data_noSpec,Adult_data_Ant_noSpec) Data_onto_phylo_MLpri <- smartbind(Data_noSpec,Adult_data_MLpri_noSpec) Data_onto_phylo_MLsec <- smartbind(Data_noSpec,Adult_data_MLsec_noSpec) res.FCPCA_ontphy_Ant_FL <- FCPCA(Data_onto_phylo_Ant[,Forel], Data_onto_phylo_Ant$species, graph=F) res.FCPCA_ontphy_Ant_HL <- FCPCA(Data_onto_phylo_Ant[,Hindl], Data_onto_phylo_Ant$species, graph=F) res.FCPCA_ontphy_MLpri_FL <- FCPCA(Data_onto_phylo_MLpri[,Forel], Data_onto_phylo_MLpri$species, graph=F) res.FCPCA_ontphy_MLpri_HL <- FCPCA(Data_onto_phylo_MLpri[,Hindl], Data_onto_phylo_MLpri$species, graph=F) res.FCPCA_ontphy_MLsec_FL <- FCPCA(Data_onto_phylo_MLsec[,Forel], Data_onto_phylo_MLsec$species, graph=F) res.FCPCA_ontphy_MLsec_HL <- FCPCA(Data_onto_phylo_MLsec[,Hindl], Data_onto_phylo_MLsec$species, graph=F) #ontogenetic cPCs res.FCPCA_FL <- FCPCA(Data[,Forel], Data$species, graph=F) res.FCPCA_HL <- FCPCA(Data[,Hindl], Data$species, graph=F) #phylogenetic pPCs res.pPCA_Ant_FL <- phyl.pca(tree_Ant, Adult_data_Ant[,Forel], method="lambda", mode="cov") res.pPCA_MLpri_FL <- phyl.pca(tree_MLpri, Adult_data_MLpri[,Forel], method="lambda", mode="cov") res.pPCA_MLsec_FL <- phyl.pca(tree_MLsec, Adult_data_MLsec[,Forel], method="lambda", mode="cov") res.pPCA_Ant_HL <- phyl.pca(tree_Ant, Adult_data_Ant[,Hindl], method="lambda", mode="cov") res.pPCA_MLpri_HL <- phyl.pca(tree_MLpri, Adult_data_MLpri[,Hindl], method="lambda", mode="cov") res.pPCA_MLsec_HL <- phyl.pca(tree_MLsec, Adult_data_MLsec[,Hindl], method="lambda", mode="cov") #non-phylogenetic pPCs res.PCA_Ant_FL <- prcomp(Adult_data_Ant[,Forel]) res.PCA_MLpri_FL <- prcomp(Adult_data_MLpri[,Forel]) res.PCA_MLsec_FL <- prcomp(Adult_data_MLsec[,Forel]) res.PCA_Ant_HL <- prcomp(Adult_data_Ant[,Hindl]) res.PCA_MLpri_HL <- prcomp(Adult_data_MLpri[,Hindl]) res.PCA_MLsec_HL <- prcomp(Adult_data_MLsec[,Hindl]) Data_Ant_new <- Data_MLpri_new <- Data_MLsec_new <- Data[-c(1:nrow(Data)),] Adult_data_Ant_new <- Adult_data_Ant_new2 <- Adult_data_Ant[-c(1:nrow(Adult_data_Ant)),] Adult_data_MLpri_new <- Adult_data_MLpri_new2 <- Adult_data_MLpri[-c(1:nrow(Adult_data_MLpri)),] Adult_data_MLsec_new <- Adult_data_MLsec_new2 <- Adult_data_MLsec[-c(1:nrow(Adult_data_MLsec)),] ##Ant #Create dataset with shared allometry for onto: X <- Data FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.FCPCA_FL$loadings.common %*% t(res.FCPCA_ontphy_Ant_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.FCPCA_HL$loadings.common %*% t(res.FCPCA_ontphy_Ant_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Data_Ant_new <- rbind(Data_Ant_new,X) #Create dataset with shared allometry for phylo (pPCA): X <- Adult_data_Ant FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_Ant_FL$a,nrow(Adult_data_Ant)),nrow(Adult_data_Ant),6,byrow=TRUE)) %*% res.pPCA_Ant_FL$Evec %*% t(res.FCPCA_ontphy_Ant_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_Ant_HL$a,nrow(Adult_data_Ant)),nrow(Adult_data_Ant),7,byrow=TRUE)) %*% res.pPCA_Ant_HL$Evec %*% t(res.FCPCA_ontphy_Ant_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_Ant_new <- rbind(Adult_data_Ant_new,X) #Create dataset with shared allometry for phylo (standard PCA): X <- Adult_data_Ant FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_Ant_FL$rotation %*% t(res.FCPCA_ontphy_Ant_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_Ant_HL$rotation %*% t(res.FCPCA_ontphy_Ant_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_Ant_new2 <- rbind(Adult_data_Ant_new2,X) ##MLpri #Create dataset with shared allometry for onto: X <- Data FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.FCPCA_FL$loadings.common %*% t(res.FCPCA_ontphy_MLpri_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.FCPCA_HL$loadings.common %*% t(res.FCPCA_ontphy_MLpri_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Data_MLpri_new <- rbind(Data_MLpri_new,X) #Create dataset with shared allometry for phylo (pPCA): X <- Adult_data_MLpri FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_MLpri_FL$a,nrow(Adult_data_MLpri)),nrow(Adult_data_MLpri),6,byrow=TRUE)) %*% res.pPCA_MLpri_FL$Evec %*% t(res.FCPCA_ontphy_MLpri_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_MLpri_HL$a,nrow(Adult_data_MLpri)),nrow(Adult_data_MLpri),7,byrow=TRUE)) %*% res.pPCA_MLpri_HL$Evec %*% t(res.FCPCA_ontphy_MLpri_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLpri_new <- rbind(Adult_data_MLpri_new,X) #Create dataset with shared allometry for phylo (standard PCA): X <- Adult_data_MLpri FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_MLpri_FL$rotation %*% t(res.FCPCA_ontphy_MLpri_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_MLpri_HL$rotation %*% t(res.FCPCA_ontphy_MLpri_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLpri_new2 <- rbind(Adult_data_MLpri_new2,X) ##MLsec #Create dataset with shared allometry for onto: X <- Data FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.FCPCA_FL$loadings.common %*% t(res.FCPCA_ontphy_MLsec_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.FCPCA_HL$loadings.common %*% t(res.FCPCA_ontphy_MLsec_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Data_MLsec_new <- rbind(Data_MLsec_new,X) #Create dataset with shared allometry for phylo (pPCA): X <- Adult_data_MLsec FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_MLsec_FL$a,nrow(Adult_data_MLsec)),nrow(Adult_data_MLsec),6,byrow=TRUE)) %*% res.pPCA_MLsec_FL$Evec %*% t(res.FCPCA_ontphy_MLsec_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_MLsec_HL$a,nrow(Adult_data_MLsec)),nrow(Adult_data_MLsec),7,byrow=TRUE)) %*% res.pPCA_MLsec_HL$Evec %*% t(res.FCPCA_ontphy_MLsec_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLsec_new <- rbind(Adult_data_MLsec_new,X) #Create dataset with shared allometry for phylo (standard PCA): X <- Adult_data_MLsec FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_MLsec_FL$rotation %*% t(res.FCPCA_ontphy_MLsec_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_MLsec_HL$rotation %*% t(res.FCPCA_ontphy_MLsec_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLsec_new2 <- rbind(Adult_data_MLsec_new2,X) ``` ### Bootstrapping from datasets with shared multivariate allometry ```{r, echo=F} Ont_Phy_coef_sim_FL <- array(NA, dim=c(9,n.Forel,n.sim)) Ont_Phy_coef_sim_HL <- array(NA, dim=c(9,n.Hindl,n.sim)) rownames(Ont_Phy_coef_sim_FL) <- rownames(Ont_Phy_coef_sim_HL) <- c("Ontogeny - Ant","Ontogeny - MLpri","Ontogeny - MLsec","Phylogeny - Ant (pPCA)","Phylogeny - MLpri (pPCA)","Phylogeny - MLsec (pPCA)","Phylogeny - Ant (standard PCA)","Phylogeny - MLpri (standard PCA)","Phylogeny - MLsec (standard PCA)") for(k in 1:n.sim){ Data_Ant_new_sub <- sample_frac(Data_Ant_new, 1, replace=T) temp <- Data_Ant_new_sub %>% group_by(species) %>% summarise(Unique_Elements = n_distinct(Individ_ID)) if ( min(temp[,2]) <= 2){ k-1 }else{ Ont_Phy_coef_sim_FL[1,,k] <- FCPCA(Data_Ant_new_sub[,Forel], Data_Ant_new_sub$species, graph=F)$loadings.common[,1] Ont_Phy_coef_sim_HL[1,,k] <- FCPCA(Data_Ant_new_sub[,Hindl], Data_Ant_new_sub$species, graph=F)$loadings.common[,1] } } for(k in 1:n.sim){ Data_MLpri_new_sub <- sample_frac(Data_MLpri_new, 1, replace=T) temp <- Data_MLpri_new_sub %>% group_by(species) %>% summarise(Unique_Elements = n_distinct(Individ_ID)) if ( min(temp[,2]) <= 2){ k-1 }else{ Ont_Phy_coef_sim_FL[2,,k] <- FCPCA(Data_MLpri_new_sub[,Forel], Data_MLpri_new_sub$species, graph=F)$loadings.common[,1] Ont_Phy_coef_sim_HL[2,,k] <- FCPCA(Data_MLpri_new_sub[,Hindl], Data_MLpri_new_sub$species, graph=F)$loadings.common[,1] } } for(k in 1:n.sim){ Data_MLsec_new_sub <- sample_frac(Data_MLsec_new, 1, replace=T) temp <- Data_MLsec_new_sub %>% group_by(species) %>% summarise(Unique_Elements = n_distinct(Individ_ID)) if ( min(temp[,2]) <= 2){ k-1 }else{ Ont_Phy_coef_sim_FL[3,,k] <- FCPCA(Data_MLsec_new_sub[,Forel], Data_MLsec_new_sub$species, graph=F)$loadings.common[,1] Ont_Phy_coef_sim_HL[3,,k] <- FCPCA(Data_MLsec_new_sub[,Hindl], Data_MLsec_new_sub$species, graph=F)$loadings.common[,1] } } for(k in 1:n.sim){ Adult_data_Ant_new2_sub <- sample_frac(Adult_data_Ant_new2, 1, replace=T) Adult_data_MLpri_new2_sub <- sample_frac(Adult_data_MLpri_new2, 1, replace=T) Adult_data_MLsec_new2_sub <- sample_frac(Adult_data_MLsec_new2, 1, replace=T) Ont_Phy_coef_sim_FL[7,,k] <- prcomp(Adult_data_Ant_new2_sub[,Forel])$rotation[,1] Ont_Phy_coef_sim_HL[7,,k] <- prcomp(Adult_data_Ant_new2_sub[,Hindl])$rotation[,1] Ont_Phy_coef_sim_FL[8,,k] <- prcomp(Adult_data_MLpri_new2_sub[,Forel])$rotation[,1] Ont_Phy_coef_sim_HL[8,,k] <- prcomp(Adult_data_MLpri_new2_sub[,Hindl])$rotation[,1] Ont_Phy_coef_sim_FL[9,,k] <- prcomp(Adult_data_MLsec_new2_sub[,Forel])$rotation[,1] Ont_Phy_coef_sim_HL[9,,k] <- prcomp(Adult_data_MLsec_new2_sub[,Hindl])$rotation[,1] } Ont_Phy_coef_sim_FL[4,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant_new[tree_Ant$tip.label, Forel]), tree_Ant, lambda.method = "ML", b = n.sim)$PC1 Ont_Phy_coef_sim_FL[5,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri_new[tree_MLpri$tip.label, Forel]), tree_MLpri, lambda.method = "ML", b = n.sim)$PC1 Ont_Phy_coef_sim_FL[6,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec_new[tree_MLsec$tip.label, Forel]), tree_MLsec, lambda.method = "ML", b = n.sim)$PC1 Ont_Phy_coef_sim_HL[4,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant_new[tree_Ant$tip.label, Hindl]), tree_Ant, lambda.method = "ML", b = n.sim)$PC1 Ont_Phy_coef_sim_HL[5,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri_new[tree_MLpri$tip.label, Hindl]), tree_MLpri, lambda.method = "ML", b = n.sim)$PC1 Ont_Phy_coef_sim_HL[6,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec_new[tree_MLsec$tip.label, Hindl]), tree_MLsec, lambda.method = "ML", b = n.sim)$PC1 Angles_FL_sim <- Angles_HL_sim <- array(NA, dim=c(6,1,n.sim)) rownames(Angles_FL_sim) <- rownames(Angles_HL_sim) <- c("Ont-Ant (pPCA)","Ont-MLpri (pPCA)","Ont-MLsec (pPCA)","Ont-Ant (standard PCA)","Ont-MLpri (standard PCA)","Ont-MLsec (standard PCA)") for(k in 1:n.sim){ Angles_FL_sim[1,,k] <- angle(abs(Ont_Phy_coef_sim_FL[1,,k]),abs(Ont_Phy_coef_sim_FL[4,,k])) Angles_HL_sim[1,,k] <- angle(abs(Ont_Phy_coef_sim_HL[1,,k]),abs(Ont_Phy_coef_sim_HL[4,,k])) Angles_FL_sim[2,,k] <- angle(abs(Ont_Phy_coef_sim_FL[2,,k]),abs(Ont_Phy_coef_sim_FL[5,,k])) Angles_HL_sim[2,,k] <- angle(abs(Ont_Phy_coef_sim_HL[2,,k]),abs(Ont_Phy_coef_sim_HL[5,,k])) Angles_FL_sim[3,,k] <- angle(abs(Ont_Phy_coef_sim_FL[3,,k]),abs(Ont_Phy_coef_sim_FL[6,,k])) Angles_HL_sim[3,,k] <- angle(abs(Ont_Phy_coef_sim_HL[3,,k]),abs(Ont_Phy_coef_sim_HL[6,,k])) Angles_FL_sim[4,,k] <- angle(abs(Ont_Phy_coef_sim_FL[1,,k]),abs(Ont_Phy_coef_sim_FL[7,,k])) Angles_HL_sim[4,,k] <- angle(abs(Ont_Phy_coef_sim_HL[1,,k]),abs(Ont_Phy_coef_sim_HL[7,,k])) Angles_FL_sim[5,,k] <- angle(abs(Ont_Phy_coef_sim_FL[2,,k]),abs(Ont_Phy_coef_sim_FL[8,,k])) Angles_HL_sim[5,,k] <- angle(abs(Ont_Phy_coef_sim_HL[2,,k]),abs(Ont_Phy_coef_sim_HL[8,,k])) Angles_FL_sim[6,,k] <- angle(abs(Ont_Phy_coef_sim_FL[3,,k]),abs(Ont_Phy_coef_sim_FL[9,,k])) Angles_HL_sim[6,,k] <- angle(abs(Ont_Phy_coef_sim_HL[3,,k]),abs(Ont_Phy_coef_sim_HL[9,,k])) } Angles_FL_sim_95CI <- apply(Angles_FL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) Angles_HL_sim_95CI <- apply(Angles_HL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) paste("95% CI of simulated (parallel) angles between onto and phylo axis - forelimb") Angles_FL_sim_95CI paste("95% CI of simulated (parallel) angles between onto and phylo axis - hindlimb") Angles_HL_sim_95CI #get P-values P_theta_ontphy_FL <- P_theta_ontphy_HL <- array(NA, dim=c(6,1)) rownames(P_theta_ontphy_FL) <- rownames(P_theta_ontphy_HL) <- c("Ont-Ant (pPCA)","Ont-MLpri (pPCA)","Ont-MLsec (pPCA)","Ont-Ant (standard PCA)","Ont-MLpri (standard PCA)","Ont-MLsec (standard PCA)") for(i in 1:6){ P_theta_ontphy_FL[i,1] <- 1-sum(as.vector(na.omit(Angles_FL_sim[i,1,1:n.sim])) < Angles_FL_obs[i,1])/length(as.vector(na.omit(Angles_FL_sim[i,1,1:n.sim]))) P_theta_ontphy_HL[i,1] <- 1-sum(as.vector(na.omit(Angles_HL_sim[i,1,1:n.sim])) < Angles_HL_obs[i,1])/length(as.vector(na.omit(Angles_HL_sim[i,1,1:n.sim]))) } print("P-values and angles theta - FL") P_theta_ontphy_FL print("P-values and angles theta - HL") P_theta_ontphy_HL ``` ### Plotting (Figure 4 and S3) ```{r, echo=F, Warning=F} plots <- list() theta_FL <- as.data.frame(t(Angles_FL_sim[1:6,,1:n.sim])) theta_HL <- as.data.frame(t(Angles_HL_sim[1:6,,1:n.sim])) colnames(theta_FL) <- colnames(theta_HL) <- c("Ant","MLpri","MLsec","Ant2","MLpri2","MLsec2") paste("Method: pPCA") paste("top row: forelimb") paste("bottom row: hindlimb") plots[[1]] <- ggplot(theta_FL, aes(x=Ant)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[1,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[1,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[2]] <- ggplot(theta_FL, aes(x=MLpri)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[2,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[2,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[3]] <- ggplot(theta_FL, aes(x=MLsec)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[3,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[3,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[4]] <- ggplot(theta_HL, aes(x=Ant)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[1,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[1,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[5]] <- ggplot(theta_HL, aes(x=MLpri)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[2,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[2,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[6]] <- ggplot(theta_HL, aes(x=MLsec)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[3,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[2,1], size=1, colour="black") + xlim(0,10) + theme_bw() paste("Method: standard PCA") paste("top row: forelimb") paste("bottom row: hindlimb") plots[[7]] <- ggplot(theta_FL, aes(x=Ant2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[4,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[4,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[8]] <- ggplot(theta_FL, aes(x=MLpri2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[5,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[5,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[9]] <- ggplot(theta_FL, aes(x=MLsec2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_FL_sim_95CI[6,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_FL_obs[6,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[10]] <- ggplot(theta_HL, aes(x=Ant2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[4,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[4,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[11]] <- ggplot(theta_HL, aes(x=MLpri2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[5,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[5,1], size=1, colour="black") + xlim(0,10) + theme_bw() plots[[12]] <- ggplot(theta_HL, aes(x=MLsec2)) + geom_histogram(aes(y = ..count..), binwidth = 0.15, fill = "darkgrey") + geom_vline(xintercept = Angles_HL_sim_95CI[6,], size=.4, colour="blue")+ geom_vline(xintercept = Angles_HL_obs[6,1], size=1, colour="black") + xlim(0,10) + theme_bw() #pdf("C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/Plots/OntoPhylo2.pdf", width=10, height=5) grid.newpage() pushViewport(viewport(layout = grid.layout(nrow = 2, ncol = 3))) # A helper function to define a region on the layout define_region <- function(row, col){ viewport(layout.pos.row = row, layout.pos.col = col) } # Arrange the plots print(plots[[1]], vp = define_region(row = 1, col = 1)) print(plots[[2]], vp = define_region(row = 1, col = 2)) print(plots[[3]], vp = define_region(row = 1, col = 3)) print(plots[[4]], vp = define_region(row = 2, col = 1)) print(plots[[5]], vp = define_region(row = 2, col = 2)) print(plots[[6]], vp = define_region(row = 2, col = 3)) #dev.off() #pdf("C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/Plots/OntoPhylo_alternative2.pdf", width=10, height=5) #change name! grid.newpage() pushViewport(viewport(layout = grid.layout(nrow = 2, ncol = 3))) # Arrange the plots print(plots[[7]], vp = define_region(row = 1, col = 1)) print(plots[[8]], vp = define_region(row = 1, col = 2)) print(plots[[9]], vp = define_region(row = 1, col = 3)) print(plots[[10]], vp = define_region(row = 2, col = 1)) print(plots[[11]], vp = define_region(row = 2, col = 2)) print(plots[[12]], vp = define_region(row = 2, col = 3)) #dev.off() ``` # Part IV: Alignment among evolutionary allometries of the three groups (Ant, MLpri and MLsec) ### Observed angles ```{r, echo=F, Warning=F} paste("Angle between phylo axis - forelimb") Angles_phy_FL_obs <- Angles_phy_HL_obs <- array(NA, dim=c(6,1)) rownames(Angles_phy_FL_obs) <- rownames(Angles_phy_HL_obs) <- c("Ant - MLpri (pPCA)","Ant - MLsec (pPCA)","MLpri - MLsec (pPCA)","Ant - MLpri (standard PCA)","Ant - MLsec (standard PCA)","MLpri - MLsec (standard PCA)") paste("Observed angles between evolutionary axis - forelimb") Angles_phy_FL_obs[1,1] <- angle(abs(Ont_Phy_coef_boot_FL[2,1:6,1]),abs(Ont_Phy_coef_boot_FL[3,1:6,1])) #4.44 Angles_phy_FL_obs[4,1] <- angle(abs(Ont_Phy_coef_boot_FL[5,1:6,1]),abs(Ont_Phy_coef_boot_FL[6,1:6,1])) #4.60 Angles_phy_FL_obs[2,1] <- angle(abs(Ont_Phy_coef_boot_FL[2,1:6,1]),abs(Ont_Phy_coef_boot_FL[4,1:6,1])) #2.88 Angles_phy_FL_obs[5,1] <- angle(abs(Ont_Phy_coef_boot_FL[5,1:6,1]),abs(Ont_Phy_coef_boot_FL[7,1:6,1])) #2.68 Angles_phy_FL_obs[3,1] <- angle(abs(Ont_Phy_coef_boot_FL[3,1:6,1]),abs(Ont_Phy_coef_boot_FL[4,1:6,1])) #5.35 Angles_phy_FL_obs[6,1] <- angle(abs(Ont_Phy_coef_boot_FL[6,1:6,1]),abs(Ont_Phy_coef_boot_FL[7,1:6,1])) #5.89 Angles_phy_FL_obs paste("Observed angles between evolutionary axis - hindlimb") Angles_phy_HL_obs[1,1] <- angle(abs(Ont_Phy_coef_boot_HL[2,1:7,1]),abs(Ont_Phy_coef_boot_HL[3,1:7,1])) #4.63 Angles_phy_HL_obs[4,1] <- angle(abs(Ont_Phy_coef_boot_HL[5,1:7,1]),abs(Ont_Phy_coef_boot_HL[6,1:7,1])) #3.86 Angles_phy_HL_obs[2,1] <- angle(abs(Ont_Phy_coef_boot_HL[2,1:7,1]),abs(Ont_Phy_coef_boot_HL[4,1:7,1])) #2.64 Angles_phy_HL_obs[5,1] <- angle(abs(Ont_Phy_coef_boot_HL[5,1:7,1]),abs(Ont_Phy_coef_boot_HL[7,1:7,1])) #2.97 Angles_phy_HL_obs[3,1] <- angle(abs(Ont_Phy_coef_boot_HL[3,1:7,1]),abs(Ont_Phy_coef_boot_HL[4,1:7,1])) #3.49 Angles_phy_HL_obs[6,1] <- angle(abs(Ont_Phy_coef_boot_HL[6,1:7,1]),abs(Ont_Phy_coef_boot_HL[7,1:7,1])) #3.35 Angles_phy_HL_obs ``` ### Create dataset with shared allometry (rotate original datasets) ```{r, echo=F, Warning=F} #run PCA on all phylo together #fuse datasets Data_phylo_all <- smartbind(Adult_data_Ant,Adult_data_MLpri,Adult_data_MLsec) Data_phylo_all$Region_Det <- as.factor(gsub("Lesser Antilles Northern","Greater Antilles", Data_phylo_all$Region_Det)) Data_phylo_all$Region_Det <- droplevels(Data_phylo_all$Region_Det) rownames(Data_phylo_all) <- Data_phylo_all$species tree_all <- drop.tip(tree, name.check(tree, Data_phylo_all)$tree_not_data) #all common PCs res.cPCA_all_FL <- FCPCA(Data_phylo_all[,Forel], Data_phylo_all$Region_Det, graph=F) res.cPCA_all_HL <- FCPCA(Data_phylo_all[,Hindl], Data_phylo_all$Region_Det, graph=F) Adult_data_Ant_new <- Adult_data_Ant_new2 <- Adult_data_Ant[-c(1:nrow(Adult_data_Ant)),] Adult_data_MLpri_new <- Adult_data_MLpri_new2 <- Adult_data_MLpri[-c(1:nrow(Adult_data_MLpri)),] Adult_data_MLsec_new <- Adult_data_MLsec_new2 <- Adult_data_MLsec[-c(1:nrow(Adult_data_MLsec)),] ##Ant #Create dataset with shared allometry for phylo: X <- Adult_data_Ant FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_Ant_FL$a,nrow(Adult_data_Ant)),nrow(Adult_data_Ant),6,byrow=TRUE)) %*% res.pPCA_Ant_FL$Evec %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_Ant_HL$a,nrow(Adult_data_Ant)),nrow(Adult_data_Ant),7,byrow=TRUE)) %*% res.pPCA_Ant_HL$Evec %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_Ant_new <- rbind(Adult_data_Ant_new,X) #Create dataset with shared allometry for phylo without pPC: X <- Adult_data_Ant FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_Ant_FL$rotation %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_Ant_HL$rotation %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_Ant_new2 <- rbind(Adult_data_Ant_new2,X) ##MLpri #Create dataset with shared allometry for phylo: X <- Adult_data_MLpri FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_MLpri_FL$a,nrow(Adult_data_MLpri)),nrow(Adult_data_MLpri),6,byrow=TRUE)) %*% res.pPCA_MLpri_FL$Evec %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_MLpri_HL$a,nrow(Adult_data_MLpri)),nrow(Adult_data_MLpri),7,byrow=TRUE)) %*% res.pPCA_MLpri_HL$Evec %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLpri_new <- rbind(Adult_data_MLpri_new,X) #Create dataset with shared allometry for phylo without pPC: X <- Adult_data_MLpri FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_MLpri_FL$rotation %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_MLpri_HL$rotation %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLpri_new2 <- rbind(Adult_data_MLpri_new2,X) ##MLsec #Create dataset with shared allometry for phylo: X <- Adult_data_MLsec FL <- X[,Forel] FL_new <- as.matrix(FL-matrix(rep(res.pPCA_MLsec_FL$a,nrow(Adult_data_MLsec)),nrow(Adult_data_MLsec),6,byrow=TRUE)) %*% res.pPCA_MLsec_FL$Evec %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- as.matrix(HL-matrix(rep(res.pPCA_MLsec_HL$a,nrow(Adult_data_MLsec)),nrow(Adult_data_MLsec),7,byrow=TRUE)) %*% res.pPCA_MLsec_HL$Evec %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLsec_new <- rbind(Adult_data_MLsec_new,X) #Create dataset with shared allometry for phylo without pPC: X <- Adult_data_MLsec FL <- X[,Forel] FL_new <- scale(as.matrix(FL), center=T, scale=F) %*% res.PCA_MLsec_FL$rotation %*% t(res.cPCA_all_FL$loadings.common) HL <- X[,Hindl] HL_new <- scale(as.matrix(HL), center=T, scale=F) %*% res.PCA_MLsec_HL$rotation %*% t(res.cPCA_all_HL$loadings.common) X[,Forel] <- FL_new X[,Hindl] <- HL_new Adult_data_MLsec_new2 <- rbind(Adult_data_MLsec_new2,X) ``` ### Bootstrapping from datasets with shared multivariate allometry ```{r, echo=F} Phy_coef_sim_FL <- array(NA, dim=c(6,n.Forel,n.sim)) Phy_coef_sim_HL <- array(NA, dim=c(6,n.Hindl,n.sim)) rownames(Phy_coef_sim_FL) <- rownames(Phy_coef_sim_HL) <- c("Ant (pPCA)","MLpri (pPCA)","MLsec (pPCA)","Ant (standard PCA)","MLpri (standard PCA)","MLsec (standard PCA)") Phy_coef_sim_FL[1,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant_new[tree_Ant$tip.label, Forel]), tree_Ant, lambda.method = "ML", b = n.sim)$PC1 Phy_coef_sim_FL[2,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri_new[tree_MLpri$tip.label, Forel]), tree_MLpri, lambda.method = "ML", b = n.sim)$PC1 Phy_coef_sim_FL[3,1:6,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec_new[tree_MLsec$tip.label, Forel]), tree_MLsec, lambda.method = "ML", b = n.sim)$PC1 Phy_coef_sim_HL[1,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_Ant_new[tree_Ant$tip.label, Hindl]), tree_Ant, lambda.method = "ML", b = n.sim)$PC1 Phy_coef_sim_HL[2,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLpri_new[tree_MLpri$tip.label, Hindl]), tree_MLpri, lambda.method = "ML", b = n.sim)$PC1 Phy_coef_sim_HL[3,1:7,1:n.sim] <- sim.ecov(as.matrix(Adult_data_MLsec_new[tree_MLsec$tip.label, Hindl]), tree_MLsec, lambda.method = "ML", b = n.sim)$PC1 for(k in 1:n.sim){ Adult_data_Ant_new2_sub <- sample_frac(Adult_data_Ant_new2, 1, replace=T) Adult_data_MLpri_new2_sub <- sample_frac(Adult_data_MLpri_new2, 1, replace=T) Adult_data_MLsec_new2_sub <- sample_frac(Adult_data_MLsec_new2, 1, replace=T) Phy_coef_sim_FL[4,,k] <- prcomp(Adult_data_Ant_new2_sub[,Forel])$rotation[,1] Phy_coef_sim_HL[4,,k] <- prcomp(Adult_data_Ant_new2_sub[,Hindl])$rotation[,1] Phy_coef_sim_FL[5,,k] <- prcomp(Adult_data_MLpri_new2_sub[,Forel])$rotation[,1] Phy_coef_sim_HL[5,,k] <- prcomp(Adult_data_MLpri_new2_sub[,Hindl])$rotation[,1] Phy_coef_sim_FL[6,,k] <- prcomp(Adult_data_MLsec_new2_sub[,Forel])$rotation[,1] Phy_coef_sim_HL[6,,k] <- prcomp(Adult_data_MLsec_new2_sub[,Hindl])$rotation[,1] } Angles_phy_FL_sim <- Angles_phy_HL_sim <- array(NA, dim=c(6,1,n.sim)) rownames(Angles_phy_FL_sim) <- rownames(Angles_phy_HL_sim) <- c("Ant - MLpri (pPCA)","Ant - MLsec (pPCA)","MLpri - MLsec (pPCA)","Ant - MLpri (standard PCA)","Ant - MLsec (standard PCA)","MLpri - MLsec (standard PCA)") for(k in 1:n.sim){ Angles_phy_FL_sim[1,,k] <- angle(abs(Phy_coef_sim_FL[1,,k]),abs(Phy_coef_sim_FL[2,,k])) Angles_phy_HL_sim[1,,k] <- angle(abs(Phy_coef_sim_HL[1,,k]),abs(Phy_coef_sim_HL[2,,k])) Angles_phy_FL_sim[2,,k] <- angle(abs(Phy_coef_sim_FL[1,,k]),abs(Phy_coef_sim_FL[3,,k])) Angles_phy_HL_sim[2,,k] <- angle(abs(Phy_coef_sim_HL[1,,k]),abs(Phy_coef_sim_HL[3,,k])) Angles_phy_FL_sim[3,,k] <- angle(abs(Phy_coef_sim_FL[2,,k]),abs(Phy_coef_sim_FL[3,,k])) Angles_phy_HL_sim[3,,k] <- angle(abs(Phy_coef_sim_HL[2,,k]),abs(Phy_coef_sim_HL[3,,k])) Angles_phy_FL_sim[4,,k] <- angle(abs(Phy_coef_sim_FL[4,,k]),abs(Phy_coef_sim_FL[5,,k])) Angles_phy_HL_sim[4,,k] <- angle(abs(Phy_coef_sim_HL[4,,k]),abs(Phy_coef_sim_HL[5,,k])) Angles_phy_FL_sim[5,,k] <- angle(abs(Phy_coef_sim_FL[4,,k]),abs(Phy_coef_sim_FL[6,,k])) Angles_phy_HL_sim[5,,k] <- angle(abs(Phy_coef_sim_HL[4,,k]),abs(Phy_coef_sim_HL[6,,k])) Angles_phy_FL_sim[6,,k] <- angle(abs(Phy_coef_sim_FL[5,,k]),abs(Phy_coef_sim_FL[6,,k])) Angles_phy_HL_sim[6,,k] <- angle(abs(Phy_coef_sim_HL[5,,k]),abs(Phy_coef_sim_HL[6,,k])) } Angles_FL_sim_95CI <- apply(Angles_phy_FL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) Angles_HL_sim_95CI <- apply(Angles_phy_HL_sim, c(1,2), quantile, na.rm=TRUE, probs=0.95) paste("95% CI of simulated (parallel) angles between evolutionary axis - forelimb") Angles_FL_sim_95CI paste("95% CI of simulated (parallel) angles between evolutionary axis - hindlimb") Angles_HL_sim_95CI #get P-values P_theta_phy_FL <- P_theta_phy_HL <- array(NA, dim=c(6,1)) rownames(P_theta_phy_FL) <- rownames(P_theta_phy_HL) <- c("Ant - MLpri (pPCA)","Ant - MLsec (pPCA)","MLpri - MLsec (pPCA)","Ant - MLpri (standard PCA)","Ant - MLsec (standard PCA)","MLpri - MLsec (standard PCA)") for(i in 1:6){ P_theta_phy_FL[i,1] <- 1-sum(as.vector(na.omit(Angles_phy_FL_sim[i,1,1:n.sim])) < Angles_phy_FL_obs[i,1])/length(as.vector(na.omit(Angles_phy_FL_sim[i,1,1:n.sim]))) P_theta_phy_HL[i,1] <- 1-sum(as.vector(na.omit(Angles_phy_HL_sim[i,1,1:n.sim])) < Angles_phy_HL_obs[i,1])/length(as.vector(na.omit(Angles_phy_HL_sim[i,1,1:n.sim]))) } print("P-values and angles theta - FL") P_theta_phy_FL print("P-values and angles theta - HL") P_theta_phy_HL #save.image(file = "C:/Users/Nathalie/Dropbox/stuff/MSc_Projects/ElianeVanDerCruyssen_2020/OntoPhylo_Final.RData") ```