3D nmds

data(dune)
dis <- vegdist(dune)
m <- monoMDS(dis, model = "loc", k = 3)
m
## 
## Call:
## monoMDS(dist = dis, k = 3, model = "loc") 
## 
## Local non-metric Multidimensional Scaling
## 
## 20 points, dissimilarity 'bray', call 'vegdist(x = dune)'
## 
## Dimensions: 3 
## Stress:     0.04788406 
## Stress type 1, weak ties
## Scores scaled to unit root mean square, rotated to principal components
## Stopped after 136 iterations: Stress nearly unchanged (ratio > sratmax)
data("dune.env")
Data <- cbind(dune.env, m$points)
plot_ly(
  data = Data, x = ~MDS1, y = ~MDS2, z = ~MDS3, 
  color = ~Use, colors = viridis::viridis(n = length(unique(Data$Use)))
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'Community'))
        )