André Faust
Reputation: 31
How to summarize NMDS results into single points with error bars?
I have conducted an NMDS analysis and plotted the results. However, for the left plot, I would like to obtain a plot where the points from my various treatments are "summarized" into a single point with error bars, instead of what I have gotten. I provide images to better explain what I want plus my code just in case. The left plot in the "desired outcome image" was slapped together in paint, but I have seen these sort of plots in the literature (e.g.: Wild et al., 2014: "Input of easily available organic C and N stimulates microbial decomposition of soil organic matter in arctic permafrost soil"). Also they can be made with R (as stated in the paper), but I just cannot figure out how or find a tutorial/thread. TLDR: Instead of left plot in "my own results" I want left plot in "desired outcome". Anyone out there who can provide me with an answer?
Data:
structure(list(Treatment = c("cs", "cs", "cs", "cs", "ce", "ce",
"ce", "ce", "1e", "1e", "1e", "1e", "4e", "4e", "4e", "4e"),
`i15:0` = c(0.038194995, 0.070483554, 0.066509405, 0.02656837,
0.024501448, 0.123952048, 0.158812766, 0.048075503, 0.042707725,
0.173164856, 0.180104263, 0.162216335, 0.115837661, 0.321511027,
0.071114723, 0.194906795), `a15:0` = c(0.044211789, 0.068811108,
0.073337229, 0.027031903, 0.022587532, 0.110265989, 0.160241456,
0.04339058, 0.049255133, 0.170517668, 0.199264385, 0.154776239,
0.155124629, 0.307943496, 0.108737048, 0.194099006), `i16:0` = c(0.01782632,
0.031463009, 0.026320789, 0.014625177, 0.022190061, 0.051997237,
0.046475526, 0.021942647, 0.051267309, 0.093449378, 0.069918965,
0.079624038, 0.06122343, 0.146494324, 0.052256874, 0.087334431
), `16:1ω7c` = c(0.039181929, 0.064326029, 0.086634185,
0.038006212, 0.022568139, 0.111887339, 0.165744837, 0.067105897,
0.137095761, 0.251219335, 0.289299471, 0.276622201, 0.27988588,
0.359676703, 0.139182249, 0.422137325), `16:1ω7t` = c(0.006220504,
0.010372077, 0.011011986, 0.004331443, 0.004117218, 0, 0.016843949,
0.006185401, 0.01803384, 0.028295985, 0.02584535, 0.022196905,
0.031798788, 0.035570765, 0.018700339, 0.029134652), `16:1ω5c` = c(0.006957877,
0.013528688, 0.01707799, 0.00905064, 0.004203058, 0, 0.027683594,
0.019134836, 0.025887336, 0.057077469, 0.044940654, 0.075541206,
0.044327836, 0.07670695, 0.030417356, 0.127591822), `16:0` = c(0.161037157,
0.217587364, 0.183307479, 0.111910577, 0.120445222, 0.283660219,
0.333206677, 0.159581943, 0.381897382, 0.546141255, 0.379301617,
0.495703963, 0.595253076, 0.82653097, 0.315155249, 0.661990203
), `10Me 16:0` = c(0.008269657, 0.014867341, 0.013705167,
0.012176098, 0.00835949, 0, 0.024733271, 0.017993387, 0.035947535,
0.045106349, 0.042875244, 0.065672086, 0.045513819, 0.06571743,
0.032058128, 0.066599308), `i17:0` = c(0.019685743, 0.03821051,
0.027937903, 0.019666127, 0.020166071, 0.046030564, 0.039561999,
0.026039132, 0.065939957, 0.098216633, 0.05983684, 0.09703123,
0.07018488, 0.147983747, 0.052353567, 0.088069882), `a17:0` = c(0.015062966,
0.022975609, 0.019264877, 0.011306938, 0.021646964, 0, 0.031982809,
0.015408173, 0.063935976, 0.073056821, 0.05054158, 0.05557835,
0.052089068, 0.10369428, 0.043037189, 0.052491098), `17:1ω7c` = c(0.00428813,
0.005924667, 0.006807987, 0.003405714, 0.003100204, 0, 0.013317004,
0.00564397, 0.016137123, 0.025780419, 0.023742601, 0.023814462,
0.019288199, 0.028512027, 0.014184457, 0.028922394), `cy17:0` = c(0.007978745,
0.020837575, 0.01825151, 0.011241851, 0.00562496, 0.037428369,
0.037580648, 0.017905727, 0.048752828, 0.083793406, 0.074501831,
0.07364191, 0.063787883, 0.117775824, 0.018038811, 0.111461244
), `17:0` = c(0.019832995, 0.012288693, 0.016991794, 0.008906124,
0.015236496, 0.015629559, 0.033730574, 0.0123141, 0.047954072,
0.028504626, 0.036356188, 0.035751034, 0.057116806, 0.042311441,
0.030330047, 0.040258649), `10Me 17:0` = c(0.001967074, 0.005045711,
0.00419998, 0.004370348, 0.003049929, 0.008151985, 0.007665705,
0.005935455, 0.015951364, 0.016685053, 0.014753528, 0.026085769,
0.016103559, 0.024073947, 0.010236068, 0.027865828), `18:2ω6c` = c(0.006594482,
0.007138977, 0.011574679, 0.004210641, 0.005679057, 0.027457586,
0.017027548, 0.009400472, 0.030160413, 0.065461815, 0.0287921,
0.029084061, 0.03284687, 0.050610582, 0.022530081, 0.037026116
), `18:1ω7c` = c(0.063524628, 0.088011652, 0.064329236,
0.064397171, 0.054856424, 0.166204663, 0.135438015, 0.115515402,
0.243383684, 0.376587956, 0.191241645, 0.440923087, 0.268201012,
0.425619782, 0.157976459, 0.439641379), `18:1ω7t` = c(0.004339952,
0.009364455, 0.005604494, 0.003420479, 0.004732609, 0, 0.007931135,
0.004397696, 0.026273248, 0.022871272, 0.012048718, 0.009896217,
0.017956939, 0.024215907, 0.010083998, 0.009818137), `18:0` = c(0.152940658,
0.073439855, 0.092698797, 0.057528911, 0.122326908, 0.067357563,
0.206714901, 0.078515138, 0.316168351, 0.134563875, 0.195944637,
0.23662163, 0.496532244, 0.195335661, 0.195507236, 0.224913003
), `10Me 18:0` = c(0.003363204, 0.010884316, 0.00630025,
0.007155183, 0.004602433, 0, 0.009008464, 0.007156422, 0.017567226,
0.027625515, 0.014018314, 0.025328229, 0.014619648, 0.038015042,
0.019599232, 0.025250972), `cy19:0` = c(0.002473057, 0.012073369,
0.00509506, 0.005664806, 0.002999546, 0.002440478, 0.007511822,
0.008154366, 0.020524718, 0.043994752, 0.013101494, 0.02818864,
0.016691308, 0.049769211, 0.005023214, 0.02574527), `20:0` = c(0.129819646,
0.055626381, 0.08380711, 0.053866344, 0.134982677, 0.047650998,
0.13721842, 0.082087929, 0.196526612, 0.067974935, 0.135653207,
0.126493615, 0.296495625, 0.095695561, 0.153097109, 0.11146806
)), class = "data.frame", row.names = c(NA, -16L))
Data workaround (wa):
structure(list(Treatment = c("cs", "cs", "cs", "cs", "ce", "ce",
"ce", "ce", "1e", "1e", "1e", "1e", "4e", "4e", "4e", "4e", NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA)), class = "data.frame", row.names = c(NA, -37L
))
My R-code:
library(ggplot2)
library(ggvegan)
library(ggpubr)
#Load data and workaround for legend
data <- read.csv2("comSFdata.csv", header = TRUE, check.names = FALSE)
wa <- read.csv2("wa.csv", header = TRUE, check.names = FALSE)
#Create new dataframe w/o 1. row
com = data[,2:ncol(data)]
#Turn data from dataframe to matrix format
mcom = as.matrix(com)
#Run NMDS - set.seed for reproducibility
set.seed(123)
nmds1 <- metaMDS(mcom, distance = "bray", k = 2, maxit = 999, trymin = 500, trymax = 500)
#Prepare data for plotting
fort <- fortify(nmds1)
#Add workaround for custom legend
fort['Treatment'] <- wa['Treatment']
#Generate 2 plots:
p1 <- ggplot() +
#Customizes samples; alpha = opaqueness
geom_point(data = subset(fort, score == 'sites'),
mapping = aes(x = NMDS1, y = NMDS2, colour = Treatment),
alpha = 1) +
#Code axis range and force rendering of objects "outside" the range
coord_cartesian(ylim = c(-0.365, 0.365), xlim = c(-1.05, 1.05), clip = 'off') +
#Remove legend
theme(legend.position="none") +
#Add custom legend
annotate("text", x = 1.1, y = 0.38, size = 4, fontface = "bold", hjust = 1, colour = "#C77CFF",
label = ("Control start")) +
annotate("text", x = 1.1, y = 0.355, size = 4, fontface = "bold", hjust = 1, colour = "#00BFC4",
label = ("Control end")) +
annotate("text", x = 1.1, y = 0.33, size = 4, fontface = "bold", hjust = 1, colour = "#F8766D",
label = ("1% substrate addition")) +
annotate("text", x = 1.1, y = 0.305, size = 4, fontface = "bold", hjust = 1, colour = "#7CAE00",
label = ("4% substrate addition")) +
#Generates coordinate system axis at the origin
geom_abline(intercept = 0, slope = 0, linetype = "dashed", linewidth = 0.8, colour = "gray") +
geom_vline(aes(xintercept = 0), linetype = "dashed", linewidth = 0.8, colour = "gray") +
#Removes grids & background & adds a frame
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title below grid customization so it is in the foreground
annotate("text", x = -1.1, y = 0.375, size = 5.5, fontface = "bold", hjust = 0, colour = "black",
label = ("(a) Slump Floor (SF)")) +
#Customize left plot axis
theme(axis.ticks.length = unit(-3, "mm"),
axis.ticks = element_line(linewidth = 1.05),
axis.text.y.right = element_blank())
p2 <- ggplot() +
#Add custom x-axis label
geom_point(data = subset(fort, score == 'species'),
mapping = aes(x = NMDS1, y = NMDS2), alpha = 0) +
#Code axis range and force rendering of objects "outside" the range
coord_cartesian(ylim = c(-0.365, 0.365), xlim = c(-1.05, 1.05), clip = 'off') +
geom_segment(data = subset(fort, score == 'species'),
#Adds & customizes vectors for PLFAs
mapping = aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
arrow = arrow(length = unit(0, "npc"),
type = "closed"),
colour = "black", size = 0.8, alpha = 1)+
#Enable custom label colours and set to bold text
geom_text(data = subset(fort, score == 'species'),
mapping = aes(label = label, colour = label, x = NMDS1 * 1.1, y = NMDS2 * 1.1),
alpha = 1, size = 3, fontface = "bold") +
#Remove legend
theme(legend.position="none") +
#Add custom legend
annotate("text", x = 1.1, y = 0.38, size = 4, fontface = "bold", hjust = 1, colour = "black",
label = ("Non-specific")) +
annotate("text", x = 1.1, y = 0.355, size = 4, fontface = "bold", hjust = 1, colour = "#ff0000",
label = ("Gram-positive bacteria")) +
annotate("text", x = 1.1, y = 0.33, size = 4, fontface = "bold", hjust = 1, colour = "#ff8800",
label = ("Gram-negative bacteria")) +
annotate("text", x = 1.1, y = 0.305, size = 4, fontface = "bold", hjust = 1, colour = "#00b300",
label = ("Fungi")) +
#Customize label colours
scale_colour_manual(values=c("i15:0" = "#ff0000", "a15:0" = "#ff0000", "i16:0" = "#ff0000", "16:1ω7c" = "#ff8800", "16:1ω7t" = "#ff8800", "16:1ω5c" = "#ff8800", "16:0" = "black", "10Me 16:0" = "#ff0000", "i17:0" = "#ff0000", "a17:0" = "#ff0000", "17:1ω7c" = "#ff8800", "cy17:0" = "#ff8800", "17:0" = "black", "10Me 17:0" = "#ff0000", "18:2ω6c" = "#00b300", "18:1ω7c" = "#ff8800", "18:1ω7t" = "#ff8800", "18:0" = "black", "10Me 18:0" = "#ff0000", "cy19:0" = "#ff8800", "20:0" = "black")) +
geom_abline(intercept = 0, slope = 0, linetype = "dashed", linewidth = 0.8, colour = "gray") +
geom_vline(aes(xintercept = 0), linetype = "dashed", linewidth = 0.8, colour = "gray") +
#Removes grids & background & adds a frame
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title below grid customization so it is in the foreground
annotate("text", x = -1.1, y = 0.375, size = 5.5, fontface = "bold", hjust = 0, colour = "black",
label = ("(b) Slump Floor (SF)")) +
#Customize right plot axis
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.length = unit(-3, "mm"),
axis.ticks = element_line(linewidth = 1.05))
#Generate multi-panel plot with 2 columns & add/customize labels
ggarrange(p1, p2, widths = c(2.2, 2), ncol = 2, nrow = 1)```
Upvotes: 2
Views: 88
Answers (2)
André Faust
Reputation: 31
I found a solution thanks mainly to @Teodor (and a bit of ChatGPT). My code is quite similar but still somewhat different. Also, there is a lot of customization concerning my plots following the solution. As I oriented myself very closely to a graph from a scientific publication, I decided to post it in case someone wants the code for a (hopefully) publication-worthy looking plot (Fig. 1). I also attached the code for a simple version of the same graph but using ordibar to generate error bars that take covariance into account as proposed by @JariOksanen (Fig. 2).
Teodor's approach:
library(vegan)
library(dplyr)
library(ggplot2)
library(ggvegan)
library(ggrepel)
library(ggpubr)
#Load data
data <- read.csv2("comSFdata.csv", header = TRUE, check.names = FALSE)
#Create new dataframe w/o 1. column
com = data[,2:ncol(data)]
#Turn data from dataframe to matrix format
mcom = as.matrix(com)
#Run NMDS - set.seed for reproducibility
set.seed(123)
nmds1 <- metaMDS(mcom, distance = "bray", k = 2, trymax = 10000)
#Prepare data for, and calculate means + se --> Summarize into single points
nmds_coords <- as.data.frame(scores(nmds1, display = "sites"))
nmds_coords['Treatment'] <- data['Treatment']
nmds_summary <- nmds_coords %>%
group_by(Treatment) %>%
summarise(
NMDS1_mean = mean(NMDS1),
NMDS2_mean = mean(NMDS2),
NMDS1_sd = sd(NMDS1),
NMDS2_sd = sd(NMDS2),
NMDS1_se = NMDS1_sd / sqrt(n()),
NMDS2_se = NMDS2_sd / sqrt(n()))
#Prepare data for 2. plot
fort <- fortify(nmds1)
#Generate 2 plots:
p1 <- ggplot(data = nmds_summary, mapping = aes(x = NMDS1_mean, y = NMDS2_mean, colour = Treatment, shape = Treatment)) +
#Customizes samples & plot error bars
geom_errorbar(mapping = aes(xmin = NMDS1_mean - NMDS1_se, xmax = NMDS1_mean + NMDS1_se), width = 0.05, colour = "black", size = 1) +
geom_errorbar(mapping = aes(ymin = NMDS2_mean - NMDS2_se, ymax = NMDS2_mean + NMDS2_se), width = 0.14383561643835616438356164383562, colour = "black", size = 1) +
geom_point(size = 6, alpha = 1) +
scale_colour_manual(values = c("#C77CFF", "#00BFC4", "#F8766D", "#7CAE00")) +
scale_shape_manual(values = c(16, 16, 15, 17)) +
#Set axis range and force rendering of objects "outside" the range
coord_cartesian(ylim = c(-0.4, 0.4), xlim = c(-1.15, 1.15), clip = 'off') +
#Remove legend
theme(legend.position="none") +
#Add custom legend
annotate("text", x = 1.1, y = 0.365, size = 4, fontface = "bold", hjust = 1, colour = "#C77CFF",
label = ("Control start")) +
annotate("text", x = 1.1, y = 0.34, size = 4, fontface = "bold", hjust = 1, colour = "#00BFC4",
label = ("Control end")) +
annotate("text", x = 1.1, y = 0.315, size = 4, fontface = "bold", hjust = 1, colour = "#F8766D",
label = ("1% substrate addition")) +
annotate("text", x = 1.1, y = 0.29, size = 4, fontface = "bold", hjust = 1, colour = "#7CAE00",
label = ("4% substrate addition")) +
#Generates coordinate system axis at the origin
geom_abline(intercept = 0, slope = 0, linetype = "dotted", linewidth = 0.65, colour = "black") +
geom_vline(aes(xintercept = 0), linetype = "dotted", linewidth = 0.7125, colour = "black") +
#Removes grids & background & adds a frame
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title as label so it is in the foreground
annotate("label", x = -1.15, y = 0.36, size = 5.5, label.size = 0, fontface = "bold", hjust = 0, colour = "black",
label = ("(a) Slump Floor (SF)")) +
#Add custom axis titles
xlab("NMDS1") + ylab("NMDS2") +
#Customize left plot axis
theme(axis.ticks.length = unit(-3, "mm"),
axis.ticks = element_line(linewidth = 1.05),
axis.title.x = element_text(color="black", size=14, face="bold", vjust = -0.6),
axis.title.y.left = element_text(color="black", size=14, face="bold", vjust = 10),
#Set custom margin for better scale readability
plot.margin = margin(10, 10, 10, 35)) +
scale_x_continuous(expand = c(0, 0), breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""), sec.axis = sec_axis(~ ., breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""))) +
scale_y_continuous(expand = c(0, 0), breaks = c(-0.2, 0 , 0.2), labels = c("", "",""), sec.axis = sec_axis(~ ., breaks = c(-0.2, 0 , 0.2), labels = c("", "",""))) +
#Set custom x-axis labels
annotate("text", x = -0.95, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-1")) +
annotate("text", x = -0.43, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-0.5")) +
annotate("text", x = 0.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0")) +
annotate("text", x = 0.57, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0.5")) +
annotate("text", x = 1.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("1")) +
#Set custom y-axis labels
annotate("text", x = -1.2, y = 0.4, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0.4")) +
annotate("text", x = -1.2, y = 0.205, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0.2")) +
annotate("text", x = -1.2, y = 0, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0")) +
annotate("text", x = -1.2, y = -0.195, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-0.2")) +
annotate("text", x = -1.2, y = -0.395, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-0.4"))
p2 <- ggplot() +
#Add custom x-axis label
geom_point(data = subset(fort, score == 'species'),
mapping = aes(x = NMDS1, y = NMDS2), alpha = 0) +
#Set axis range and force rendering of objects "outside" the range
coord_cartesian(ylim = c(-0.4, 0.4), xlim = c(-1.15, 1.15), clip = 'off') +
geom_segment(data = subset(fort, score == 'species'),
#Adds & customizes vectors for PLFAs
mapping = aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
arrow = arrow(length = unit(0, "npc"), type = "closed"),
colour = "black", size = 0.8, alpha = 1)+
#Enable custom label colours and set to bold text
geom_text_repel(data = subset(fort, score == 'species'),
mapping = aes(label = label, colour = label, x = NMDS1 * 1.025, y = NMDS2 * 1.025),
alpha = 1, size = 3, fontface = "bold") +
#Remove legend
theme(legend.position="none") +
#Add custom legend
annotate("text", x = 1.1, y = 0.365, size = 4, fontface = "bold", hjust = 1, colour = "black",
label = ("Non-specific")) +
annotate("text", x = 1.1, y = 0.34, size = 4, fontface = "bold", hjust = 1, colour = "#ff0000",
label = ("Gram-positive bacteria")) +
annotate("text", x = 1.1, y = 0.315, size = 4, fontface = "bold", hjust = 1, colour = "#ff8800",
label = ("Gram-negative bacteria")) +
annotate("text", x = 1.1, y = 0.29, size = 4, fontface = "bold", hjust = 1, colour = "#00b300",
label = ("Fungi")) +
#Customize label colours
scale_colour_manual(values=c("i15:0" = "#ff0000", "a15:0" = "#ff0000", "i16:0" = "#ff0000", "16:1ω7c" = "#ff8800", "16:1ω7t" = "#ff8800", "16:1ω5c" = "#ff8800", "16:0" = "black", "10Me 16:0" = "#ff0000", "i17:0" = "#ff0000", "a17:0" = "#ff0000", "17:1ω7c" = "#ff8800", "cy17:0" = "#ff8800", "17:0" = "black", "10Me 17:0" = "#ff0000", "18:2ω6c" = "#00b300", "18:1ω7c" = "#ff8800", "18:1ω7t" = "#ff8800", "18:0" = "black", "10Me 18:0" = "#ff0000", "cy19:0" = "#ff8800", "20:0" = "black")) +
#Generates coordinate system axis at the origin
geom_abline(intercept = 0, slope = 0, linetype = "dotted", linewidth = 0.65, colour = "black") +
geom_vline(aes(xintercept = 0), linetype = "dotted", linewidth = 0.7125, colour = "black") +
#Removes grids & background & adds a frame
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title below grid customization so it is in the foreground
annotate("label", x = -1.15, y = 0.36, size = 5.5, label.size = 0, fontface = "bold", hjust = 0, colour = "black",
label = ("(b) Slump Floor (SF)")) +
#Customize right plot axis
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.length = unit(-3, "mm"),
axis.ticks = element_line(linewidth = 1.05),
axis.title.x = element_text(color="black", size=14, face="bold", vjust = -0.6),
#Set custom margin for better scale readability
plot.margin = margin(10, 10, 10, 10)) +
scale_x_continuous(expand = c(0, 0), breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""), sec.axis = sec_axis(~ ., breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""))) +
scale_y_continuous(expand = c(0, 0), breaks = c(-0.2, 0 , 0.2), labels = c("", "",""), sec.axis = sec_axis(~ ., breaks = c(-0.2, 0 , 0.2), labels = c("", "",""))) +
#Set custom x-axis labels
annotate("text", x = -0.95, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-1")) +
annotate("text", x = -0.43, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("-0.5")) +
annotate("text", x = 0.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0")) +
annotate("text", x = 0.57, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("0.5")) +
annotate("text", x = 1.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
label = ("1"))
#Generate multi-panel plot with 2 columns & add/customize labels
ggarrange(p1, p2, ncol = 2, nrow = 1, widths = c(1.06, 1), align = "h")
JariO's approach:
library(vegan)
library(dplyr)
library(ggplot2)
library(ggvegan)
library(ggrepel)
library(ggpubr)
#Load data
data <- read.csv2("comSFdata.csv", header = TRUE, check.names = FALSE)
#Create new dataframe w/o 1. column
com = data[,2:ncol(data)]
#Turn data from dataframe to matrix format
mcom = as.matrix(com)
#Run NMDS - set.seed for reproducibility
set.seed(123)
nmds1 <- metaMDS(mcom, distance = "bray", k = 2, trymax = 10000)
#Prepare data for plotting
fort <- fortify(nmds1)
#Add workaround for adding back Treatment groups
fort['Treatment'] <- wa['Treatment']
#Remove rows without Treatment info
newfort <- na.omit(fort)
label <- c("#C77CFF","#C77CFF","#C77CFF","#C77CFF","#00BFC4","#00BFC4","#00BFC4","#00BFC4","#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00")
label2 <- c("#C77CFF","#00BFC4","#F8766D","#7CAE00")
#Ordibar
ordiplot(nmds1, display = "sites", type = "none", main = "Ordibar NMDS")
ordibar(nmds1, display = "sites", groups = data$Treatment, col = label2, length = 0.2)
points(newfort[, 3], newfort[, 4], pch = 16, cex = 1.5, col = label)
ordiellipse(nmds1, groups = data$Treatment, draw = "lines", col = label2, alpha = 1)
Upvotes: 1
Teodor
Reputation: 38
Is this close to what you need?
The code:
# Load packages
library(tidyverse)
library(ggplot2)
# Make some fake NMDS results:
sampsize <- 10
treatments <- c("start", "end", "1%", "4%")
treatment <- rep(treatments, each = sampsize)
nmds1means <- c(-0.3, 0, 0.3, 0.6)
nmds2means <- c(0.2, -0.7, -0.1, -0.4)
nmds1 <- sapply(nmds1means, FUN = function(x){rnorm(n = sampsize, mean = x, sd = 0.12)}) %>% as.vector
nmds2<- sapply(nmds2means, FUN = function(x){rnorm(n = sampsize, mean = x, sd = 0.2)}) %>% as.vector
df <- data.frame(treatment, nmds1, nmds2)
# Plot as "normal" NMDS plot:
ggplot(data = df, mapping = aes(x = nmds1, y = nmds2)) +
geom_point()
# Prepare data for summarised NMDS plot:
df2 <- df %>%
dplyr::group_by(treatment) %>%
dplyr::summarise(nmds1_mean = mean(nmds1),
nmds1_se = sd(nmds1),
nmds2_mean = mean(nmds2),
nmds2_se = sd(nmds2))
# Plot summarised NMDS plot:
ggplot(data = df2, mapping = aes(x = nmds1_mean, y = nmds2_mean, shape = treatment)) +
geom_point(size = 4) +
scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_errorbar(mapping = aes(ymin = nmds2_mean - nmds2_se, ymax = nmds2_mean + nmds2_se)) +
geom_errorbar(mapping = aes(xmin = nmds1_mean - nmds1_se, xmax = nmds1_mean + nmds1_se)) +
lims(x = c(-1, 1), y = c(-1, 1)) +
xlab("NMDS1") +
ylab("NMDS2")
Upvotes: 1
Related Questions
- R & ggplot2: How to get arrows under the axis label?
- How does one control the appearance (e.g. line size, line type, colour) of mqgam plots produced using plot.mgamViz from the "mgcViz" package?
- Change the size of X/ Y labels in a ggplot with unusual arguments
- Add a legend for multiple geom_segments while maintaining a "solid" linetype in ggplot2
- Subscript a title in a Graph (ggplot2) with label of another file
- Removing axis labelling for one geom when multiple geoms are present
- Why does geom_text() throw coercion errors when hjust and vjust are strings?
- ggplot2 polar plot axis label location