Packages used

library(mdatools)
library(vegan)
library(ggplot2)
library(plotly)
library(gapminder)
library(readr)
library(knitr)
library(rsconnect)

Import and organize data

phloem            <- read_csv("phloem_avg_no_outliers.csv")
phloem$Familytype <- as.factor(phloem$Familytype)
phloem            <- subset(phloem, Familytype=="Half")

Second-derivative transformation

expl <- as.numeric(phloem$BV)
spec <- phloem[, 6:262]

derv <- prep.savgol(as.matrix(spec), width=15, porder=2, dorder=2)
fam  <- as.factor(phloem$Family)

NMDS + output

dist      <- vegdist(derv, method = 'euclidian')
peak.NMDS <- metaMDS(dist, distance = "euclidian", k = 2, try = 1000)
## Run 0 stress 0.02573384 
## Run 1 stress 0.06446386 
## Run 2 stress 0.06630614 
## Run 3 stress 0.04577365 
## Run 4 stress 0.05060327 
## Run 5 stress 0.04763242 
## Run 6 stress 0.07126577 
## Run 7 stress 0.05757068 
## Run 8 stress 0.03845386 
## Run 9 stress 0.06767453 
## Run 10 stress 0.06917945 
## Run 11 stress 0.05792765 
## Run 12 stress 0.08011542 
## Run 13 stress 0.04264554 
## Run 14 stress 0.04694564 
## Run 15 stress 0.07157927 
## Run 16 stress 0.05385447 
## Run 17 stress 0.4153467 
## Run 18 stress 0.06215306 
## Run 19 stress 0.05896037 
## Run 20 stress 0.0606855 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     19: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
stressplot(peak.NMDS)

peak.NMDS$stress
## [1] 0.02573384
adonis2(dist ~ fam,  permutations = 9999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = dist ~ fam, permutations = 9999)
##           Df SumOfSqs      R2      F Pr(>F)
## fam       39   1.3072 0.27276 1.1059  0.302
## Residual 115   3.4852 0.72724              
## Total    154   4.7924 1.00000
## Significance indicates that the families are significantly grouped by distance

scores <- as.data.frame(scores(peak.NMDS))
all.scores   <- cbind.data.frame(expl, fam, scores)
colnames(all.scores) <- c("IBV", 'family', 'NMDS1', 'NMDS2')



centroids <- aggregate(cbind(NMDS1, NMDS2) ~ family*IBV, all.scores, mean)
f         <- function(z)sd(z)/sqrt(length(z))
se        <- aggregate(cbind(se.x=NMDS1, se.y=NMDS2) ~ family*IBV, all.scores, f)
centroids <- merge(centroids, se, by = "family")
centroids <- centroids[, c(1:4, 6,7)]
colnames(centroids) <- c("family","IBV",  "NMDS1", "NMDS2", "se.1", "se.2")
p <- ggplot(centroids, aes(NMDS1, NMDS2, label = family)) + geom_point(aes(col = IBV)) + theme_bw() +
  geom_point(data = centroids, aes(col=IBV)) +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=centroids, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=centroids, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = IBV)) +
  geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label = family), alpha=.5) +
  labs(title = "Half sibs, phloem") + scale_color_viridis_c() 
## Warning in geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label =
## family), : Ignoring unknown aesthetics: label
ggplotly(p) ## interactive plot

Since there’s a lot of overlap, I made the same figured zoomed in, which excludes some of the points but gives you a better look at the center.

p.zoom <- ggplot(centroids, aes(NMDS1, NMDS2, label = family)) + geom_point(aes(col = IBV)) + theme_bw() +
  geom_point(data = centroids, aes(col=IBV)) +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=centroids, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=centroids, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = IBV)) +
  geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label = family), alpha=.5) +
  labs(title = "Half sibs, phloem") + scale_color_viridis_c()  + ylim(-.11, .10) + xlim(-.27, .36)
## Warning in geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label =
## family), : Ignoring unknown aesthetics: label
ggplotly(p.zoom)

Also plotting just the average centroid with standard error for clarity.

p.cent <- ggplot(centroids, aes(NMDS1, NMDS2, label = family)) + geom_point(aes(col = IBV)) + theme_bw() +
  geom_point(data = centroids, aes(col=IBV)) +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=centroids, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=centroids, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = IBV)) +
  labs(title = "Half sibs, phloem") + scale_color_viridis_c() 

ggplotly(p.cent)