nMDS in R vs PRIMER

Question

I used PRIMER v6 to plot my species data before following the steps shown in this video (https://www.youtube.com/watch?v=3I0zpl4v5wo&list=PL-W1dvsmEZ3KISwdig6X18dkQppJfcrYW&index=8) and get the plot below:

Now I'm trying to do the same thing using R with the codes that I copied from https://jkzorz.github.io/2019/06/06/NMDS.html:

dat <- Sp_Data[1:15,2:8]
dat <- dat %>% filter_all(any_vars(. != 0))
m_dat <- as.matrix(dat, labels = T)
mds <- metaMDS(m_dat, distance = "bray")

# ggplot
#extract NMDS scores (x and y coordinates)
data.scores = as.data.frame(scores(mds))
#add columns to data frame 
pt <- Sp_Data[which(rowSums(Sp_Data[1:15,2:8])>0),]
data.scores$Condition = pt$Condition
head(data.scores)

library(ggplot2)

xx = ggplot(data.scores, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 4, aes( shape = Condition))+ 
  theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"), 
        axis.text.x = element_text(colour = "black", face = "bold", size = 12), 
        legend.text = element_text(size = 12, face ="bold", colour ="black"), 
        legend.position = "right", axis.title.y = element_text(face = "bold", size = 14), 
        axis.title.x = element_text(face = "bold", size = 14, colour = "black"), 
        legend.title = element_text(size = 14, colour = "black", face = "bold"), 
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
        legend.key=element_blank()) + 
  labs(x = "NMDS1", y = "NMDS2", shape = "Condition")  + 
  scale_colour_manual(values = c("#009E73", "#E69F00")) 

xx

U1 is actually UL, and U2 is UR

My question is, how can the two plots be so different? The stress values are different too. Using R, it is 0.03435446, much lower than that using PRIMER

Here's the the compressed Sp_Data:

    structure(list(SiteConditionReplicate = c("JU1One", "JU1Two", 
"JU1Three", "JU1Four", "JU1Five", "JSOne", "JSTwo", "JSThree", 
"JSFour", "JSFive", "JU2One", "JU2Two", "JU2Three", "JU2Four", 
"JU2Five", "PU1One", "PU1Two", "PU1Three", "PU1Four", "PU1Five", 
"PSOne", "PSTwo", "PSThree", "PSFour", "PSFive", "PU2One", "PU2Two", 
"PU2Three", "PU2Four", "PU2Five"), `Camptandrium sexdentatum` = c(0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0), Diogenidae = c(0, 4, 0, 0, 1, 53, 0, 
2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
0, 0), `Nassarius globasus` = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
    `Tellina sp` = c(0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), `Marcia opima` = c(0, 
    0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Nereidae = c(0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 17, 18, 17, 12, 4, 3, 12, 
    1, 6, 0, 0, 0, 0, 0), Maldanidae = c(4, 5, 3, 2, 2, 1, 0, 
    0, 0, 0, 1, 0, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0), ...9 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA), Site = c("J", "J", "J", "J", "J", 
    "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "P", "P", 
    "P", "P", "P", "P", "P", "P", "P", "P", "P", "P", "P", "P", 
    "P"), Condition = c("U1", "U1", "U1", "U1", "U1", "S", "S", 
    "S", "S", "S", "U2", "U2", "U2", "U2", "U2", "U1", "U1", 
    "U1", "U1", "U1", "S", "S", "S", "S", "S", "U2", "U2", "U2", 
    "U2", "U2"), Replicate = c("One", "Two", "Three", "Four", 
    "Five", "One", "Two", "Three", "Four", "Five", "One", "Two", 
    "Three", "Four", "Five", "One", "Two", "Three", "Four", "Five", 
    "One", "Two", "Three", "Four", "Five", "One", "Two", "Three", 
    "Four", "Five")), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -30L))

Answer 1

This is not a complete solution, but it may get you on the right track. You can simplify your code a bit:

dat <- Sp_Data[1:15,2:8]
idx <- which(rowSums(dat)>0)   # Which rows are not all 0's
mds <- metaMDS(dat[idx, ], distance = "bray")  # No need to convert dat
code <- Sp_Data$Condition[idx]   # vector of Condition codes

This will run the analysis and save a vector of the Condition values. I can't follow your ggplot code since the first line throws an error:

data.scores = as.data.frame(scores(mds))
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  : 
  arguments imply differing number of rows: 13, 7

You are trying to convert a list to a data frame, but the number of rows is not the same as the number of columns. You probably want

data.scores <- scores(mds, tidy=TRUE)

There are still problems because Condition only applies to rows, not columns. Perhaps you do not want to plot the columns? In that case use

scores(mds)$sites