library(RRPP) library(phytools) library(geiger) library(phylolm) library(ggplot2) library(forcats) data <- read.csv("data/data_hbER_220627.csv") rownames(data) <- data[,1] tree <- read.nexus("data/squirrel_mcctree.nex") #remove species in phylogeny that aren't present in dataset tree_pr <- treedata(phy = tree, data = data, sort=TRUE)$phy #prune data to check if all species are in tree #data_mean_pr <- data.frame(treedata(phy = tree_pr, data = data, sort=TRUE)$data) size <- (as.numeric(as.character(data$geomean))) #save size as vector names(size) <- rownames(data) #give species names for size vector lnsize <- log(size) #natural log transform size vector (for stats purposes) hbER <- (as.numeric(as.character(data$hbER))) #save hbER as vector names(hbER) <- rownames(data) #give species names for hbER vector lnhbER <- log(hbER) #natural log transform hbER vector (for stats purposes) ecology4 <- data$ecotype4 #save ecology as factor names(ecology4) <- rownames(data) #give species names for ecology factor col_ecology4 <- c("#073b4c", "#118ab2", "#ffd166", "#06d6a0") #set colors for ecological groups: blue = flying; orange = ground; green = tree headER <- (as.numeric(as.character(data$headER))) names(headER) <- rownames(data) lnheadER <- log(headER) AEI_C <- (as.numeric(as.character(data$AEI_C))) names(AEI_C) <- rownames(data) lnAEI_C <- log(AEI_C) AEI_T <- (as.numeric(as.character(data$AEI_T))) names(AEI_T) <- rownames(data) lnAEI_T <- log(AEI_T) AEI_L <- (as.numeric(as.character(data$AEI_L))) names(AEI_L) <- rownames(data) lnAEI_L <- log(AEI_L) AEI_S <- (as.numeric(as.character(data$AEI_S))) names(AEI_S) <- rownames(data) lnAEI_S <- log(AEI_S) rib <- (as.numeric(as.character(data$avgRibLength))) names(rib) <- rownames(data) #give species names for hbER vector lnrib <- log(rib) lnrib_sc_raw <- (phylolm(lnrib~lnsize, phy = tree_pr, model = "lambda"))$residuals lnrib_sc <- lnrib_sc_raw[match(names(lnrib),names(lnrib_sc_raw))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!! # plots of lnhbER and each component plot(lnhbER ~ lnheadER, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) plot(lnhbER ~ lnAEI_C, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) plot(lnhbER ~ lnAEI_T, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) plot(lnhbER ~ lnAEI_L, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) plot(lnhbER ~ lnAEI_S, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) plot(lnhbER ~ lnrib_sc, pch = c(24, 21, 22, 23)[factor(ecology4)], bg = col_ecology4[factor(ecology4)], las =1, cex = 1.75, cex.axis=1.5) #### Mulitple Regression - ER components #### ERcomp_rrpp_mccphy <- lm.rrpp(lnhbER ~ lnheadER + lnAEI_C + lnAEI_T + lnAEI_L + lnAEI_S + lnrib_sc, Cov=vcv.phylo(tree_pr), iter = 999) anova(ERcomp_rrpp_mccphy) #### chipmunk Mulitple Regression - ER components #### data_chipmunk <- filter(data, ecotype4 == "chipmunk") rownames(data_chipmunk) <- data_chipmunk[,1] tree_pr_chipmunk <- treedata(phy = tree_pr, data = data_chipmunk)$phy lnhbER_chipmunk <- log(as.numeric(as.character(data_chipmunk$hbER))) names(lnhbER_chipmunk) <- rownames(data_chipmunk) lnsize_chipmunk <- log(as.numeric(as.character(data_chipmunk$geomean))) names(lnsize_chipmunk) <- rownames(data_chipmunk) lnheadER_chipmunk <- log(as.numeric(as.character(data_chipmunk$headER))) names(lnheadER_chipmunk) <- rownames(data_chipmunk) lnAEI_C_chipmunk <- log(as.numeric(as.character(data_chipmunk$AEI_C))) names(lnAEI_C_chipmunk) <- rownames(data_chipmunk) lnAEI_T_chipmunk <- log(as.numeric(as.character(data_chipmunk$AEI_T))) names(lnAEI_T_chipmunk) <- rownames(data_chipmunk) lnAEI_L_chipmunk <- log(as.numeric(as.character(data_chipmunk$AEI_L))) names(lnAEI_L_chipmunk) <- rownames(data_chipmunk) lnAEI_S_chipmunk <- log(as.numeric(as.character(data_chipmunk$AEI_S))) names(lnAEI_S_chipmunk) <- rownames(data_chipmunk) lnrib_chipmunk <- log(as.numeric(as.character(data_chipmunk$avgRibLength))) names(lnrib_chipmunk) <- rownames(data_chipmunk) lnrib_sc_raw_chipmunk <- (phylolm(lnrib_chipmunk~lnsize_chipmunk, phy = tree_pr_chipmunk, model = "lambda"))$residuals lnrib_sc_chipmunk <- lnrib_sc_raw_chipmunk[match(names(lnrib_chipmunk),names(lnrib_sc_raw_chipmunk))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!! plot(lnrib_sc_raw_chipmunk~lnsize_chipmunk) #check out plots to see how they can lead to different results if they are in different order!! plot(lnrib_sc_chipmunk~lnsize_chipmunk) plot(lnhbER_chipmunk ~ lnheadER_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_chipmunk ~ lnAEI_C_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_chipmunk ~ lnAEI_T_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_chipmunk ~ lnAEI_L_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_chipmunk ~ lnAEI_S_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_chipmunk ~ lnrib_sc_chipmunk, pch = c( 22), las =1, cex = 1.75, cex.axis=2) ERcomp_rrpp_mccphy_chipmunk <- lm.rrpp(lnhbER_chipmunk ~ lnheadER_chipmunk + lnAEI_C_chipmunk + lnAEI_T_chipmunk + lnAEI_L_chipmunk + lnAEI_S_chipmunk + lnrib_sc_chipmunk, Cov=vcv.phylo(tree_pr_chipmunk), iter = 999) anova(ERcomp_rrpp_mccphy_chipmunk) #### flying Mulitple Regression - ER components #### data_flying <- filter(data, ecotype4 == "flying") rownames(data_flying) <- data_flying[,1] tree_pr_flying <- treedata(phy = tree_pr, data = data_flying)$phy lnhbER_flying <- log(as.numeric(as.character(data_flying$hbER))) names(lnhbER_flying) <- rownames(data_flying) lnsize_flying <- log(as.numeric(as.character(data_flying$geomean))) names(lnsize_flying) <- rownames(data_flying) lnheadER_flying <- log(as.numeric(as.character(data_flying$headER))) names(lnheadER_flying) <- rownames(data_flying) lnAEI_C_flying <- log(as.numeric(as.character(data_flying$AEI_C))) names(lnAEI_C_flying) <- rownames(data_flying) lnAEI_T_flying <- log(as.numeric(as.character(data_flying$AEI_T))) names(lnAEI_T_flying) <- rownames(data_flying) lnAEI_L_flying <- log(as.numeric(as.character(data_flying$AEI_L))) names(lnAEI_L_flying) <- rownames(data_flying) lnAEI_S_flying <- log(as.numeric(as.character(data_flying$AEI_S))) names(lnAEI_S_flying) <- rownames(data_flying) lnrib_flying <- log(as.numeric(as.character(data_flying$avgRibLength))) names(lnrib_flying) <- rownames(data_flying) lnrib_sc_raw_flying <- (phylolm(lnrib_flying~lnsize_flying, phy = tree_pr_flying, model = "lambda"))$residuals lnrib_sc_flying <- lnrib_sc_raw_flying[match(names(lnrib_flying),names(lnrib_sc_raw_flying))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!! plot(lnrib_sc_raw_flying~lnsize_flying) #check out plots to see how they can lead to different results if they are in different order!! plot(lnrib_sc_flying~lnsize_flying) plot(lnhbER_flying ~ lnheadER_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_flying ~ lnAEI_C_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_flying ~ lnAEI_T_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_flying ~ lnAEI_L_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_flying ~ lnAEI_S_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_flying ~ lnrib_sc_flying, pch = c( 22), las =1, cex = 1.75, cex.axis=2) ERcomp_rrpp_mccphy_flying <- lm.rrpp(lnhbER_flying ~ lnheadER_flying + lnAEI_C_flying + lnAEI_T_flying + lnAEI_L_flying + lnAEI_S_flying + lnrib_sc_flying, Cov=vcv.phylo(tree_pr_flying), iter = 999) anova(ERcomp_rrpp_mccphy_flying) #### Ground Mulitple Regression - ER components #### data_ground <- filter(data, ecotype4 == "ground") rownames(data_ground) <- data_ground[,1] tree_pr_ground <- treedata(phy = tree_pr, data = data_ground)$phy lnhbER_ground <- log(as.numeric(as.character(data_ground$hbER))) names(lnhbER_ground) <- rownames(data_ground) lnsize_ground <- log(as.numeric(as.character(data_ground$geomean))) names(lnsize_ground) <- rownames(data_ground) lnheadER_ground <- log(as.numeric(as.character(data_ground$headER))) names(lnheadER_ground) <- rownames(data_ground) lnAEI_C_ground <- log(as.numeric(as.character(data_ground$AEI_C))) names(lnAEI_C_ground) <- rownames(data_ground) lnAEI_T_ground <- log(as.numeric(as.character(data_ground$AEI_T))) names(lnAEI_T_ground) <- rownames(data_ground) lnAEI_L_ground <- log(as.numeric(as.character(data_ground$AEI_L))) names(lnAEI_L_ground) <- rownames(data_ground) lnAEI_S_ground <- log(as.numeric(as.character(data_ground$AEI_S))) names(lnAEI_S_ground) <- rownames(data_ground) lnrib_ground <- log(as.numeric(as.character(data_ground$avgRibLength))) names(lnrib_ground) <- rownames(data_ground) lnrib_sc_raw_ground <- (phylolm(lnrib_ground~lnsize_ground, phy = tree_pr_ground, model = "lambda"))$residuals lnrib_sc_ground <- lnrib_sc_raw_ground[match(names(lnrib_ground),names(lnrib_sc_raw_ground))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!! plot(lnrib_sc_raw_ground~lnsize_ground) #check out plots to see how they can lead to different results if they are in different order!! plot(lnrib_sc_ground~lnsize_ground) plot(lnhbER_ground ~ lnheadER_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_ground ~ lnAEI_C_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_ground ~ lnAEI_T_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_ground ~ lnAEI_L_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_ground ~ lnAEI_S_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_ground ~ lnrib_sc_ground, pch = c( 22), las =1, cex = 1.75, cex.axis=2) ERcomp_rrpp_mccphy_ground <- lm.rrpp(lnhbER_ground ~ lnheadER_ground + lnAEI_C_ground + lnAEI_T_ground + lnAEI_L_ground + lnAEI_S_ground + lnrib_sc_ground, Cov=vcv.phylo(tree_pr_ground), iter = 999) anova(ERcomp_rrpp_mccphy_ground) #### tree Mulitple Regression - ER components #### data_tree <- filter(data, ecotype4 == "tree") rownames(data_tree) <- data_tree[,1] tree_pr_tree <- treedata(phy = tree_pr, data = data_tree)$phy lnhbER_tree <- log(as.numeric(as.character(data_tree$hbER))) names(lnhbER_tree) <- rownames(data_tree) lnsize_tree <- log(as.numeric(as.character(data_tree$geomean))) names(lnsize_tree) <- rownames(data_tree) lnheadER_tree <- log(as.numeric(as.character(data_tree$headER))) names(lnheadER_tree) <- rownames(data_tree) lnAEI_C_tree <- log(as.numeric(as.character(data_tree$AEI_C))) names(lnAEI_C_tree) <- rownames(data_tree) lnAEI_T_tree <- log(as.numeric(as.character(data_tree$AEI_T))) names(lnAEI_T_tree) <- rownames(data_tree) lnAEI_L_tree <- log(as.numeric(as.character(data_tree$AEI_L))) names(lnAEI_L_tree) <- rownames(data_tree) lnAEI_S_tree <- log(as.numeric(as.character(data_tree$AEI_S))) names(lnAEI_S_tree) <- rownames(data_tree) lnrib_tree <- log(as.numeric(as.character(data_tree$avgRibLength))) names(lnrib_tree) <- rownames(data_tree) lnrib_sc_raw_tree <- (phylolm(lnrib_tree~lnsize_tree, phy = tree_pr_tree, model = "lambda"))$residuals lnrib_sc_tree <- lnrib_sc_raw_tree[match(names(lnrib_tree),names(lnrib_sc_raw_tree))] ### very important - must make sure lnrib_sc is in the same order as your other variables!!! plot(lnrib_sc_raw_tree~lnsize_tree) #check out plots to see how they can lead to different results if they are in different order!! plot(lnrib_sc_tree~lnsize_tree) plot(lnhbER_tree ~ lnheadER_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_tree ~ lnAEI_C_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_tree ~ lnAEI_T_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_tree ~ lnAEI_L_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_tree ~ lnAEI_S_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) plot(lnhbER_tree ~ lnrib_sc_tree, pch = c( 22), las =1, cex = 1.75, cex.axis=2) ERcomp_rrpp_mccphy_tree <- lm.rrpp(lnhbER_tree ~ lnheadER_tree + lnAEI_C_tree + lnAEI_T_tree + lnAEI_L_tree + lnAEI_S_tree + lnrib_sc_tree, Cov=vcv.phylo(tree_pr_tree), iter = 999) anova(ERcomp_rrpp_mccphy_tree)