##--------------------------------------------------------------------------------------------------------------
## R-Scripts used for analysis in part II ......."population genetics and clonal richness!!""
#---------------------------------------------------------------------------------------------------------------
# Title: Environment not dispersal limitation drives clonal composition of arctic Daphnia 
# in a recently deglaciated area 
#-------------------------------------------------------------------------------
CR<-subset(haile,select=cbind(1:42))  			                          ## data frame of all 42 clones

#-------------------------------------------------------------------------------
# Data frame of (9) Environmental variables( TP, TN, Chl-a, Cond, pH, Temp, Alt and Dep)
ENV<- subset(haile, select = cbind (43:51))	                            
ENV
dim(ENV)

#-------------------------------------------------------------------------------
## Data frame of Geographic location with UTM reading collected in the field
SPACE<- subset(haile, select = cbind (53:54))		                          
SPACE
dim(SPACE)

## Data of distance (in km) from the Ice sheet terminus in the SE direction
Dist<- subset(haile, select = cbind (distance1))	                      
Dist

## Data of distance (in km) from tongue of Ice sheet terminal in NE direction
Dist2<- subset(haile, select = cbind (distance3))	                      
Dist2

#-------------------------------------------------------------------------------
age1<- subset(haile, select = cbind (Young,Old))		                   ## age of habitat (clusters)
age1

#----------------------------------------------------------------------------------------------------------------------
age2<- subset(haile, select = cbind (distance1,distance3))		          
## age of habitat (clustered with the two distance measures) another way to express age gradient with idea
age2

age3<- subset(haile, select = cbind (Young,Old,distance1,distance3))		## age of habitat (clusters)
age3

#-------------------------------------------------------------------------------
hab<- subset(haile, select = cbind (Lake,Pond))		                     ## habitat (Pond or Lake)
hab
head(hab)                                                                    ## Display only the first few lines

#------------------------------------------------------------------------------ 
# Figure 3. Clonal composition in lakes (A) and ponds (B).
#------------------------------------------------------------------------------ 
# Map of All the Lakes and ponds
 x11(width=22, height=22)
plot(SPACE,asp=1,type="p", col="blue",main="Geographic Coordinates of the Study Sites",
xlab="X coordinates(m)",ylab="Y coordinates(m)")
text(546900,7666950, "Young cluster", cex=1.2, col="Black")
text(544000,7660000, "Old cluster", cex=1.2, col="Black")
text(SPACE, row.names(SPACE), cex=0.4, col="blue",pos=1)


#-----------------------------------------------------------------------------------------------------------------------
##----------------------------------------------------------------------------------------------------------------------
CR<-subset(haile,select=cbind(1:42))  			## All 42 MLGs
CR
# Factor out environmental variables (ENV) which are 
# 9 variables measured for each site( i.e TP, TN, Chl-a, Cond, pH, Temp, Alt and Dep)
env<- subset(haile, select = cbind (43:51))	      
env

# Factor out the geographic location with UTM reading collected in the field (SPACE)
spa<- subset(haile, select = cbind (53:54))		
dim(spa)
#spatial predictors from the X-Y reading centred X-Y coordinates
spa.c <- scale(spa, center=TRUE, scale=FALSE)     
spa.c

#-------------------------------------------------------------------------------
# Factor out habitat age clusters based on dummy variables and distances to icesheet and
age<- subset(haile, select = cbind (Young,Old))		
age

## Factor out habitat age clusters (continous)
age2<- subset(haile, select = cbind (distance1,distance3))		
age2

# Factor out habitat age cluster dummy varaibles and surrogate as AGE
age3<- subset(haile, select = cbind (Young,Old,distance1,distance3))		
age3

#-------------------------------------------------------------------------------
## Factor out habitat type (Pond or Lake) using dummy  varaibles
hab<- subset(haile, select = cbind (Lake,Pond))	
hab
head(hab) 								
#--------------------------------------------------------------------------------------------------------------#
# Transform the clonal abaundance data using Hellinger transformation
# to create a hellinger transformed data matrix use vegan package
# Load the package Vegan
library(vegan)
spe.hel <- decostand (CR, "hellinger")        	  
spe.hel

