-
Notifications
You must be signed in to change notification settings - Fork 1
/
plotNMDS.R
98 lines (82 loc) · 4.04 KB
/
plotNMDS.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
PlotNMDS = function(NMDS1, # results of the NMDS you want to plot
data, # Environmental data (not species matrix)
textp = T, # do you want to plot text for the species centroids?
group, # Group you want to put hulls around
group2= NA, # Group you want to put spiders around
taxa = NA) { # Vector of relative percentages of a given taxa from your species matrix, used to scale the points.
require(vegan)
#set up colors
cols = c("red", "blue", "green", "purple", "yellow",
"pink1", "brown", "orangered", "black", "cyan")
#set up sapes
shapes = c(1,2,3,4,5,6)
#set up group
names(data)[(which(names(data)== group))] = "group"
#number of groups
data = droplevels(data)
m = length(levels(data$group))
#set parameters so the legend is outside the plot area
par(xpd=NA,oma=c(3,0,0,0))
# plot the sample points
plot(NMDS1)
if (is.na(taxa[1])) points(NMDS1, pch=as.numeric(data$group), cex=1, col=cols[as.numeric(data$group)]) else {
# scale the points by relative abundance of a certain taxa
scalef = max(taxa) - min(taxa)
points(NMDS1, pch=as.numeric(data$group), cex=.1 + taxa/mean(taxa), col="blue")
}
# plot text for species, if desired
if(textp) text(NMDS1, dis="species", cex=.8)
# draw a hull around each group
for (i in 1:m) {
with(data, ordiellipse(NMDS1, kind = "sd", conf = 0.95, groups = group,
show = levels(group)[i], col=cols[i], lwd=2, lty = 1:length(group)))
}
#second group with spiders
if (is.na(group2)==F) {
names(data)[(which(names(data)== group2))] = "group2"
for (j in 1:length(levels(data$group2))) {
with(data, ordispider(NMDS1, group = group2, show = levels(group2)[j], col=cols[j + length(levels(group))], lwd=1))
}
legend(par("usr")[1],par("usr")[3], legend = c(levels(data$group),levels(data$group2)),
col = cols[1:(length(c(levels(data$group), levels(data$group2))))],
lwd=c(rep(2, length(levels(data$group))), rep(1, length(levels(data$group2)))))
} else legend(par("usr")[1],par("usr")[3], legend = levels(data$group),
col = cols[1:length(levels(data$group))],
lty = 1:length(levels(data$group)),
pch = shapes[1:length(levels(data$group))], lwd=2, bty = "o")
# legend(par("usr")[1]+1,par("usr")[3], legend = levels(data$group),
# pch = shapes[1:length(levels(data$group))], col = "blue", bty = "n")
}
PlotNMDS2 = function(NMDS1, # results of the NMDS you want to plot
data, # Environmental data (not species matrix)
textp = T, # do you want to plot text for the species centroids?
lines, #variables to use for ordisurf
group #group to change color/shape
) {
require(vegan)
#set up colors
cols = c("red", "blue", "green", "purple", "yellow",
"pink1", "brown", "orangered", "black", "cyan")
#set up sapes
shapes = c(1,2,3,4,5,6)
#set up group
names(data)[(which(names(data)== group))] = "group"
names(data)[(which(names(data)== lines))] = "lines"
#number of groups
data = droplevels(data)
#set parameters so the legend is outside the plot area
par(xpd=NA,oma=c(3,0,0,0))
# plot the sample points
plot(NMDS1, type="n", shrink = T, cex=1, xlim=c(-1, 1))
points(NMDS1, pch=as.numeric(data$group), cex=1, col=cols[as.numeric(data$group)])
# plot text for species, if desired
if(textp) text(NMDS1, dis="species", cex=.8)
#
with(data, ordisurf(NMDS1~lines, data = data, add=T, col = "black"))
legend(par("usr")[1],par("usr")[3], legend = levels(data$group),
col = cols[1:length(levels(data$group))],
lty = 1:length(levels(data$group)),
pch = shapes[1:length(levels(data$group))], lwd=2, bty = "o")
# legend(par("usr")[1]+1,par("usr")[3], legend = levels(data$group),
# pch = shapes[1:length(levels(data$group))], col = "blue", bty = "n")
}