# Carnivore designations: INV=Hypocarnivore, Mixed=mesocarnivore, Vert=Hypercarnivore #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 marine mammals #dietnowhale<-data[which(data$group!="whale"),] #take out the flying mammals #terrestrial_carnivores<-data[which(data$group!="whale" & data$group!="bat" & data$group!="otter"),] #data<-terrestrial_carnivores 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<- Hyper.intercept_list<- HyperInt.tstat_list<- HyperInt.pval_list<- Hypo.intercept_list<- HypoInt.tstat_list<- HypoInt.pval_list<- Meso.intercept_list<- MesoInt.tstat_list<- MesoInt.pval_list<- Hyper.slope_list<- HyperSlope.tstat_list<- HyperSlope.pval_list<- Hypo.slope_list<- HypoSlope.tstat_list<- HypoSlope.pval_list<- Meso.slope_list<- MesoSlope.tstat_list<- MesoSlope.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)) { print(i) 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. specific_diet<-data$specific_diet names(specific_diet) <- rownames(data) Diet<-data$Diet #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 diet") #legend("topright", legend= sort(unique(specific_diet)), col=c("green","blue","purple", "yellow"), pch=19, pt.cex = 3) #Now we can plot the character traits as colored 'pies' at the tips: #cols<-setNames(c("green","blue","purple", "yellow"), sort(unique(specific_diet))) #tiplabels(pie=to.matrix(specific_diet,sort(unique(specific_diet))),piecol=cols,cex=0.5) #Continuous character maps This maps continuous traits on the phylogenetic tree. With these graphics we are able to visualize #how the dimensions of body size and or reproductive index are distributed among species. #making named vector logbodymass<-data$lnmass names(logbodymass) <- rownames(data) #making a contour map obj<-contMap(pruned.trees, logbodymass, plot=FALSE) #plot(obj,type="phylogram",leg.txt="lnmass",lwd=6, mar=c(4,2,4,2), outline=FALSE) #title("carnivore lnmass") #same code but for PC1..RSI PC1_RSI <-data$PC1 names(PC1_RSI) <- rownames(data) obj2<-contMap(pruned.trees, PC1_RSI, plot=FALSE) #plot(obj2,type="phylogram",leg.txt="PC1",lwd=6, mar=c(4,2,4,2), outline=FALSE) #title("Carnivore PC1") #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") #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 * specific_diet, data=carnivore) summary(categories.pgls) categories1.pgls <- pgls(PC1~lnmass + specific_diet, 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? Hypo<-subset(data,specific_diet=="INV",select = c("Taxon","lnmass","PC1")) Meso<-subset(data,specific_diet=="Mixed",select = c("Taxon","lnmass","PC1")) Hyper<-subset(data,specific_diet=="Vert",select = c("Taxon","lnmass","PC1")) # generate pgls for 3 diet subcategories Hypo_Carnivores<- comparative.data(pruned.trees, Hypo, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) Meso_Carnivores<- comparative.data(pruned.trees, Meso, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) Hyper_Carnivores<- comparative.data(pruned.trees, Hyper, names.col=Taxon, vcv=TRUE, na.omit=FALSE, warn.dropped=TRUE) #-----?? data dropped? #Invert pgls Hypo.pgls<-pgls(PC1~lnmass, data=Hypo_Carnivores, lambda="ML") summary(Hypo.pgls) ipgls<-summary(Hypo.pgls) #Mixed pgls Meso.pgls<-pgls(PC1~lnmass, data=Meso_Carnivores, lambda="ML") summary(Meso.pgls) mpgls<-summary(Meso.pgls) #Vert pgls Hyper.pgls<-pgls(PC1~lnmass, data=Hyper_Carnivores, lambda="ML") summary(Hyper.pgls) vpgls<-summary(Hyper.pgls) #create models for the 3 diet categories model2<-lm(PC1~lnmass, data=Hypo) summary(model2) model3<-lm(PC1~lnmass, data=Meso) summary(model3) model4<-lm(PC1~lnmass, data=Hyper) summary(model4) #Plot original uninformed/pgls regression comparissons for 3 diet subcategories #Invertebrate feeders #plot(PC1~lnmass, data=Hypo) #abline(model2, col="green") #title("Original regression Hypo_Carnivores") #plot PGLS model(InV_feeders) #plot(PC1~lnmass, data=Hypo) #abline(Hypo.pgls, col="green") #title("PGLS_Hypo_carnivores") #Mixed Feeders #plot(PC1~lnmass, data=Meso) # abline(model3, col="blue") # title("Original regression Meso_Carnivores") #plot PGLS model(Mixed_feeders) #plot(PC1~lnmass, data=Meso) #abline(Meso.pgls, col="blue") #title("PGLS_Meso_Carnivores") #Vert Feeders #plot(PC1~lnmass, data=Hyper) #abline(model4, col="red") #title("Original regression Hyper_Carnivores") #plot PGLS model(Mixed_feeders) #plot(PC1~lnmass, data=Hyper) #abline(Hyper.pgls, col="red") # title("PGLS_Hyper_Carnivores") #plot all PGLS on same axis # plot(PC1~lnmass,data=data) # abline(carnivore.pgls,col="white")+ # abline(Hypo.pgls,col="green")+ # abline(Meso.pgls,col="blue")+ # abline(Hyper.pgls,col="red")+ # title("pgls_carnivore_regression") # legend("topright", legend= sort(unique(specific_diet)), col=c("green","blue","red"), 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(model3,col="purple")+ # abline(model4,col="blue")+ # title("uninformed_carnivore_regression") #legend("topright", legend= sort(unique(specific_diet)), col=c("green","blue","red"), 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] Hyper.intercept_list[i]<-vpgls$coefficients[1,1] HyperInt.tstat_list[i]<-vpgls$coefficients[1,3] HyperInt.pval_list[i]<-vpgls$coefficients[1, 4] Hypo.intercept_list[i]<- ipgls$coefficients[1,1] HypoInt.tstat_list[i]<-ipgls$coefficients[1,3] HypoInt.pval_list[i]<-ipgls$coefficients[1, 4] Meso.intercept_list[i]<- mpgls$coefficients[1,1] MesoInt.tstat_list[i]<-mpgls$coefficients[1,3] MesoInt.pval_list[i]<-mpgls$coefficients[1, 4] #update the other slopes like this vert slope Hyper.slope_list[i]<-vpgls$coefficients[2,1] HyperSlope.tstat_list[i]<-vpgls$coefficients[2,3] HyperSlope.pval_list[i]<-vpgls$coefficients[2,4] Hypo.slope_list[i]<- ipgls$coefficients[2,1] HypoSlope.tstat_list[i]<-ipgls$coefficients[2,3] HypoSlope.pval_list[i]<-ipgls$coefficients[2, 4] Meso.slope_list[i]<- mpgls$coefficients[2,1] MesoSlope.tstat_list[i]<-mpgls$coefficients[2,3] MesoSlope.pval_list[i]<-mpgls$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,1),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 Vert carnivore slope, intercept, and P-values for each hist(Hyper.slope_list,main="Hypercarnivore slopes",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(-0.23,-0.22),ylim=c(0,6),col="red",breaks=100) Median_Hyperslope<-median(Hyper.slope_list) Mean_Hyperslope<-mean(Hyper.slope_list) #Vert_slope P- values hist(HyperSlope.pval_list,main="Hyper_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.0011,0.0017),ylim=c(0,6),col="red",breaks=100) Median_Hyperslopval<-median(HyperSlope.pval_list) Mean_Hyperslopval<-mean(HyperSlope.pval_list) # Now Vert carnivore intercepts hist(Hyper.intercept_list,main="Hypercarnivore Intercepts",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(1.06,1.2),ylim=c(0,6),col="red",breaks=100) Median_Hyper_Intercept<-median(Hyper.intercept_list) Mean_Hyper_Intercept<-mean(Hyper.intercept_list) #Vert carnivore Intercept P-values hist(HyperInt.pval_list,main="Hyper_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.13,0.17),ylim=c(0,6),col="red",breaks=100) Median_Hyperinpval<-Median_Vertinpval<-median(HyperInt.pval_list) Mean_Hyperinpval<-mean(HyperInt.pval_list) # do the same for Invertebrate carnivores hist(Hypo.slope_list,main="Hypocarnivore slopes",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(-0.23,-0.24),ylim=c(0,6),col="blue",breaks=100) Median_Hyposlope<-median(Hypo.slope_list) Mean_Hyposlope<-mean(Hypo.slope_list) #InVert_slope P- values hist(HypoSlope.pval_list,main="Hypo_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(8.75E-12,1.75E-9),ylim=c(0,15),col="blue",breaks=100) Median_Hyposlopval<-median(HypoSlope.pval_list) Mean_Hyposlopval<-mean(HypoSlope.pval_list) # Now InVert carnivore intercepts hist(Hypo.intercept_list,main="Hypocarnivore Intercepts",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(1.4,1.48),ylim=c(0,10),col="blue",breaks=100) Median_Hypo_Intercept<-median(Hypo.intercept_list) Mean_Hypo_Intercept<-mean(Hypo.intercept_list) #InVert carnivore Intercept P-values hist(HypoInt.pval_list,main="Hypo_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(0.004,0.01),ylim=c(0,8),col="blue",breaks=100) Median_Hypointpval<-median(HypoInt.pval_list) Mean_Hypointpval<-mean(HypoInt.pval_list) #same thing for mixed carnivores hist(Meso.slope_list,main="Mesocarnivore slopes",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(-0.4858943,-0.4858943),ylim=c(0,8),col="purple",breaks=100) Meso_Median_slope<-median(Meso.slope_list) Meso_Mean_slope<-mean(Meso.slope_list) #Mixed_slope P- values hist(MesoSlope.pval_list,main="Meso_slope p-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(1.761702e-12,1.761702e-12),ylim=c(0,101),col="purple",breaks=100) Meso_Median_slopval<-median(MesoSlope.pval_list) Meso_Mean_slopval<-mean(MesoSlope.pval_list) # Now Mixed carnivore intercepts hist(Meso.intercept_list,main="Mesocarnivore Intercepts",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(2.904124,2.904124),ylim=c(0,6),col="purple",breaks=100) Meso_Median_Intercept<-median(Meso.intercept_list) Meso_Mean_Intercept<-mean(Meso.intercept_list) #Mixed carnivore Intercept P-values hist(MesoInt.pval_list,main="Meso_intercept P-values",border="black",xlab="P_values ",ylab="number_of_occurences",xlim=c(1.356481e-08,1.356481e-08),ylim=c(0,10),col="purple",breaks=100) Meso_Int_MedPval<-median(MesoInt.pval_list) Meso_INT_MeaPval<-mean(MesoInt.pval_list) #making a df of the summary results SummaryResults <- data.frame(Median_Panco1, Mean_Panco1, Median_Panco2, Mean_Panco2, Median_Hyperslope, Mean_Hyperslope, Median_Hyperslopval, Mean_Hyperslopval, Median_Hyper_Intercept, Mean_Hyper_Intercept, Median_Hyperinpval, Mean_Hyperinpval, Median_Hyposlope, Mean_Hyposlope, Median_Hyposlopval, Mean_Hyposlopval, Median_Hypo_Intercept, Mean_Hypo_Intercept, Median_Hypointpval, Mean_Hypointpval, Meso_Median_slope, Meso_Mean_slope, Meso_Median_slopval, Meso_Mean_slopval, Meso_Median_Intercept, Meso_Mean_Intercept, Meso_Int_MedPval, Meso_INT_MeaPval) #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_Hyperslope <- read_in$Median_Hyperslope Mean_Hyperslope <- read_in$Mean_Hyperslope Median_Hyperslopval<-read_in$Median_Hyperslopval Mean_Hyperslopval<-read_in$Mean_Hyperslopval Median_Hyper_Intercept<-read_in$Median_Hyper_Intercept Mean_Hyper_Intercept<-read_in$Mean_Hyper_Intercept Median_Hyperinpval<-read_in$Median_Hyperinpval Mean_Hyperinpval<-read_in$Mean_Hyperinpval Median_Hyposlope<-read_in$Median_Hyposlope Mean_Hyposlope<-read_in$Mean_Hyposlope Median_Hyposlopval<-read_in$Median_Hyposlopval Mean_Hyposlopval<-read_in$Mean_Hyposlopval Median_Hypo_Intercept<-read_in$Median_Hypo_Intercept Mean_Hypo_Intercept<-read_in$Mean_Hypo_Intercept Median_Hypointpval<-read_in$Median_Hypointpval Mean_Hypointpval<-read_in$Mean_Hypointpval Meso_Median_slope<-read_in$Meso_Median_slope Meso_Mean_slope<-read_in$Meso_Mean_slope Meso_Median_slopval<-read_in$Meso_Median_slopval Meso_Mean_slopval<-read_in$Meso_Mean_slopval Meso_Median_Intercept<-read_in$Meso_Median_Intercept Meso_Mean_Intercept<-read_in$Meso_Mean_Intercept Meso_Int_MedPval<-read_in$Meso_Int_MedPval Meso_INT_MeaPval<-read_in$Meso_INT_MeaPval #plotting the final regression plot(PC1_RSI~logbodymass, data=data,main="Carnivore_Diets",xlab="lnmass",ylab="RSI") #now add the lines abline(carnivore.pgls, col="grey") abline(a=Median_Hypo_Intercept,b=Median_Hyposlope,col="green") abline(a=Median_Hyper_Intercept,b=Median_Hyperslope, col="red") abline(a=Meso_Median_Intercept,b=Meso_Median_slope,col="blue") legend("topright", legend= sort(unique(specific_diet)), col=c("green","blue","red"), pch=19, pt.cex = 3) #adding colors to the points qqplot1<-qplot(logbodymass,PC1_RSI,data=data,colour=Diet,xlab="lnmass",ylab="PC1_RSI",main="All Carnivore Diets",geom=c("point")) qqplot1+geom_abline(intercept=Median_Hypo_Intercept,slope=Median_Hyposlope,colour="green")+geom_abline(intercept=Median_Hyper_Intercept,slope=Median_Hyperslope,colour="red")+geom_abline(intercept=Meso_Median_Intercept,slope=Meso_Median_slope,colour="blue") # adding distinct colors to the lines in gg plot df<-as.data.frame(data) ggplot(data=df, aes (x=lnmass, y= PC1, group=Diet, color=Diet))+ geom_abline(intercept = Median_Hypo_Intercept, slope = Median_Hyposlope,color="green")+geom_abline(intercept=Median_Hyper_Intercept,slope=Median_Hyperslope,colour="red")+geom_abline(intercept=Meso_Median_Intercept,slope=Meso_Median_slope,colour="blue")+ 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="ln(mass)",title = " Carnivore Diet By Preferred Prey Organisms") # grabbing some data for the 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) #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_HyperInpval<-sd(HyperInt.pval_list)/sqrt(length(HyperInt.pval_list)) SE_HyperSlopval<-sd(HyperSlope.pval_list)/sqrt(length(HyperSlope.pval_list)) SE_Hyperslope<-sd(Hyper.slope_list)/sqrt(length(Hyper.slope_list)) SE_Hyperint<-sd(Hyper.intercept_list)/sqrt(length(Hyper.intercept_list)) SE_HypoInpval<-sd(HypoInt.pval_list)/sqrt(length(HypoInt.pval_list)) SE_HypoSlopval<-sd(HypoSlope.pval_list)/sqrt(length(HypoSlope.pval_list)) SE_Hyposlope<-sd(Hypo.slope_list)/sqrt(length(Hypo.slope_list)) SE_hypoint<-sd(Hypo.intercept_list)/sqrt(length(Hypo.intercept_list)) SE_MesoInpval<-sd(MesoInt.pval_list)/sqrt(length(MesoInt.pval_list)) SE_MesoSlopval<-sd(MesoSlope.pval_list)/sqrt(length(MesoSlope.pval_list)) SE_Mesoslope<-sd(Meso.slope_list)/sqrt(length(Meso.slope_list)) SE_Mesoint<-sd(Meso.intercept_list)/sqrt(length(Meso.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) Hyposlopetstat<-median(HypoSlope.tstat_list) Hypointstat<-median(HypoInt.tstat_list) Mesosloptstat<-median(MesoSlope.tstat_list) mesointstat<-median(MesoInt.tstat_list) Hyperslopetstat<-median(HyperSlope.tstat_list) Hyperintstat<-median(HyperInt.tstat_list) #Plotting lines for each iteration (input line # between 1 and 101 for slopes and intercepts) df<-as.data.frame(data) ggplot(data=df, aes (x=lnmass, y= PC1, group=Diet, color=Diet))+ geom_abline(intercept = Hypo.intercept_list[1], slope = Hypo.slope_list[1],color="green")+geom_abline(intercept=Hyper.intercept_list[1],slope=Hyper.slope_list[1],colour="red")+geom_abline(intercept=Meso.intercept_list[1],slope=Meso.slope_list[1],colour="blue")+ 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="ln(mass)",title = " Carnivore Diet By Preferred Prey Organisms")