#-------------------------------------------------------------------------------#-----------------------------------#
##*******************************************************************************************************
# Table 2 : 
# Global efffect A: effect of Environmental variables with out taking space into account tested in RDA
#  RDA of Hellinger-transformed clonal composition data (spe.hel) constrained 
#  by all the environmental variables (env) reported in Table 2  and used for Fig 4 too
#---------------------------------------------------------------------------#---------------------------------------#
#--------------------------------------------------------------------------------------------------------------------
## RDA with all explanatory environmental variables
spe.rda.all <- rda(spe.hel ~ ., data=env)	            
spe.rda.all

### test if the relationship between environement and spe abu is significant and calculate ##the global adjusted R^2
RsquareAdj(spe.rda.all)
(R2a.all <- RsquareAdj(spe.rda.all)$adj.r.squared) 
anova(spe.rda.all,step=1000)   
plot(spe.rda.all)   	

#******************************************************************************************************
# using Forward selection in R package  packfor's forward.sel()   
# Only depth, Cond, Tp and TN were found to be significant
library(packfor)
forward.sel(spe.hel, env, adjR2thresh=R2a.all)		
(step.forward <- ordistep(rda(spe.hel ~ 1, data=env),scope=formula(spe.rda.all ), direction="forward", pstep=1000))
summary(step.forward)					
# Forward selection using a double stopping criterion (Blanchet et al. 2008a)
(step.forward2 <- ordiR2step(rda(spe.hel ~ 1, data=env),scope=formula(spe.rda.all ), direction="forward", pstep=1000))
summary(step.forward2)						

#------------------------------------------------------------------------------------------------------------------
# Fowrward selection (FStotal) results of the environmental variables for Table 2
(rda.envspe <- rda(spe.hel ~ ., env))			  			## 
anova(rda.envspe,step=1000)
(R2a.rda.envspe<- RsquareAdj(rda.envspe)$adj.r.squared)
(envspe.fwd <- forward.sel(spe.hel,env, adjR2thresh=R2a.rda.envspe))
envspe.sign <- sort(envspe.fwd[,1])				                   
envspe.sign<- env[,c(envspe.sign)]
envspe.sign                            ##"Cond" "Dep" "TN"  "TP", are  significant environmental varaibles

#------------------------------------------------------------------------------------------------------------------
rda.envspe.forw=rda(spe.hel,envspe.sign)
summary(rda.envspe.forw)
anova(rda.envspe.forw,step=1000)
RsquareAdj(rda.envspe.forw)


#-----------------------------------------------------------------------------------------------------------------#
#******************************************************************************************************
# Global efffect A: spatial descriptors(variables (SPACE) and 
# PCNM Trend-surface analysis of clonal data
# This shows the the linear combination of spatial variables, in a projection preserving 
#------------------------------------------------------------------------------------------------------------------#
# principal coordinates of neighbor matrices(PCNM)
spa.d<-dist(spa)
spanning<-spantree(spa.d)

#chosing the threshold as the maximum distance between two points
dmin=max(spanning$dist)
spa.d[spa.d>dmin]<- 4*dmin					                   	      	  #truncate the distance matrix
spa.PCoA=cmdscale(spa.d, k=nrow(spa)-1, eig=T)	        	      	#PCoA of truncated distance matrix
nb.ev=length(which(spa.PCoA$eig>0.0000001))				                #Count the positive eigenvalues 
spe.PCNM=as.data.frame(spa.PCoA$points[1:nrow(spa), 1:nb.ev])	    #construct a data frame with PCNM variables
summary(spe.PCNM)

#--------------------------------------------------------------------------------------------------------------------------------
# calculate Moran's I for each Eigenvectors for the threshold distance and 
# Part one Construct the PCNM variables automatically using PCNM()
library(PCNM)
spe.PCNM.auto=PCNM(spa.d)
summary(spe.PCNM.auto)
source("plot.spantree.R")
plot.spantree(spe.PCNM.auto$spanning, spa)				    #plot the minimum spanning tree used to find the truncation distance
plot(spe.PCNM.auto$spanning, spa)
dmin= spe.PCNM.auto$thresh				                    ###truncation distance
nb.ev=length(spe.PCNM.auto$values) 			              ###nbr of eigenvalues
# Moran's I of the PCNM variables (in the first distance class 0 to truncation threshold).
spe.PCNM.auto$expected_Moran						              # Expected value of I, no spatial correlation
spe.PCNM.auto$Moran_I								                  #Moran's I for each PCNM variables
select=which(spe.PCNM.auto$Moran_I$Positive == T)			#Eigenfunction with positive spatial correlation
spe.PCNM.pos=as.data.frame(spe.PCNM.auto$vectors)[,select]
spe.PCNM.pos

