Mercurial > hg > plosone_underreview

comparison scripts_R/PlotOutliersCountry.R @ 91:d3e05cd49feb branch-tests

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
trying to plot wrt PLOSONE guidelines
author Maria Panteli <m.x.panteli@gmail.com>
date Mon, 02 Oct 2017 15:32:51 +0100
parents 103f7411c3ad
children
comparison
equal deleted inserted replaced
89:8a2d56880050 91:d3e05cd49feb
10 library(ape) 10 library(ape)
11 library(cluster) 11 library(cluster)
12 12
13 df = read.csv("../data/results/cluster_freq.csv") 13 df = read.csv("../data/results/cluster_freq.csv")
14 data = df[,2:dim(df)[2]] 14 data = df[,2:dim(df)[2]]
15 levels(df$labels)[which(levels(df$labels)=="Democratic Republic of the Congo")]="DR Congo"
16 df$labels[which(df$labels=="Democratic Republic of the Congo")] = "DR Congo"
15 rownames(data) <- df$labels 17 rownames(data) <- df$labels
16 distMahal = as.dist(apply(data, 1, function(i) mahalanobis(data, i, cov = cov(data),tol=1e-18))) 18 distMahal = as.dist(apply(data, 1, function(i) mahalanobis(data, i, cov = cov(data),tol=1e-18)))
17 hc=hclust(distMahal, method="average") 19 hc=hclust(distMahal, method="average")
18 mypal = c("#000000", "#9B0000", "#9B0000", "#9B0000", "#9B0000") 20 mypal = c("#000000", "#9B0000", "#9B0000", "#9B0000", "#9B0000")
19 clus5 = cutree(hc, 4) 21 clus5 = cutree(hc, 4)
20 pdf('../data/results/hierarchical_cluster.pdf') 22 pdf('../data/results/hierarchical_cluster.pdf', pointsize=12)
21 par(mar=c(1,1,1,1)) 23 par(mar=c(1,1,1,1))
22 plot(as.phylo(hc),type="fan",tip.color=mypal[clus5], cex=.5, label.offset=.5) 24 plot(as.phylo(hc),type="fan",tip.color=mypal[clus5], cex=.5, label.offset=.5)
23 dev.off() 25 dev.off()
26 postscript('../data/results/hierarchical_cluster.eps', pointsize=12)
27 par(mar=c(1,1,1,1))
28 plot(as.phylo(hc),type="fan",tip.color=mypal[clus5], cex=.5, label.offset=.5)
29 dev.off()