diff scripts_R/PlotOutliersCountry.R @ 77:bde45ce0eeab branch-tests

plots and figures for results

--- a/scripts_R/PlotOutliersCountry.R	Fri Sep 22 16:30:36 2017 +0100
+++ b/scripts_R/PlotOutliersCountry.R	Fri Sep 22 18:02:59 2017 +0100
@@ -1,19 +1,23 @@
 source("MetadataPlots.R")
 
-PlotCountryOutliers(df=read.csv("data/global_outliers.csv",header=TRUE), output="data/global_outliers.pdf")
-PlotCountryOutliers(df=read.csv("data/global_outliers_rhy.csv",header=TRUE), output="data/global_outliers_rhy.pdf")
-PlotCountryOutliers(df=read.csv("data/global_outliers_mel.csv",header=TRUE), output="data/global_outliers_mel.pdf")
-PlotCountryOutliers(df=read.csv("data/global_outliers_mfc.csv",header=TRUE), output="data/global_outliers_mfc.pdf")
-PlotCountryOutliers(df=read.csv("data/global_outliers_chr.csv",header=TRUE), output="data/global_outliers_chr.pdf")
-PlotCountryOutliers(df=read.csv("data/spatial_outliers.csv",header=TRUE), output="data/spatial_outliers.pdf")
-#PlotCountryOutliers(df=read.csv("data/global_outliers_rhy_1band.csv",header=TRUE))
+PlotCountryOutliers(df=read.csv("../data/results/global_outliers.csv",header=TRUE), output="../data/results/global_outliers.pdf")
+PlotCountryOutliers(df=read.csv("../data/results/global_outliers_rhy.csv",header=TRUE), output="../data/results/global_outliers_rhy.pdf")
+PlotCountryOutliers(df=read.csv("../data/results/global_outliers_mel.csv",header=TRUE), output="../data/results/global_outliers_mel.pdf")
+PlotCountryOutliers(df=read.csv("../data/results/global_outliers_mfc.csv",header=TRUE), output="../data/results/global_outliers_mfc.pdf")
+PlotCountryOutliers(df=read.csv("../data/results/global_outliers_chr.csv",header=TRUE), output="../data/results/global_outliers_chr.pdf")
+PlotCountryOutliers(df=read.csv("../data/results/spatial_outliers.csv",header=TRUE), output="../data/results/spatial_outliers.pdf")
 
-require(graphics)
-par(mfrow=c(2,2))
-g1<-PlotCountryOutliers(df=read.csv("data/global_outliers_rhy.csv",header=TRUE))
-g2<-PlotCountryOutliers(df=read.csv("data/global_outliers_mel.csv",header=TRUE))
-g3<-PlotCountryOutliers(df=read.csv("data/global_outliers_mfc.csv",header=TRUE))
-g4<-PlotCountryOutliers(df=read.csv("data/global_outliers_chr.csv",header=TRUE))
-#do.call(addMapLegend, c(g3,labelFontSize=0.7, legendWidth=0.5, tcl=0.3, legendMar = 7, legendLabels="all",horizontal=T, legendIntervals="page"))
-#legend("bottomleft", legend = c(paste(seq(100,1,-10),'%'), 'missing countries'), fill = c(heat.colors(10, alpha = 1), 'grey'), cex = 0.56, bty = "n")
-legend("right", legend = c(paste(seq(90,0,-10),'-',seq(100,10,-10),'%'), 'NA'), fill = c(heat.colors(10, alpha = 1), 'grey'), cex = 0.56, bty = "o",bg="white",box.lwd=0,box.col="white")
+library(ape)
+library(cluster) 
+
+df = read.csv("../data/results/cluster_freq.csv")
+data = df[,2:dim(df)[2]]
+rownames(data) <- df$labels
+distMahal = as.dist(apply(data, 1, function(i) mahalanobis(data, i, cov = cov(data),tol=1e-18)))
+hc=hclust(distMahal, method="average")
+mypal = c("#000000", "#9B0000", "#9B0000", "#9B0000", "#9B0000")
+clus5 = cutree(hc, 5)
+pdf('../data/results/hierarchical_cluster.pdf')
+par(mar=c(1,1,1,1))
+plot(as.phylo(hc),type="fan",tip.color=mypal[clus5], cex=.5, label.offset=.5)
+dev.off()