#--------------------------------------------------------------------------------------------------------------------------------
# part two. Run the global PCNM analysis on the  data
(spe.PCNM.rda <- rda(spe.hel, spe.PCNM.pos))
anova(spe.PCNM.rda, step=1000)
RsquareAdj(spe.PCNM.rda)

## 3. Since the analysis is significant, compute the adjusted R2 and and run a forward selection of the PCNM variables
(spe.R2a <- RsquareAdj(spe.PCNM.rda)$adj.r.squared)
spe.R2a 
(spe.PCNM.fwd <- forward.sel(spe.hel, as.matrix(spe.PCNM.pos), adjR2thresh=spe.R2a))
spe.PCNM.fwd

# According to the R2adj criterion, if we retain 2 PCNMs (dbMEMs), we get a model with 
(nb.sig.PCNM <- nrow(spe.PCNM.fwd))	# Number of signif. PCNM
# Identity of significant PCNMs in increasing order
(PCNM.sign <- sort(spe.PCNM.fwd[,2]))
# Write the significant PCNMs to a new object
PCNM.read <- spe.PCNM.pos[,c(PCNM.sign)]
PCNM.read

# New RDA using the 2 spatial distance reetained
PCNM.read.rda2=rda(spe.hel, PCNM.read)
summary(PCNM.read.rda2)
RsquareAdj(PCNM.read.rda2)

#overall test and test of canonical axes
anova(PCNM.read.rda2, step=1000)

#------------------------------------------------------------------------------------------------------------------
# Table 2 partial RDA analysis[Variation  partitioning]between the environmental and spatial matrices 
# with out age as independent variable.. Figure S1 A
#------------------------------------------------------------------------------------------------------------------
# In a partial RDA analysis variation partitioning Used to produce Figure S1
# total variation explained(21.1%) by explanatory variables
# unclaimed variation given in residuals
(spe.varpar<-varpart(spe.hel,PCNM.read,envspe.sign))		
plot(spe.varpar,digit=1)				                      	
spe.varpar
plot(spe.varpar,digit=1)

#------------------------------------------------------------------------------------------------------------------
# Unique effect of envspe.sign...... {only the signifcant variables selected by forward selection}
rda.envspe<-rda(spe.hel,envspe.sign,cbind(PCNM.read))
anova(rda.envspe,step=10000,perm.max=10000,model="reduced")
RsquareAdj(rda.envspe)$adj.r.squared

#------------------------------------------------------------------------------------------------------------------
# Unique effect of spa......only the significant dbMEM seceted by the forward selection}
rda.spaspe<-rda(spe.hel,PCNM.read,cbind(envspe.sign))
anova(rda.spaspe,step=10000,perm.max=10000,model="reduced")
RsquareAdj(rda.spaspe)$adj.r.squared

#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
# Figure 5. clone abundance - environmental variables retained by forward selection - age3 - PCNM variation partitioning
#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
# when AGE is included in the explanation
varpar.spe11<-varpart(spe.hel,envspe.sign,PCNM.read, age3)		## total variation explained
varpar.spe11
plot(varpar.spe11,digit=1) 						#### unclaimed variation given in residuals

# Fraction [a], pure PCNM.read
anova(rda(spe.hel, PCNM.read, cbind(age3, envspe.sign)))

# Fraction [b], pure environmental
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read, age3)))

# Fraction [c], age3
anova(rda(spe.hel, age3, cbind(PCNM.read, envspe.sign)))

#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
# when age is repersented by age cluster dummy varaibles and distance as surrogate messure is indluded
# as presented in FIG.5 of the main text.....TABLE s4 of the suplementary information
(spe.varpart03 <- varpart(spe.hel, PCNM.read,envspe.sign, age3))
par(mfrow=c(1,2))
plot(spe.varpart03, digits=1)

