#clear the memory rm(list=ls()) #If you have never installed these packages remove the # to run the install.packages code #install.packages("ape") #install.packages("geiger") #install.packages("phytools") #install.packages("dplyr") #install.packages("caper") #Now we can turn the packages on library(ape) library(geiger) library(phytools) library(dplyr) library(caper) #library(evomap) library(ggplot2) #3)Reading in Data and Phylogeny------------------------------------ #this gets rid of factors. options(stringsAsFactors = FALSE) #reading in data data <- read.csv("Carnivore_data_corrected.csv") #take out the otter locomotionnootter<-data[which(data$group!="otter"),] data<-locomotionnootter #take out additional locomotionnoweasel<-data[which(data$group!="kusimans"),] data<-locomotionnoweasel specname<-data$Taxon #reading in phylogenetic tree Fixed_tree<-read.nexus("Fritztree.rs200k.100trees.tre") #class(Fixed_tree)<-"phylo" #prune tree to these taxa species<-c(specname) rownames(data)<-species #Creating space in which to save the results rk_signal<- rk_p_low<- rk_p_upp<- ml_list<- CI_upper <- CI_lower <- intercept_list<- int.tstat_list<- int.pval_list<- mass.slope_list<- mass.tstat_list<- mass.pval_list<- Aquatic.intercept_list<- AquaticInt.tstat_list<- AquaticInt.pval_list<- Terrestrial.intercept_list<- TerrestrialInt.tstat_list<- TerrestrialInt.pval_list<- Volant.intercept_list<- VolantInt.tstat_list<- VolantInt.pval_list<- Fossorial.intercept_list<- FossorialInt.tstat_list<- FossorialInt.pval_list<-Aquatic.slope_list<- AquaticSlope.tstat_list<- AquaticSlope.pval_list<- Terrestrial.slope_list<- TerrestrialSlope.tstat_list<- TerrestrialSlope.pval_list<- Volant.slope_list<- VolantSlope.tstat_list<- VolantSlope.pval_list<- Fossorial.slope_list<- FossorialSlope.tstat_list<- FossorialSlope.pval_list<-multiRsq_list<- multiRsq_list <- multiadjRsq_list <- PanCov.res1 <- PanCov.res2 <- matrix(0,length(Fixed_tree)) #the loop starts here. for(i in 1: length(Fixed_tree)) { tryCatch({ pruned.trees<-drop.tip(Fixed_tree[[i]], setdiff(Fixed_tree[[i]]$tip.label, species)) class(pruned.trees)<-"phylo" #match trees to the data tree<-treedata(pruned.trees,data,sort=T,warnings=T)$phy #4)Getting data and tree ready------------------------------------- #Let's first specify row names to allow matching of tree and data rownames(data) <- data$Taxon #5)Plotting the tree!------------------------------------------------ plot.phylo(pruned.trees, type="phylogram", show.tip=TRUE, main="carnivore tree") #If it seems smushed you can always click on teh zoom button to see the plot larger #Discrete character maps #Discrete character maps allow us to visualize the occurrence of discrete traits, like the presence/absence of a certain phenotypes (like color) or behavior (e.g., burrowing, paternal care, diet ect.), on a phylogenetic tree. This visualization technique is helpful when trying to see how traits are spread amoung your groups of animals. #We need to make a vector of our character data named with the corresponding species names. Locomotion<-data$Locomotion #Next we will plot the tree with a title and a legend plot(pruned.trees, type="phylogram", cex=0.9, lwd=2, label.offset=4) title("carnivore Locomotion") legend("topright", legend= sort(unique(Locomotion)), col=c("green","blue","purple", "yellow"), pch=19, pt.cex = 3) #Continuous character maps This maps continuous traits, (e.g., body size or M1 #length) on the phylogenetic tree. With these graphics we are able to visualize #how the dimensions of body size and or M1 length are distributed among species. #making named vector logbodymass<-data$lnmass names(logbodymass) <- rownames(data) #same code but for PC1..RSI PC1_RSI <-data$PC1 names(PC1_RSI) <- rownames(data) #what do you think of these trees that you have plotted? Are there lineages that #are particularly small? Does body mass and M1L track each other? Can you tell #by looking at the trees? #6)Incorporating Phylogenetic signal!--------- ok so what if after looking at #our trees we wanted to know if body mass and M1 length have a relationship but #we are concerned that the relationships amoung the animals might be a problem. #We can correct our regression to account for this information (you can do this #for more than just regressions, I have used this method to correct my data when #I was testing for differences in body mass related to diet which requires tests #like ANOVA's) #First let's perform a non-phylogenetic regression to see what happens. lm=linear model model1<-lm(PC1_RSI~logbodymass, data=data) summary(model1) #plotting the regression plot(PC1_RSI~logbodymass, data=data) #now add the line abline(model1, col="green") #What do you think? What is your regression like? #Now we can use a Phylogenetic generalized least squares (PGLS) analysis. It can #fit a linear model, taking into account phylogenetic non-independence between #data points. #We will use our caper object that combines phylogeny and data: carnivore<- comparative.data(pruned.trees, data, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) # Modify code to grab the appropriate subset of data? #We will now fit a model using the function pgls: #have to let lambda be the same between the two, so cannot ML fit it; using the default of lambda = 1 carnivore.pgls<-pgls(PC1~lnmass, data=carnivore) summary(carnivore.pgls) cpgls<-summary(carnivore.pgls) #PGLS for diet categories categories.pgls <- pgls(PC1~lnmass * Locomotion, data=carnivore) summary(categories.pgls) categories1.pgls <- pgls(PC1~lnmass + Locomotion, data=carnivore) summary(categories1.pgls) #this anova compares the overall model to the model with the diet categories fit separately pancova.result <- anova(carnivore.pgls, categories.pgls, categories1.pgls) pancova.result #need to assign to variable for reporting out from loop: PanCov.res1[i] <- pancova.result$`Pr(>F)`[2] PanCov.res2[i] <- pancova.result$`Pr(>F)`[3] #setting our plotting window for two graphs par(mfrow=c(1,2)) #plot PGLS model plot(PC1_RSI~logbodymass, data=data) abline(carnivore.pgls, col="green") title("PGLS") #Plot original uninformed regression plot(PC1_RSI~logbodymass, data=data) abline(model1, col="green") title("Original regression") #Incorporate code to "pick out" carnivores based on diet category and plot them against the general trend? Terrestrial<-subset(data,Locomotion=="Terrestrial",select = c("Taxon","lnmass","PC1")) Aquatic<-subset(data,Locomotion=="Aquatic",select = c("Taxon","lnmass","PC1")) Volant<-subset(data,Locomotion=="Volant",select = c("Taxon","lnmass","PC1")) Fossorial<-subset(data,Locomotion=="Fossorial",select = c("Taxon","lnmass","PC1")) # generate pgls for 3 diet subcategories Terrestrial_Carnivores<- comparative.data(pruned.trees, Terrestrial, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) Aquatic_Carnivores<- comparative.data(pruned.trees, Aquatic, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) Volant_Carnivores<- comparative.data(pruned.trees, Volant, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) Fossorial_Carnivores<- comparative.data(pruned.trees, Fossorial, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) #-----?? data dropped? #Terrestrial pgls Terrestrial.pgls<-pgls(PC1~lnmass, data=Terrestrial_Carnivores, lambda="ML") summary(Terrestrial.pgls) Tpgls<-summary(Terrestrial.pgls) #Aquatic pgls Aquatic.pgls<-pgls(PC1~lnmass, data=Aquatic_Carnivores, lambda="ML") summary(Aquatic.pgls) Spgls<-summary(Aquatic.pgls) #Volant pgls Volant.pgls<-pgls(PC1~lnmass, data=Volant_Carnivores, lambda="ML") summary(Volant.pgls) Vpgls<-summary(Volant.pgls) #Fossorial pgls Fossorial.pgls<-pgls(PC1~lnmass, data=Fossorial_Carnivores, lambda="ML") summary(Fossorial.pgls) Fpgls<-summary(Fossorial.pgls) #create models for the 3 diet categories model1<-lm(PC1~lnmass, data=Fossorial) summary(model1) model2<-lm(PC1~lnmass, data=Terrestrial) summary(model2) model4<-lm(PC1~lnmass, data=Aquatic) summary(model4) model3<-lm(PC1~lnmass, data=Volant) summary(model3) #Plot original uninformed/pgls regression comparissons for 3 diet subcategories #Terrestrial Carnivores plot(PC1~lnmass, data=Terrestrial) abline(model2, col="red") title("Original regression Terrestrial Carnivores") #plot PGLS model(Terrestrial_Carnivores) plot(PC1~lnmass, data=Terrestrial) abline(Terrestrial.pgls, col="red") title("PGLS_Terrestrial_Carnivores") #Swimmer Feeders plot(PC1~lnmass, data=Aquatic) abline(model4, col="blue") title("Original regression Aquatic Carnivores") #plot PGLS model(Aquatic_ Carnivores) plot(PC1~lnmass, data=Aquatic) abline(Aquatic.pgls, col="blue") title("PGLS_Aquatic_Carnivores") #Volant Feeders plot(PC1~lnmass, data=Volant) abline(model2, col="pink") title("Original regression Volant_Carnivores") #plot PGLS model(InV_feeders) plot(PC1~lnmass, data=Volant) abline(Volant.pgls, col="pink") title("PGLS_Volant_Carnivores") # Fossorial Feeders plot(PC1~lnmass, data=Fossorial) abline(model2, col="grey") title("Original regression Fossorial_Carnivores") #plot PGLS model(InV_feeders) plot(PC1~lnmass, data=Fossorial) abline(Volant.pgls, col="grey") title("PGLS_Fossorial_Carnivores") #plot all PGLS on same axis plot(PC1~lnmass,data=data) abline(carnivore.pgls,col="white")+ abline(Terrestrial.pgls,col="red")+ abline(Aquatic.pgls,col="blue")+ abline(Volant.pgls,col="pink")+ abline(Fossorial.pgls,col="grey")+ title("pgls_carnivore_regression") legend("topright", legend= sort(unique(Locomotion)), col=c("blue","green","grey","red","pink"), pch=19, pt.cex = 3) #plot uninformed regression for all groups on same axis plot(PC1~lnmass,data=data) abline(model1,col="white")+ abline(model2,col="red")+ abline(model4,col="blue")+ abline(model3,col="pink")+ abline(model1,col="grey")+ title("uninformed_carnivore_regression") legend("topright", legend= sort(unique(Locomotion)), col=c("blue","green","grey","red","pink"), pch=19, pt.cex = 3) #Assign data locations mass.slope_list[i]<-cpgls$coefficients[2,1] mass.tstat_list[i]<-cpgls$coefficients[2,3] mass.pval_list[i]<-cpgls$coefficients[2,4] intercept_list[i] <- cpgls$coefficients[1,1] int.tstat_list[i]<-cpgls$coefficients[1,3] int.pval_list[i]<-cpgls$coefficients[1, 4] Aquatic.intercept_list[i]<-Spgls$coefficients[1,1] AquaticInt.tstat_list[i]<-Spgls$coefficients[1,3] AquaticInt.pval_list[i]<-Spgls$coefficients[1, 4] Terrestrial.intercept_list[i]<- Tpgls$coefficients[1,1] TerrestrialInt.tstat_list[i]<-Tpgls$coefficients[1,3] TerrestrialInt.pval_list[i]<-Tpgls$coefficients[1, 4] Volant.intercept_list[i]<- Vpgls$coefficients[1,1] VolantInt.tstat_list[i]<-Vpgls$coefficients[1,3] VolantInt.pval_list[i]<-Vpgls$coefficients[1, 4] Fossorial.intercept_list[i]<- Fpgls$coefficients[1,1] FossorialInt.tstat_list[i]<-Fpgls$coefficients[1,3] FossorialInt.pval_list[i]<-Fpgls$coefficients[1, 4] #update the other slopes like this vert slope Aquatic.slope_list[i]<-Spgls$coefficients[2,1] AquaticSlope.tstat_list[i]<-Spgls$coefficients[2,3] AquaticSlope.pval_list[i]<-Spgls$coefficients[2,4] Volant.slope_list[i]<-Vpgls$coefficients[2,1] VolantSlope.tstat_list[i]<-Vpgls$coefficients[2,3] VolantSlope.pval_list[i]<-Vpgls$coefficients[2,4] Fossorial.slope_list[i]<-Fpgls$coefficients[2,1] FossorialSlope.tstat_list[i]<-Fpgls$coefficients[2,3] FossorialSlope.pval_list[i]<-Fpgls$coefficients[2,4] Terrestrial.slope_list[i]<- Tpgls$coefficients[2,1] TerrestrialSlope.tstat_list[i]<-Tpgls$coefficients[2,3] TerrestrialSlope.pval_list[i]<-Tpgls$coefficients[2, 4]}) } #end of loop # generate histograms of P values from pancova 1 hist(PanCov.res1,main="pancova1_results",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0,0.4),ylim=c(0,10),col="orange",breaks=100) Median_Panco1<-median(PanCov.res1) Mean_Panco1<-mean(PanCov.res1) #generate histogram of P values from pancova 2 hist(PanCov.res2,main="pancova2_results",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0,1),ylim=c(0,10),col="pink",breaks=100) Median_Panco2<-median(PanCov.res2) Mean_Panco2<-mean(PanCov.res2) # Plot histograms for Aquatic carnivore slope, intercept, and P-values for each hist(Aquatic.slope_list,main="Aquatic carnivore slopes",border="black",xlab="Slope_values ",ylab="number_of_occurences",xlim=c(-0.63,-0.43),ylim=c(0,6),col="blue",breaks=100) Median_Aquaslope<-median(Aquatic.slope_list) Mean_Aquaslope<-mean(Aquatic.slope_list) #Aquatic_slope P- values hist(AquaticSlope.pval_list,main="Aquatic_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.002,0.02),ylim=c(0,6),col="blue",breaks=100) Median_Aquaslopval<-median(AquaticSlope.pval_list) Mean_Aquaslopval<-mean(AquaticSlope.pval_list) # Now Aquatic carnivore intercepts hist(Aquatic.intercept_list,main="Aquatic carnivore Intercepts",border="black",xlab="Intercept_values ",ylab="number_of_occurences",xlim=c(2,5),ylim=c(0,6),col="blue",breaks=100) Median_Aqua_Intercept<-median(Aquatic.intercept_list) Mean_Aqua_Intercept<-mean(Aquatic.intercept_list) #Aquatic carnivore Intercept P-values hist(AquaticInt.pval_list,main="Aquatic carnivore intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.08,0.3),ylim=c(0,6),col="blue",breaks=100) Median_Aquainpval<-median(AquaticInt.pval_list) Mean_Aquainpval<-mean(AquaticInt.pval_list) # do the same for Terrestrial carnivores hist(Terrestrial.slope_list,main="Terrestrial carnivore slopes",border="black",xlab="Slope_values ",ylab="number_of_occurences",xlim=c(-0.223,-0.2227),ylim=c(0,6),col="red",breaks=100) Median_Terrestrialslope<-median(Terrestrial.slope_list) Mean_Terrestrialslope<-mean(Terrestrial.slope_list) #Terrestrial_slope P- values hist(TerrestrialSlope.pval_list,main="terrestrial_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(9.503287e-12,8.861623e-11),ylim=c(0,15),col="red",breaks=100) Median_Terrestrialslopval<-median(TerrestrialSlope.pval_list) Mean_Terrestrialslopval<-mean(TerrestrialSlope.pval_list) # Now Terrestrial carnivore intercepts hist(Terrestrial.intercept_list,main="Terrestrial carnivore Intercepts",border="black",xlab="values ",ylab="number_of_occurences",xlim=c(1.536,1.54),ylim=c(0,15),col="red",breaks=100) Median_Terrestrial_Intercept<-median(Terrestrial.intercept_list) Mean_Terrestrial_Intercept<-mean(Terrestrial.intercept_list) #Terrestrial carnivore Intercept P-values hist(TerrestrialInt.pval_list,main="Terrestrial_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.0008,0.0036),ylim=c(0,8),col="red",breaks=100) Median_Terrestrialintpval<-median(TerrestrialInt.pval_list) Mean_Terrestrialintpval<-mean(TerrestrialInt.pval_list) #same thing for Volant carnivores hist(Volant.slope_list,main="Volant carnivore slopes",border="black",xlab="values ",ylab="number_of_occurences",xlim=c(-0.13,-0.078),ylim=c(0,8),col="purple",breaks=100) Volant_Median_slope<-median(Volant.slope_list) Volant_Mean_slope<-mean(Volant.slope_list) # Volant_slope P- values hist(VolantSlope.pval_list,main="Volant_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.2,0.45),ylim=c(0,15),col="purple",breaks=100) Volant_Median_slopval<-median(VolantSlope.pval_list) Volant_Mean_slopval<-mean(VolantSlope.pval_list) # Now Volant carnivore intercepts hist(Volant.intercept_list,main="Volant carnivore Intercepts",border="black",xlab="values ",ylab="number_of_occurences",xlim=c(-0.0009,0.07),ylim=c(0,6),col="purple",breaks=100) Volant_Median_Intercept<-median(Volant.intercept_list) Volant_Mean_Intercept<-mean(Volant.intercept_list) #Volant carnivore Intercept P-values hist(VolantInt.pval_list,main="Volant_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.8,1.01),ylim=c(0,10),col="purple",breaks=100) Volant_Int_MedPval<-median(VolantInt.pval_list) Volant_INT_MeaPval<-mean(VolantInt.pval_list) # Same thing for Fossorial animals hist(Fossorial.slope_list,main="Fossorial carnivore slopes",border="black",xlab="values ",ylab="number_of_occurences",xlim=c(-0.152,-0.15),ylim=c(0,6),col="grey",breaks=100) Median_Fossorialslope<-median(Fossorial.slope_list) Mean_Fossorialslope<-mean(Fossorial.slope_list) # Fossorial_slope P- values hist(FossorialSlope.pval_list,main="Fossorial_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.0083,0.0084),ylim=c(0,6),col="grey",breaks=100) Median_Fossorialslopval<-median(FossorialSlope.pval_list) Mean_Fossorialslopval<-mean(FossorialSlope.pval_list) # Now Fossorial carnivore intercepts hist(Fossorial.intercept_list,main="Fossorial carnivore Intercepts",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.9,1),ylim=c(0,6),col="grey",breaks=100) Median_Fossorial_Intercept<-median(Fossorial.intercept_list) Mean_Fossorial_Intercept<-mean(Fossorial.intercept_list) #Fossorial carnivore Intercept P-values hist(FossorialInt.pval_list,main="Fossorial_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.0006,0.0007),ylim=c(0,6),col="grey",breaks=100) Median_Fossorialinpval<-Median_Fossorialinpval<-median(FossorialInt.pval_list) Mean_Fossorialinpval<-mean(FossorialInt.pval_list) #getting data for overall model median_carnslopval<-median(mass.pval_list) median_carninpval<-median(int.pval_list) median_carnslope<-median(mass.slope_list) median_carnin<-median(intercept_list) #making a df of the summary results SummaryResults <- data.frame(Median_Panco1, Mean_Panco1, Median_Panco2, Mean_Panco2, Median_Aquaslope, Mean_Aquaslope, Median_Aquaslopval, Mean_Aquaslopval, Median_Aqua_Intercept, Mean_Aqua_Intercept, Median_Aquainpval, Mean_Aquainpval, Median_Terrestrialslope, Mean_Terrestrialslope, Median_Terrestrialslopval, Mean_Terrestrialslopval, Median_Terrestrial_Intercept, Mean_Terrestrial_Intercept, Median_Terrestrialintpval, Mean_Terrestrialintpval, Volant_Median_slope, Volant_Mean_slope, Volant_Median_slopval, Volant_Mean_slopval, Volant_Median_Intercept, Volant_Mean_Intercept, Volant_Int_MedPval, Volant_INT_MeaPval, Median_Fossorialslope, Mean_Fossorialslope, Median_Fossorialslopval, Mean_Fossorialslopval, Median_Fossorial_Intercept, Mean_Fossorial_Intercept, Median_Fossorialinpval, Mean_Fossorialinpval) #writing out the df of summary results to allow skipping the loop write.csv(SummaryResults, file ="SummaryResults.csv", row.names = FALSE) #Start here if you skip the loop--------------------- #reading in the summary results if you have skipped the loop start here read_in <- read.csv("SummaryResults.csv") #just assigning all of the results back to their appropriate variables Median_Panco1 <- read_in$Median_Panco1 Mean_Panco1 <- read_in$Mean_Panco1 Median_Panco2 <- read_in$Median_Panco2 Mean_Panco2 <- read_in$Mean_Panco2 Median_Aquaslope <- read_in$Median_Aquaslope Mean_Aquaslope <- read_in$Mean_Aquaslope Median_Aquaslopval<-read_in$Median_Aquaslopval Mean_Aquaslopval<-read_in$Mean_Aquaslopval Median_Aqua_Intercept<-read_in$Median_Aqua_Intercept Mean_Aqua_Intercept<-read_in$Mean_Aqua_Intercept Median_Aquainpval<-read_in$Median_Aquainpval Mean_Aquainpval<-read_in$Mean_Aquainpval Median_Terrestrialslope<-Median_Terrestrialslope Mean_Terrestrialslope<-read_in$Mean_Terrestrialslope Median_Terrestrialslopval<-read_in$Median_Terrestrialslopval Mean_Terrestrialslopval<-read_in$Mean_Terrestrialslopval Median_Terrestrial_Intercept<-read_in$Median_Terrestrial_Intercept Mean_Terrestrial_Intercept<-read_in$Mean_Terrestrial_Intercept Median_Terrestrialintpval<-read_in$Median_Terrestrialintpval Mean_Terrestrialintpval<-read_in$Mean_Terrestrialintpval Volant_Median_slope<-read_in$Volant_Median_slope Volant_Mean_slope<-read_in$Volant_Mean_slope Volant_Median_slopval<-read_in$Volant_Median_slopval Volant_Mean_slopval<-read_in$Volant_Mean_slopval Volant_Median_Intercept<-read_in$Volant_Median_Intercept Volant_Mean_Intercept<-read_in$Volant_Mean_Intercept Volant_Int_MedPval<-read_in$Volant_Int_MedPval Volant_INT_MeaPval<-read_in$Volant_INT_MeaPval Median_Fossorialslope <- read_in$Median_Fossorialslope Mean_Fossorialslope <- read_in$Mean_Fossorialslope Median_Fossorialslopval<-read_in$Median_Fossorialslopval Mean_Fossoriallopval<-read_in$Mean_Fossorialslopval Median_Fossorial_Intercept<-read_in$Median_Fossorial_Intercept Mean_Fossorial_Intercept<-read_in$Mean_Fossorial_Intercept Median_Fossorialinpval<-read_in$Median_Fossorialinpval Mean_Fossorialinpval<-read_in$Mean_Fossorialinpval #plotting the final regression plot(PC1_RSI~logbodymass, data=data,main="Carnivore_Locomotion",xlab="lnmass",ylab="RSI") #now add the lines abline(carnivore.pgls, col="green") abline(a=Median_Terrestrial_Intercept,b=Median_Terrestrialslope,col="red") abline(a=Median_Aqua_Intercept,b=Median_Aquaslope, col="blue") abline(a=Volant_Median_Intercept,b=Volant_Median_slope,col="purple") abline(a=Median_Fossorial_Intercept,b=Median_Fossorialslope,col="grey") legend("topright", legend= sort(unique(specific_diet)), col=c("red","purple","blue"), pch=19, pt.cex = 3) #adding colors to the points qqplot1<-qplot(logbodymass,PC1_RSI,data=data,colour=Locomotion,xlab="lnmass",ylab="PC1_RSI",main="All Carnivore Locomotion",geom=c("point")) qqplot1+geom_abline(intercept=Median_Terrestrial_Intercept,slope=Median_Terrestrialslope,colour="red")+geom_abline(intercept=Median_Aqua_Intercept,slope=Median_Aquaslope,colour="blue")+geom_abline(intercept=Volant_Median_Intercept,slope=Volant_Median_slope,colour="purple")+geom_abline(intercept=Median_Fossorial_Intercept,slope=Median_Fossorialslope,colour="orange") # adding distinct colors to the lines in gg plot df<-as.data.frame(data) ggplot(data=df, aes (x=lnmass, y= PC1, group=Locomotion, color=Locomotion))+ geom_abline(intercept = Median_Terrestrial_Intercept, slope = Median_Terrestrialslope,color="blue")+geom_abline(intercept=Median_Aqua_Intercept,slope=Median_Aquaslope,colour="red")+geom_abline(intercept=Volant_Median_Intercept,slope=Volant_Median_slope,colour="purple")+ geom_abline(intercept = Median_Fossorial_Intercept, slope = Median_Fossorialslope,color="green")+ geom_point(aes(shape=Locomotion))+ theme_bw()+ theme_classic()+ theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ labs(y="Reproductive Strategy Index (PC1)", x="lnmass",title = " Carnivores by Mode of Locomotion") # Calculating standard error for each relevant stat #calculating standard error for each relevant statistic SE_Panco1<-sd(PanCov.res1)/sqrt(length(PanCov.res1)) SE_Panco2<-sd(PanCov.res2)/sqrt(length(PanCov.res2)) SE_AquaInpval<-sd(AquaticInt.pval_list)/sqrt(length(AquaticInt.pval_list)) SE_AquaSlopval<-sd(AquaticSlope.pval_list)/sqrt(length(AquaticSlope.pval_list)) SE_Aquaslope<-sd(Aquatic.slope_list)/sqrt(length(Aquatic.slope_list)) SE_Aquain<-sd(Aquatic.intercept_list)/sqrt(length(Aquatic.intercept_list)) SE_TerrestInpval<-sd(TerrestrialInt.pval_list)/sqrt(length(TerrestrialInt.pval_list)) SE_TerrestSlopval<-sd(TerrestrialSlope.pval_list)/sqrt(length(TerrestrialSlope.pval_list)) SE_Terrestslope<-sd(Terrestrial.slope_list)/sqrt(length(Terrestrial.slope_list)) SE_Terrestint<-sd(Terrestrial.intercept_list)/sqrt(length(Terrestrial.intercept_list)) SE_VolInpval<-sd(VolantInt.pval_list)/sqrt(length(VolantInt.pval_list)) SE_VolSlopval<-sd(VolantSlope.pval_list)/sqrt(length(VolantSlope.pval_list)) SE_Volslope<-sd(Volant.slope_list)/sqrt(length(Volant.slope_list)) SE_Volinter<-sd(Volant.intercept_list)/sqrt(length(Volant.intercept_list)) SE_FossInpval<-sd(FossorialInt.pval_list)/sqrt(length(FossorialInt.pval_list)) SE_FossSlopval<-sd(FossorialSlope.pval_list)/sqrt(length(FossorialSlope.pval_list)) SE_FossSlop<-sd(Fossorial.slope_list)/sqrt(length(Fossorial.slope_list)) SE_FossInt<-sd(Fossorial.intercept_list)/sqrt(length(Fossorial.intercept_list)) SE_Carnslopval<-sd(mass.pval_list)/sqrt(length(mass.pval_list)) SE_carninpval<-sd(int.pval_list)/sqrt(length(int.pval_list)) SE_carnslop<-sd(mass.slope_list)/sqrt(length(mass.slope_list)) SE_carnin<-sd(intercept_list)/sqrt(length(intercept_list)) # collecting t-stats for each group slope and intercept carnslopetstat<-median(mass.tstat_list) carninttstat<-median(int.tstat_list) Aquaslopetsta<-median(AquaticSlope.tstat_list) Aquaintstat<-median(AquaticInt.tstat_list) Fossorialslotstat<-median(FossorialSlope.tstat_list) Fossorialintstat<-median(FossorialInt.tstat_list) Terrestrialslotstat<-median(TerrestrialSlope.tstat_list) Terrestrialintstat<-median(TerrestrialInt.tstat_list) Volantslotstat<-median(VolantSlope.tstat_list) Volantintstat<-median(VolantInt.tstat_list) # plotting each individual iteration (input #'s 1-101 for slopes and intercepts) df<-as.data.frame(data) ggplot(data=df, aes (x=lnmass, y= PC1, group=Locomotion, color=Locomotion))+ geom_abline(intercept =Terrestrial.intercept_list[1], slope = Terrestrial.slope_list[1],color="blue")+geom_abline(intercept=Aquatic.intercept_list[1],slope=Aquatic.slope_list[1],colour="red")+geom_abline(intercept=Volant.intercept_list[1],slope=Volant.slope_list[1],colour="purple")+ geom_abline(intercept = Fossorial.intercept_list[1], slope = Fossorial.slope_list[1],color="green")+ geom_point(aes(shape=Locomotion))+ theme_bw()+ theme_classic()+ theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ labs(y="Reproductive Strategy Index (PC1)", x="lnmass",title = " Carnivores by Mode of Locomotion")