library(magrittr)
library(dplyr)
library(ggpubr)
library(MASS)
library(vegan)
library(foreign)
library(reshape2)
library(irlba)
library(mclust)
library(nnet)
library(readxl)
###data import###
setwd("insert path here")
abundance=read_excel("abundance data here", sheet=1)
abundance=cbind(abundance, freq=1)
habitat=read_excel("climate data here", sheet=1)
###user's choice###
Gmax=30  #maximum number of mixture components
NVmax=20 #maximum number of principal components fro the sparce abundance matrix
###data handling (long to wide format)
data=dcast(abundance, Plots~Species,value.var="freq")
data=as.matrix(data)
data[is.na(data)]=0
###data projection onto a low-dimensional space###
data.svd=irlba(data, nv=NVmax)
data.rid=data%*%data.svd$v
###finite mixture model###
data.clus=Mclust(data.rid, G=1:Gmax)
data.rid=cbind(data.rid, data.clus$classification, data.clus$z)
data.rid=as.data.frame(data.rid)
colnames(data.rid)=c(paste("Comp", 1:NVmax,sep=""),"classification", paste("Post", 1:dim(data.clus$z)[2],sep="")) 
idx1=order(data[,1])
idx2=order(habitat$Plot)
data.rid=data.rid[idx1,]
habitat=habitat[idx2,]
data.all=cbind(data.rid, habitat)
data.log=multinom(classification~COVARIATES, data=data.all, weights=data.clus$z[data.clus$classification])
###mixture discriminant analysis###
data.mda=mda(classification~COVARIATES, data=data.all)
###non-metric MDS###
cluster.means=apply(data.all[,1:NVmax], 2, function(z) tapply(z, data.all$classification, mean))
distance=dist(cluster.means)
dist=as.matrix(distance)
nmds=isoMDS(dist, k=2)
mds=cmdscale(dist, k=2)
colnames(nmds$points)=c("Dim.1", "Dim.2")
rownames(nmds$points)=c(paste("Group", 1:24,sep="")) 
# Plot NMDS
ggscatter(data.frame(nmds$points), x = "Dim.1", y = "Dim.2", label = rownames(nmds$points), size = 1, repel = TRUE)

colnames(mds)=c("Dim.1", "Dim.2")
rownames(mds)=c(paste("Group", 1:24,sep="")) 
# Plot MDS
ggscatter(data.frame(mds), x = "Dim.1", y = "Dim.2", label = rownames(nmds$points), size = 1, repel = TRUE)