# Total effect of all the three predictors together # (ENV, SPACE and AGE)
rda.globalspe<-rda(spe.hel,cbind(PCNM.read, envspe.sign, age3))
anova(rda.globalspe,step=10000,perm.max=10000,model="reduced")
R2a.rda.globalspe=RsquareAdj(rda.globalspe)$adj.r.squared
R2a.rda.globalspe 

# Tests of the unique fractions [a], [b], [c] and [d] and other marginal effects
# Fraction [a], pure environmental effects
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read, age3)))

# Fraction [b], pure PCNM.read space effects
anova(rda(spe.hel, PCNM.read, cbind(age3, envspe.sign)))

# Fraction [c], pure (age3)  age effects
anova(rda(spe.hel, age3, cbind(PCNM.read, envspe.sign)))


#------------------------------------------------------------------------------------------------------------------
# Table S5 and Figure S2. clone - environment significant only - age3 - PCNM variation partitioning and hab
#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
# to partition the variation in the species matrix between ENV, SPACE, AGE and Hab
# [ age using distance as surrogate measure and hab is included ]
varparnew<-varpart(spe.hel,PCNM.read,envspe.sign, age3,hab)		## total variation explained
windows(title="clone - environment - PCNM variation partitioning", 8, 6)
plot(varparnew,digit=2) 						#### unclaimed variation given in residuals
varparnew

# Total effect of all the 4 predictors together (env, space , age and hab)
rda.globalspe2<-rda(spe.hel,cbind(PCNM.read, envspe.sign, age3,hab))
anova(rda.globalspe2,step=10000,perm.max=10000,model="reduced")
R2a.rda.globalspe2=RsquareAdj(rda.globalspe2)$adj.r.squared
R2a.rda.globalspe2 

# Tests of the unique fractions [a], [b], [c] and [d]

# Fraction [a], pure environmental
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read, age3,hab)))

# Fraction [b], pure PCNM.read
anova(rda(spe.hel, PCNM.read, cbind(age3, envspe.sign,hab)))

# Fraction [c], age3
anova(rda(spe.hel, age3, cbind(PCNM.read, envspe.sign,hab)))

# Fraction [c], hab
anova(rda(spe.hel, hab, cbind(PCNM.read, envspe.sign,age3)))

# Fraction [c], SPACE | AGE+HAB        
anova(rda(spe.hel, PCNM.read, cbind(hab,age3)))

#SPACE | ENV+HAB 
anova(rda(spe.hel, PCNM.read, cbind(envspe.sign,hab)))

#SPACE | ENV+AGE 
anova(rda(spe.hel, PCNM.read, cbind(envspe.sign,age3)))

#ENV | AGE+HAB 
anova(rda(spe.hel, envspe.sign, cbind(hab,age3)))

#ENV | SPACE+HAB
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read, hab)))

#ENV | SPACE+AGE  
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read,age3)))

#AGE | SPACE+HAB 
anova(rda(spe.hel, age3, cbind(PCNM.read, envspe.sign)))

#AGE | ENV+HAB 
anova(rda(spe.hel, age3, cbind(hab, envspe.sign)))

#AGE | SPACE+ENV
anova(rda(spe.hel, age3, cbind(PCNM.read,hab)))

#HAB | ENV+AGE 
anova(rda(spe.hel, hab, cbind(envspe.sign,age3)))

#HAB | SPACE+AGE 
anova(rda(spe.hel, hab, cbind(PCNM.read, age3)))

#HAB | SPACE+ENV 
anova(rda(spe.hel, hab, cbind(PCNM.read, envspe.sign)))

#HAB | SPACE
anova(rda(spe.hel, hab, cbind(PCNM.read)))
#AGE | SPACE+ENV
anova(rda(spe.hel, hab, cbind(envspe.sign)))
#AGE | SPACE+ENV
anova(rda(spe.hel, hab, cbind(age3)))

#AGE | SPACE+ENV
anova(rda(spe.hel, PCNM.read, cbind(age3)))
#AGE | SPACE+ENV
anova(rda(spe.hel, PCNM.read, cbind(envspe.sign)))
#AGE | SPACE+ENV
anova(rda(spe.hel, PCNM.read, cbind(hab)))

#AGE | SPACE+ENV
anova(rda(spe.hel, age3, cbind(hab)))
#AGE | SPACE+ENV
anova(rda(spe.hel, age3, cbind(envspe.sign)))
#AGE | SPACE+ENV
anova(rda(spe.hel, age3, cbind(PCNM.read)))

