library(ape)
library(vegan)

### for the PCoA analysis 
## load plant matrix

veg_psv[veg_psv!=0]<-1
veg_pdc[veg_pdc!=0]<-1

D<-vegdist(veg_psv,method="bray")  ###note that bray method for incidence data is soerensen
pcoa_psv<-pcoa(D,correction="none")

E<-vegdist(veg_pdc,method="bray")  
pcoa_pdc<-pcoa(E,correction="none")

Axes_psv<-as.data.frame(pcoa_psv$vectors)
Axes_pdc<-as.data.frame(pcoa_pdc$vectors)


### load diverse abundance/incidence matrices (txt) of psv and pdc

all_psv
all_pdc
lide_psv
lide_pdc
specis_psv
specis_pdc
macros_psv
macros_pdc
micros_psv
micros_pdc

### load factors of psv and pdc

facs_psv_new
facs_pdc_new

facs_psv_new[,"Axis.1"]<-Axes_psv$Axis.1
facs_psv_new[,"Axis.2"]<-Axes_psv$Axis.2
facs_pdc_new[,"Axis.1"]<-Axes_pdc$Axis.1
facs_pdc_new[,"Axis.2"]<-Axes_pdc$Axis.2

###capscale for all groups (example for all_species)

cap_psv<-capscale(all_psv~Mean_Basalarea+forestdensity+
                          Mean_T+Mean_humi+
                          Axis.1+Axis.2+Year,data=facs_psv_new,dist="bray")

cap_pdc<-capscale(all_pdc~Mean_Basalarea+forestdensity+
                    Mean_T+Mean_humi+
                    Axis.1+Axis.2+year,data=facs_pdc_new,dist="bray")

### make plots

plot(cap_psv,type="n")
points(cap_psv,display="sites",col=facs_psv_new$color)
text(cap_psv,display="sites",col="grey",cex=1)
points(cap_psv,display="bp",col="black")
text(cap_psv,display="bp",col="grey",cex=1)

plot(cap_pdc,type="n")
points(cap_pdc,display="sites",col=facs_pdc_new$color)
text(cap_pdc,display="sites",col="grey",cex=1)
points(cap_pdc,display="bp",col="black")
text(cap_pdc,display="bp",col="grey",cex=1)


### for explained variation

anova(cap_psv, by="terms", permutations=999)

### for axes

anova(cap_psv, by="axis", permutations=999)

### for getting incidence data (but please load abundance data again for randomizations!)

all_psv[all_psv!=0]<-1
all_pdc[all_pdc!=0]<-1


### take random sample with fixed number of species

psv_1<-t(all_psv)
pdc_1<-t(all_pdc)

random_psv_50<-psv_1[sample(nrow(psv_1),size=50),]
random_psv_100<-psv_1[sample(nrow(psv_1),size=100),]
random_psv_200<-psv_1[sample(nrow(psv_1),size=200),]

ran_psv50<-t(random_psv_50)
ran_psv100<-t(random_psv_100)
ran_psv200<-t(random_psv_200)

random_pdc_50<-pdc_1[sample(nrow(pdc_1),size=50),]
random_pdc_100<-pdc_1[sample(nrow(pdc_1),size=100),]
random_pdc_200<-pdc_1[sample(nrow(pdc_1),size=200),]

ran_pdc50<-t(random_pdc_50)
ran_pdc100<-t(random_pdc_100)
ran_pdc200<-t(random_pdc_200)

### replicate with 1000 resamples for PsV
### for 50 species/sample

ran<-function(x)
{a1<-anova((capscale((t(psv_1[sample(nrow(psv_1),size=50),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+Year,data=facs_psv_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
  x<-a1$SumOfSqs
  return(x)}

results_rans<-sapply(seq_len(1000),ran)

rowMeans(results_rans)

### for 100 species/sample

ran<-function(x)
{a1<-anova((capscale((t(psv_1[sample(nrow(psv_1),size=100),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+Year,data=facs_psv_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
x<-a1$SumOfSqs
return(x)}

results_ran_100<-sapply(seq_len(1000),ran)

rowMeans(results_ran_100)

### for 200 species/sample

ran<-function(x)
{a1<-anova((capscale((t(psv_1[sample(nrow(psv_1),size=200),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+Year,data=facs_psv_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
x<-a1$SumOfSqs
return(x)}

results_ran_200<-sapply(seq_len(1000),ran)

rowMeans(results_ran_200)

### replicate with 1000 resamples for PdC
### for 50 species/sample

rans_pdc<-function(x)
{a2<-anova((capscale((t(pdc_1[sample(nrow(pdc_1),size=50),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+year,data=facs_pdc_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
x<-a2$SumOfSqs
return(x)}

results_ranspdc<-sapply(seq_len(1000),rans_pdc)

rowMeans(results_rans)

### for 100 species/sample

rans_pdc<-function(x)
{a2<-anova((capscale((t(pdc_1[sample(nrow(pdc_1),size=100),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+year,data=facs_pdc_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
x<-a2$SumOfSqs
return(x)}

results_ranpdc_100<-sapply(seq_len(1000),rans_pdc)

rowMeans(results_ranpdc_100)

### for 200 species/sample

rans_pdc200<-function(x)
{a2<-anova((capscale((t(pdc_1[sample(nrow(pdc_1),size=200),]))~Mean_Basalarea+forestdensity+Mean_T+Mean_humi+Axis.1+Axis.2+year,data=facs_pdc_new,dist="bray")),by="terms",model="direct",perm.max=9999,step=1000)
x<-a2$SumOfSqs
return(x)}

results_ranpdc_200<-sapply(seq_len(1000),rans_pdc200)

rowMeans(results_ranpdc_200)

#### build a nice graph for rē
###load bar_rsqu (table with all rē values)

# library
library(ggplot2)

# Stacked
ggplot(bar_rsqu, aes(fill=Factor, y=value, x=Site)) + 
  geom_bar(position="stack", stat="identity")

### or summed up to 100%

ggplot(bar_rsqu, aes(fill=Factor, y=value, x=Site)) + 
  geom_bar(position="fill", stat="identity")