################################################################################################################## ####### BRT Zoonotic Results (Residuals + Full) ####### Last Edit: Aug 17th, 2019 ################################################################################################################## ############################################################### #Loading Data ############################################################### #Loading in measures of variable importance vimp.full<-read.csv("Relative_Importance_Matrix_5000.csv", header=T) vimp.res<-read.csv("Relative_Importance_Matrix_5001.csv", header =T) #Loading in measures of accuarcy acc.full<-read.csv("Accuracy_Matrix_5000.csv", header=T) acc.res<-read.csv("Accuracy_Matrix_5001.csv", header=T) ############################################################### #Realtive Importance of Variables - VALUES ############################################################### #Need to transpose measures of varaible importance to make further manipulation easier #Reording so matrix is vertical rather than horzontal and column names are the variables t.vimpfull<-t(vimp.full[,c(2:200)]) colnames(t.vimpfull)<-vimp.full$X t.vimpres<-t(vimp.res[,c(2:200)]) colnames(t.vimpres)<-vimp.res$X #**********SCIENCE NOTES: I need to create a shaded polygon in the back of my plot that reflects the 95% range of realtive importance measures from the different analyses. #This means that I need to order my relative importance measures and find the lower and upper bounds. The lower bound will be 2.5% (0.025*200=5th measure) and the upper bound will be #97.5% (0.975*200=195th measure)). For the graph I need to pull out the average value, and the upper and lower bounds - for each variable. #Following function is meant to get values need to create the graph and the return them as a matrix ordering.vals<-function(orig.matrix, t.matrix,x){ #Creating an empty matrix to store the reuslts results<-matrix(nrow = length(orig.matrix$X), ncol = 3) colnames(results)<-c("Mean", "Lower", "Upper") rownames(results)<-orig.matrix$X #Ordering the measures of realtive importance from smallest to largest. #Doing this by using the apply function to loop over all columns (2) in the transposed matrix sorted<-apply(t.matrix, 2, sort) #Assigning the upper and lower measures to the empty matrix using the sorted matrix above results[,2]<-sorted[0.025*200,] results[,3]<-sorted[0.975*200,] #Calcualting means of each column and assigning to empty matrix results[,1]<-apply(t.matrix,2,mean) #Exporting matrix as a csv file write.csv(results, paste("RI_Results",x,".csv", sep=""), row.names = T) #Function returns the filled resulsts matric return(results) } #Using above function to get the values that we need plotted full<-ordering.vals(orig.matrix = vimp.full, t.matrix = t.vimpfull, x=5000) dna<-ordering.vals(orig.matrix = vimp.res, t.matrix = t.vimpres, x=5001) ############################################################### #Training and test accuracy - VALUES ############################################################### #Following function is meant to order values of accuracy to get the upper and lower bounds as well as the average ordering.acc<-function(acc.matrix,x){ #Creating an empty matrix to store the reuslts results<-matrix(nrow = 3, ncol = 3) colnames(results)<-c("Mean", "Lower", "Upper") rownames(results)<-c("x", "Training", "Test") #Ordering the accuracy measures from smallest to largest #Doing this by using the apply function to loop over all columns (2) in the transposed matrix sorted<-apply(acc.matrix, 2, sort) #Assigning the upper and lower measures to the empty matrix using the sorted matrix above results[,2]<-sorted[0.025*200,] results[,3]<-sorted[0.975*200,] #Calcualting means of each column and assigning to empty matrix results[,1]<-apply(acc.matrix,2,mean) #Exporting matrix as a csv file write.csv(results, paste("Acc_Results",x,".csv", sep=""), row.names = T) #Function returns the filled resulsts matric return(results) } #using the above function to get the values that need to be plotted for the accuracy graphy acc.full<-ordering.acc(acc.matrix = acc.full, x=5000) acc.res<-ordering.acc(acc.matrix = acc.res, x=5001) ############################################################### #Plotting Graphs - Relative Importance ############################################################### #Loading back in csv files RI.full<-read.csv("RI_Results5000.csv", header=T) RI.res<-read.csv("RI_Results5001.csv", header=T) #Plotting the measures of realtive importance. Before I plot though, I am ordering the dataframe in decreasing values #Ordering RI for plotting sRI.full<-RI.full[order(-RI.full$Mean),] sRI.res<-RI.res[order(-RI.res$Mean),] #The variables for family and full are not ordered the same. To fix this I'm going to export as a csv and order the full by hand to #match the family model write.csv(sRI.full, "RI_full_zoo.csv", row.names = T) write.csv(sRI.res, "RI_residuals.csv", row.names = T) #Loading in the dataset that I reordered by hand. There is likely a much smarter way to do this...probably by matching the two lists #and ordering that way. sRI.residuals.plot<-read.csv("RI_residuals_T.csv", header=T) #Loading in accuracy measures plotacc.full<-read.csv("Acc_Results5000.csv", header = T) plotacc.residuals<-read.csv("Acc_Results5001.csv", header=T) ############################################################### ################################################# #######RI With 2 pannels no cit residuals ################################################# #Using dev.new to open a new graphic window dev.new() #Specifying that plot will be saved as a pdf with the following name #using pdf because it can be resized pdf(file = "Full_&_NoCit_Models_Residuals_Zoonotic.pdf", 15,9) par(mfrow=c(2,1)) #Seting up plot with points for display par(mai=c(1.2,1.0,.05,0.05)) plot(sRI.full$Mean, pch=19, axes=F, ylim = c(0,70), xlab = "", ylab= "") #Variable names for the x axsis ecovars.fam<-c("No. Citations", "Family", "Group Size", "Lifespan", "Range Area", "Longitude", "Relative \nWing Loading", "Torpor", "Aspect Ratio", "Sympatry \nResiduals", "Forearm", "Latitude", "No. Birthing \nPulses", "No. Mixed \nAssociations", "Migration", "Roost Cat.", "IUCN Trend", "IUCN Status", "Diet") #Setting up the x and y axes axis(1, at=1:19, labels=F, lwd = 1, lwd.ticks = 1) text(seq(0.5,18.5, by=1), -12.0, labels = ecovars.fam, srt = 45, pos =1.0, xpd = TRUE, cex = 1.1) axis(2, cex.axis=1.1, at=seq(0,70, by=5)) #Adding in rectangles of measures of realtive importance rect(xleft = c(1:19)+0.2, xright = c(1:19)-0.2, ybottom = sRI.full$Lower, ytop = sRI.full$Upper, col=adjustcolor("grey10", alpha.f = 0.3), border=NA) #Labelling Axes box(bty="l") #Indicating which cluster solution this is text(16.5,55 , "Zoonotic Viral Families \nFamily Clusters \nwith Citations", cex=1.3) #Adding in Y axis labe mtext("Relative Importance (%)", side=2, cex=1.4, line = 3.5) #Setting up the Residuals par(mai=c(1.5,1.0,.05,0.05)) plot(sRI.residuals.plot$Mean, pch=19, axes=F, ylim = c(0,70), xlab = "", ylab= "") ecovars.fam.resid<-c("No. Citations", "Family", "Group Size \nResiduals*", "Lifespan \nResiduals*", "Range Area \nResiduals*", "Longitude \nResiduals*", "Relative \nWing Loading", "Torpor", "Aspect Ratio", "Sympatry \nResiduals", "Forearm", "Latitude \nResiduals*", "No. Birthing \nPulses", "No. Mixed \nAssociations", "Migration", "Roost Cat.", "IUCN Trend", "IUCN Status", "Diet") #Setting up the x and y axes axis(1, at=1:19, labels=F, lwd = 1, lwd.ticks = 1) text(seq(0.5,18.5, by=1), -12.0, labels = ecovars.fam.resid, srt = 45, pos =1.0, xpd = TRUE, cex = 1.1) axis(2, cex.axis=1.1, at=seq(0,70, by=5)) #Adding in rectangles of measures of realtive importance rect(xleft = c(1:19)+0.2, xright = c(1:19)-0.2, ybottom = sRI.residuals.plot$Lower, ytop = sRI.residuals.plot$Upper, col=adjustcolor("darkgreen", alpha.f = 0.3), border=NA) #Labelling Axes box(bty="l") #Indicating which cluster solution this is text(16.5,55 , "Zoonotic Viral Familes \nFamily Clusters \nCitations Residuals", cex = 1.3) #Adding in Y axis labe mtext("Relative Importance (%)", side=2, cex=1.4, line = 3.5) #Adding text for X and Y axes mtext("Bat Traits", side = 1, cex = 1.4, line = 5.7) #Turning off graphics output dev.off() ############################################################### #Citations/No Citations Accuracy ############################################################### #Using dev.new to open a new graphic window dev.new() #Specifying that plot will be saved as a pdf with the following name #using pdf because it can be resized pdf(file= "ZooCitations_NoCitations_Acc.pdf", 13, 9) #Adjusting plotting window margins par(mfrow=c(1,2)) #Variable names for the x axsis acc<-c("Training Data", "Test Data") #Seting up plot with points for display par(mai=c(0.85,1.0,.05,0.1)) plot(c(0.2,0.6),c(plotacc.full$Mean),ylim = c(0,1),xlim=c(0,1), axes=F, xlab="", ylab = "", col=NA, xaxs="i") #Adding in rectangles of measures of realtive importance rect(xleft = c(0.2,0.6)+0.02, xright = c(0.2,0.6)-0.02, ybottom = plotacc.full$Lower, ytop = plotacc.full$Upper, col=adjustcolor("grey10", alpha.f = 0.3), border=NA) points(c(0.2,0.6),c(plotacc.full$Mean),pch=19, cex=1, col="black") axis(1, at=c(0.2,0.6),labels = acc, cex.axis=1.0, xpd=T) axis(2, at=seq(0,1.0, by=0.1), cex.axis=1.0) title(ylab = "Pseudo R2 Value", line =3.0, cex.lab=1.1) #Indicating which cluster solution this is text(0.7,0.95 , "Zoonotic Viral Families \nwith Publications", cex = 1.1) #Plotting 41MYA Clusters par(mai=c(0.85,0.4,.05,0.1)) plot(c(0.2,0.6),c(plotacc.residuals$Mean),ylim = c(0,1),xlim=c(0,1), axes=F, xlab="", ylab = "", col=NA, xaxs="i") #Adding in rectangles of measures of realtive importance rect(xleft = c(0.2,0.6)+0.02, xright = c(0.2,0.6)-0.02, ybottom = plotacc.residuals$Lower, ytop = plotacc.residuals$Upper, col=adjustcolor("darkgreen", alpha.f = 0.3), border=NA) points(c(0.2,0.6),c(plotacc.residuals$Mean),pch=19, cex=1, col="black") axis(1, at=c(0.2,0.6),labels = acc, cex.axis=1.0, xpd=T) axis(2, at=seq(0,1.0, by=0.1), cex.axis=1.0) #Indicating which cluster solution this is text(0.7,0.95 , "Zoonotic Viral Familes \n Publications Residuals", cex = 1.1) #Turning off graphics output dev.off() ############################################################################################