#AGE | SPACE+ENV
anova(rda(spe.hel, envspe.sign, cbind(hab)))
#AGE | SPACE+ENV
anova(rda(spe.hel, envspe.sign, cbind(age3)))
#AGE | SPACE+ENV
anova(rda(spe.hel, envspe.sign, cbind(PCNM.read)))

#AGE | SPACE+ENV
anova(rda(spe.hel, hab, cbind(PCNM.read, envspe.sign,age3)))

(spe.varpar0<-varpart(spe.hel,PCNM.read,envspe.sign,age3,hab))		        
## total variation explained(23.1%) by explanatory variables
plot(spe.varpar0,digit=2)					               
#### unclaimed variation given in residuals
spe.varpar0

## Use fill colours
showvarparts(4, bg = c("hotpink","skyblue","yellow","green"))
plot(spe.varpar0,digit=2, bg = c("hotpink","skyblue","yellow","green"))

#------------------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------------------
#*****************************************************************************************************************
#*****************************************************************************************************************
##  Part Two: Molecular ecology........_R Analysis of  Paper two: spatial population genetics in greenland 
#  To start polySat...Locate where the file is "polysat-paper-two.txt" is and 
# Load polysat package Clark, LV and Jasieniuk, M, 2011  Molecular Ecology Resources 11(3): 562-566.
# ****************************************************************************************************************		
# This is the analsyis and data that was used to asses the genetic relatedness figure 3a and b
# The objective of the analysis
# 1. To assess the genetic relatedness among different MLGs identified from different population
# packages used are (vegan")and (polysat) and adegenet
library(vegan)
library(polysat)  
library(adegenet)
Clones_01<- read.GeneMapper("Polysat-paper-two-all-clones.txt") 	### read data of 9 loci of 4 alleles and 1420 individuals from 61 population
summary(Clones_01)							### summary of data
Samples(Clones_01)							### sample names
Loci(Clones_01)	
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Description(Clones_01) <- "Dataset for 1420 Individuals from Lakes and ponds grouped by population (n=61) and clones (n=42)"
PopNames(Clones_01) <-c("M1","M2","M3","M4","M5","M6","M7","M8","M9","M10","M11","M12","M13",
                        "M14","M15","M16","M17","M18","M19","M20","M21","M22","M23","M24","M25",
                        "M26","M27","M28","M29","M30","M31","M32","M33","M34","M35","M36","M37",
                        "M38","M39","M40","M41","M42")
PopInfo(Clones_01)<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,
                      26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42)

# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Usatnts(Clones_01) <- c(2, 2, 2, 2, 2, 3, 2, 3, 2)
PopInfo(Clones_01)
Usatnts(Clones_01)

# To edit genotypes after importing the data					
Clones_011<- editGenotypes(Clones_01, maxalleles = 4) 				#### to calculate polyploid level of multilocus genotype
Clones_011 <- estimatePloidy(Clones_011)
summary(Clones_011)
PopInfo(Clones_011) 
#  Part one: Individual-level population genetics statistic analysis [Genetic distances between individuals]
geneticDBR <- meandistance.matrix(Clones_011) #.........1.Genetic Distance matrix of Bruvo et al. (2004)
geneticDBR 

# PCA of genetic distance and Ward Hierarchical Clustering distance matrix
#------------------------------------------------------------------------------------------------------------------
pca<- cmdscale(geneticDBR)
pca
plot(pca[,1], pca[,2],
     main="Bruvo distance with meandistance.matrix",xlab="PCA1",ylab="PCA2")

D <- dist(pca[,1:2])^2
clust <- hclust(D, method="ward.D2")
plot(clust, hang = -1,main="dendrogram based on ward clustering based on the first 2 PCs", cex=1)
pop <- factor(cutree(clust, k=4))
pop
sort(pop)
library(adegenet)
x11(width=22, height=22)
s.class(pca, fac=pop , col=transp(funky(5)), cpoint=2,sub="PCA - axes 1 and 2")
mypal = c("red", "blue", "green", "black")				# vector of colors
s.class(pca, fac=pop, col=mypal, cpoint=2,sub="PCA - axes 1 and 2")

x11(width=22, height=22)
plot(clust,ylab="Bruvo's genetic distance")

