Mercurial > hg > plosone_underreview

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
plots and figures for results
author Maria Panteli <m.x.panteli@gmail.com>
date Fri, 22 Sep 2017 18:02:59 +0100
parents cc028157502a
children 103f7411c3ad
comparison
equal deleted inserted replaced
76:d17833be50ca 77:bde45ce0eeab
1 source("MetadataPlots.R") 1 source("MetadataPlots.R")
2 2
3 PlotCountryOutliers(df=read.csv("data/global_outliers.csv",header=TRUE), output="data/global_outliers.pdf") 3 PlotCountryOutliers(df=read.csv("../data/results/global_outliers.csv",header=TRUE), output="../data/results/global_outliers.pdf")
4 PlotCountryOutliers(df=read.csv("data/global_outliers_rhy.csv",header=TRUE), output="data/global_outliers_rhy.pdf") 4 PlotCountryOutliers(df=read.csv("../data/results/global_outliers_rhy.csv",header=TRUE), output="../data/results/global_outliers_rhy.pdf")
5 PlotCountryOutliers(df=read.csv("data/global_outliers_mel.csv",header=TRUE), output="data/global_outliers_mel.pdf") 5 PlotCountryOutliers(df=read.csv("../data/results/global_outliers_mel.csv",header=TRUE), output="../data/results/global_outliers_mel.pdf")
6 PlotCountryOutliers(df=read.csv("data/global_outliers_mfc.csv",header=TRUE), output="data/global_outliers_mfc.pdf") 6 PlotCountryOutliers(df=read.csv("../data/results/global_outliers_mfc.csv",header=TRUE), output="../data/results/global_outliers_mfc.pdf")
7 PlotCountryOutliers(df=read.csv("data/global_outliers_chr.csv",header=TRUE), output="data/global_outliers_chr.pdf") 7 PlotCountryOutliers(df=read.csv("../data/results/global_outliers_chr.csv",header=TRUE), output="../data/results/global_outliers_chr.pdf")
8 PlotCountryOutliers(df=read.csv("data/spatial_outliers.csv",header=TRUE), output="data/spatial_outliers.pdf") 8 PlotCountryOutliers(df=read.csv("../data/results/spatial_outliers.csv",header=TRUE), output="../data/results/spatial_outliers.pdf")
9 #PlotCountryOutliers(df=read.csv("data/global_outliers_rhy_1band.csv",header=TRUE))
10 9
11 require(graphics) 10 library(ape)
12 par(mfrow=c(2,2)) 11 library(cluster)
13 g1<-PlotCountryOutliers(df=read.csv("data/global_outliers_rhy.csv",header=TRUE)) 12
14 g2<-PlotCountryOutliers(df=read.csv("data/global_outliers_mel.csv",header=TRUE)) 13 df = read.csv("../data/results/cluster_freq.csv")
15 g3<-PlotCountryOutliers(df=read.csv("data/global_outliers_mfc.csv",header=TRUE)) 14 data = df[,2:dim(df)[2]]
16 g4<-PlotCountryOutliers(df=read.csv("data/global_outliers_chr.csv",header=TRUE)) 15 rownames(data) <- df$labels
17 #do.call(addMapLegend, c(g3,labelFontSize=0.7, legendWidth=0.5, tcl=0.3, legendMar = 7, legendLabels="all",horizontal=T, legendIntervals="page")) 16 distMahal = as.dist(apply(data, 1, function(i) mahalanobis(data, i, cov = cov(data),tol=1e-18)))
18 #legend("bottomleft", legend = c(paste(seq(100,1,-10),'%'), 'missing countries'), fill = c(heat.colors(10, alpha = 1), 'grey'), cex = 0.56, bty = "n") 17 hc=hclust(distMahal, method="average")
19 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") 18 mypal = c("#000000", "#9B0000", "#9B0000", "#9B0000", "#9B0000")
19 clus5 = cutree(hc, 5)
20 pdf('../data/results/hierarchical_cluster.pdf')
21 par(mar=c(1,1,1,1))
22 plot(as.phylo(hc),type="fan",tip.color=mypal[clus5], cex=.5, label.offset=.5)
23 dev.off()