# labels at the same level
plot(clust, hang =-1, ylab="Bruvo's genetic distance",xlab="clone",cex=.8)

# using dendrogram objects
hcd = as.dendrogram(clust)
# alternative way to get a dendrogram
plot(hcd)



 #------------------------------------------------------------------------------------------------------------------
 #------------------------------------------------------------------------------------------------------------------
 ## The indicator species analysis, the Indicator Value index of R
 #------------------------------------------------------------------------------------------------------------------
 #------------------------------------------------------------------------------------------------------------------
 # Step one read the data and calculate diversity partitioining
 Masterdata<-read.table("species richness-paper-one-new-02.txt",header=T,row.names=1)
 # column 1 t0 42 is clonal abuandance data frame (multilocusgenotype abundances)
 Masterdata1<-Masterdata[-41,]		                                       ## remove Ip-20 for some of the analysis 
 head(Masterdata)
 haile=Masterdata1
 dim(haile)
  #-------------------------------------------------------------------------------
 spe1<-subset(haile,select=cbind(1:42))  			                          ## data frame of all 42 clones
 spe1
 dim(spe1)
 
 #-------------------------------------------------------------------------------
 env<- subset(haile, select = cbind (43:51))	                            
 ## data frame of (9) Environmental variables( TP, TN, Chl-a, Cond, pH, Temp, Alt, Dep and area)
 env
 head(env)

 #-------------------------------------------------------------------------------
 spa<- subset(Masterdata, select = cbind (53:54))		                          
 ## Geographic location with UTM reading collected in the field
 spa
 head(spa)
 dim(spa)
 
 Dist<- subset(haile, select = cbind (distance1))	                      
 ## distance (in km) from the Ice sheet terminalin SE direction
 Dist
 
 Dist2<- subset(haile, select = cbind (distance3))	                      
 ## distance (in km) from tongue of Ice sheet terminal in NE direction
 Dist2

#------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 
# Species indicator values (Dufrene and Legendre)
# ***********************************************
# Divide the sites into 2 groups depending on the distance to the icesheet/glacier toungue (i.e age groupd -old vs young)
 Age<- factor(haile[, 59])
 Age
 library(indicspecies)
 citation("indicspecies")
 Indicator_MLGs = multipatt(spe1, Age,control = permControl(nperm=999))
 summary(Indicator_MLGs)  #Displaying the results
 summary(Indicator_MLGs, Indicator_MLGscomp=TRUE)
 summary(Indicator_MLGs, alpha=1)
 Indicator_MLGs$sign
 
 icesheet.D1 <- dist(data.frame(das=haile[,52], row.names=rownames(env))) 
 icesheet.kmeans1 <- kmeans( icesheet.D1, centers=2, nstart=100)
 icesheet.kmeans1$cluster
 
 # Indicator species of the sites
 (iva <- multipatt(spe1,  icesheet.kmeans1$cluster))
 iva
 summary(iva)
 # Table of the significant indicator species
 gr <- iva$maxcls[iva$pval <= 0.05]
 iv <- iva$indcls[iva$pval <= 0.05]
 pv <- iva$pval[iva$pval <= 0.05]
 fr <- apply(spe1 > 0, 2, sum)[iva$pval <= 0.05]
 fidg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
 fidg <- fidg[order(fidg$group, -fidg$indval),]
 fidg
 #-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------  
 icesheet.D2 <- dist(data.frame(das=haile[,63], row.names=rownames(env))) 
 icesheet.kmeans2 <- kmeans( icesheet.D2, centers=2, nstart=100)
 icesheet.kmeans2$cluster
# Indicator species for this typology of the sites
(iva <- indval(spe1,  icesheet.kmeans2$cluster))
 (iva2 <- multipatt(spe1,  icesheet.kmeans2$cluster))
 iva2
 summary(iva2)
# Table of the significant indicator species
gr <- iva$maxcls[iva$pval <= 0.05]
iv <- iva$indcls[iva$pval <= 0.05]
pv <- iva$pval[iva$pval <= 0.05]
fr <- apply(spe > 0, 2, sum)[iva$pval <= 0.05]
fidg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fidg <- fidg[order(fidg$group, -fidg$indval),]
fidg
# Export the result to a CSV file (to be opened in a spreadsheet)
write.csv(fidg, "IndVal-das.csv")