--- 
title: "frontier_chemical_ecology_yoneya2022 v3 (final)"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Ctrl+Shift+Enter*. 

## Loading libraries
```{r}
library(vegan)
library(ggplot2)
library(ggsci)
library(gplots)
library(iNEXT)
library(nlme)
library(rcompanion)
library(car)
set.seed(4321)
```
## Loading datasets
```{r}
uds_data5.temp <- read.csv("monthly-for-R3.csv",header=T)  #for 2 or 3 continuous sampling at day 10, 32, and 60

uds_data_all_day <-uds_data5.temp[,c(-1)]
uds_data_all_day2 <-uds_data5.temp

#converting treatment to chamber
chamber <- uds_data_all_day2$treatment
chamber[chamber == "N"] <- "C1"
chamber[chamber == "W"] <- "C2"
chamber[chamber == "T"] <- "C2"
uds_data_all_day2 <- cbind.data.frame(uds_data_all_day2, chamber)

uds_data_all_LBa <- read.csv("leaf_beetle_adult_summary.csv", header=T) #for leaf beetle adult data
uds_dataLB <- uds_data_all_day[, c(3)]
head(uds_data_all_LBa)
head(uds_data_all_day2)
```
## Note for symbols to represent plant initial conditions
- "U" and "N" were interchangeably used to represent "Uninfested plants" in chamber C1
- "E", "S", and "T" were interchangeably used to represent "Exposed plants" in chamber C2
- "D" and "W" were interchangeably used to represent "Damaged plants" in chamber C2

## Basic Analysis of the time evolution (succession) of the abundant species (for supporting information)
### Data pretreatment
```{r}
SAD_total <- data.frame(abundance=apply(uds_data_all_day2[,4:76], 2, sum))
SAD_rank <- order(SAD_total[,1], decreasing=T)
SAD_total2 <- data.frame(rank=c(1:73), ID=rownames(SAD_total)[SAD_rank], abundance=SAD_total[SAD_rank,1])
#how many species?
sum(SAD_total2$abundance>0)

plot(SAD_total2$rank, log10(SAD_total2$abundance)) #SAD, including the rank 1 speices (wl) but with log10 scale
plot(SAD_total2$rank[-1], SAD_total2$abundance[-1]) #SAD, excluding the rank 1 speices (wl)
plot(SAD_total2$rank[2:10], SAD_total2$abundance[2:10]) #SAD, from rank 2 to rank 10
#text(SAD_total2$rank[2:10], SAD_total2$abundance[2:10], labels=SAD_total2$ID[2:10])
ggplot(SAD_total2, aes(x=rank, y=log10(abundance))) + geom_point()
SAD_total2_sub <- SAD_total2[2:10,]
sad_plot <- ggplot(SAD_total2_sub, aes(x=rank, y=abundance)) + geom_point()
sad_plot <- sad_plot + geom_text(aes(label=ID), size=3,vjust=-1)
plot(sad_plot)
```
### Time evolution of the rank 1-10 taxa
```{r}
#Taxon (species) identity
SAD_total2_sub$ID
uds_data_all_day_SI <- uds_data_all_day2  #copy for SI only
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="N"] <- "U"
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="W"] <- "D"
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="T"] <- "E"
treatment_day <- as.factor(paste(uds_data_all_day_SI$treatment,uds_data_all_day_SI$day))
uds_data_all_day_SI <- cbind.data.frame(uds_data_all_day_SI, treatment_day)
plot(uds_data_all_day_SI$treatment_day ,uds_data_all_day_SI$wl)

#rank 1
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=wl, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.2) + ggtitle("Plagiodera versicolora")

#rank 2
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=spider3, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("spider group 3 (< 5mm with spiderweb)")

#rank 3
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep9, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Caloptilia sp.2")

#rank 4
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep7, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Saliciphaga caesia")

#rank 5
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep10, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Phyllocolpa sp")

#rank 6
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=aphi5, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Cavariella salicicola")

#rank 7
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep12, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("spider group 4 (< 5mm without spiderweb)")

#rank 8
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=parasitoid3, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Tachinidae sp.1 (host: P. versicolora)")

#rank 9
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep12, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Microleon longipalpis")

#rank 10
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep13, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Monema flavescens")

#rank 11
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=grasshopper1, fill=treatment))  + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Orthoptera sp.1")


```



## Leaf beetle population dynamics analysis (Fig.1)
### Incorporating nested structure
```{r}
chamber <- append(rep("C1",16),rep("C2",32))
uds_data_all_LBa <- cbind.data.frame(uds_data_all_LBa, chamber)
uds_data_all_LBa
```
### Repeated ANOVA with autoregression
Ref: http://rcompanion.org/handbook/I_09.html
```{r}
#Data decomposition
uds_data_LBa24h <- subset(uds_data_all_LBa, uds_data_all_LBa$day==1)
uds_data_LBa2_7d <- subset(uds_data_all_LBa, (uds_data_all_LBa$day >= 2 & uds_data_all_LBa$day <= 7))
uds_data_LBa10_60d <- subset(uds_data_all_LBa, (uds_data_all_LBa$day >= 10 & uds_data_all_LBa$day <= 60))
```
#### hour scale (Fig.1a)
```{r}
#short (Fig.1a)
model24h_0 = gls(adult ~ chamber + hour + chamber*hour,data=uds_data_LBa24h)
#Anova(model24h_0)
ACF_s <- ACF(model24h_0, form = ~ hour|id)
head(ACF_s)
model24h = gls(adult ~ chamber + hour + chamber*hour, correlation = corAR1(form = ~ hour|id, value = ACF_s$ACF[2]),data=uds_data_LBa24h)
Anova(model24h)
#Post-hoc test
data12h<-subset(uds_data_LBa24h, hour==12)
pairwise.t.test(data12h$adult, data12h$chamber,p.adjust.method="holm")
pairwise.t.test(data12h$adult[-c(1:16)],data12h$treatment[-c(1:16)],p.adjust.method="holm")

data24h<-subset(uds_data_LBa24h, hour==24)
pairwise.t.test(data24h$adult, data24h$chamber,p.adjust.method="holm")
pairwise.t.test(data12h$adult[-c(1:16)], data12h$treatment[-c(1:16)],p.adjust.method="holm")

#Summary statistics for graph
sum24h = groupwiseMean(adult ~ treatment + hour, data = uds_data_LBa24h, conf = 0.95, digits = 3,traditional = FALSE, percentile  = TRUE)
sum24h

#shortest graph (Fig.1a)
pd = position_dodge(1)
g_line<-ggplot(sum24h, aes(y=Mean, x=hour)) + xlim(0,25) + ylim(0,1.1) + coord_fixed(25/1.1)
g_line<-g_line+geom_line(aes(colour=treatment),size=1,position=pd) 
g_line<- g_line+ geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd) 
print(g_line)
```
#### Summary of hour scale results (Fig.1a)
- There were statistically significant differences bet. U vs (S + D), i.e., chamber effect.
- By the add-hoc test,
  - there were differences bet. U vs (S + D) at 12 h and 24 h 
  - but no differences bet. S vs D either. 

#### a week scale (Fig.1b)
```{r}
#middle (Fig.1b)
model2_7d0 = gls(adult ~ chamber + day + chamber*day,data=uds_data_LBa2_7d)
#Anova(model2_7d0)
ACF_m <- ACF(model2_7d0, form = ~ day|id)
head(ACF_m)
model2_7d = gls(adult ~ chamber + day + chamber*day, correlation = corAR1(form = ~ day|id, value = ACF_m$ACF[2]),data=uds_data_LBa2_7d)
Anova(model2_7d)
#Post-hoc test
data7d<-subset(uds_data_LBa2_7d, day==7)
pairwise.t.test(data7d$adult, data7d$chamber,p.adjust.method="holm")
pairwise.t.test(data7d$adult[-c(1:16)], data7d$treatment[-c(1:16)],p.adjust.method="holm")

#Summary statistics for graph
sum2_7d = groupwiseMean(adult ~ treatment + day, data = uds_data_LBa2_7d, conf = 0.95, digits = 3,traditional = FALSE, percentile  = TRUE)

#a week scale graph (Fig.1b)
pd = position_dodge(0.3)
g_line<-ggplot(sum2_7d, aes(y=Mean, x=day)) + xlim(1,8) + ylim(0,1.1) + coord_fixed(7/1.1)
g_line<-g_line+geom_line(aes(colour=treatment),size=1, position=pd) 
g_line<- g_line+ geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd) 
print(g_line)
```

#### Summary of a week scale results (Fig.1b)
- There was no statistically significant difference bet. U vs (S + D), i.e., chamber effect.
- By the add-hoc test,
  - there was no difference bet. U vs (S + D) at 7d
  - but no differences bet. S vs D either. 

#### monthly scale (Fig.1c)
```{r}
model10_60d0 = gls(adult ~ chamber + day + chamber*day, data=uds_data_LBa10_60d)
#Anova(model10_60d0)
ACF_l <- ACF(model10_60d0, form = ~ day|id)
head(ACF_l)
model10_60d = gls(adult ~ chamber + day + chamber*day, correlation = corAR1(form = ~ day|id, value = ACF_l$ACF[2]),data=uds_data_LBa10_60d)
Anova(model10_60d)
#Post-hoc test
data10d<-subset(uds_data_LBa10_60d, day==10)
pairwise.t.test(data10d$adult, data10d$chamber,p.adjust.method="holm")
pairwise.t.test(data10d$adult[-c(1:16)], data10d$treatment[-c(1:16)],p.adjust.method="holm")

#Summary statistics for graph
sum10_60d = groupwiseMean(adult ~ treatment + day, data = uds_data_LBa10_60d, conf = 0.95, digits = 3,traditional = FALSE, percentile  = TRUE)

#graphics
pd = position_dodge(3)
g_line<-ggplot(sum10_60d, aes(y=Mean, x=day)) + xlim(5,65) + ylim(0,4) + coord_fixed(ratio = 60/4)
g_line<-g_line+geom_line(aes(colour=treatment),size=1, position=pd) 
g_line<- g_line+geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd) 
print(g_line)
```

#### Summary of monthly scale results (Fig.1c)
- There was no statistically significant difference bet. U vs (S + D), i.e., chamber effect.
- By the add-hoc test,
  - there was no difference bet. U vs (S + D) at 10d
  - but a difference bet. S vs D at 10d. 

## Ordinations for community composition
### Preparing the function for Fig.2 & Fig.S1
```{r}
disper_uds_centroid_graph <- function(data_name, c_LB, meth, adj=NULL, type="violin", lim_min=-0.6, lim_max=0.4) {
  if(c_LB == "wLB") {
    uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    uds_each<-uds_data1[,c(-1,-2,-3, -77,-78, -79, -80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
    uds_each<-uds_data1[,c(-1,-2,-3,-4,-77,-78,-79, -80)]
  }  #When leaf beele is NOT included
  if(c_LB == "LBonly") {
    uds_data1 <- subset(data_name, abundance2 > 0)
    uds_each <- subset(uds_data1, day==date)[, c(4)]  #picking up LB only
  }
  
  sample_each <- uds_data1$sample
  #meth="chao"
  uds_bray_each0.d<-vegdist(uds_each, method=meth)
  if(c_LB == "LBonly") tp <- "centroid"
  else tp <- "median"
  uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
  uds_median <- data.frame(uds_dist_each0$centroids)  #calculate the centroid (spatial median) for repeated measures
  
  sample <- row.names(uds_median)
  #sample <- sample_each
  #Replace the sample name into the treatment
  sample[grep("1n", sample)]<- "U10"
  sample[grep("2n", sample)]<- "U32"
  sample[grep("3n", sample)]<- "U60"
  sample[grep("1t", sample)]<- "S10"
  sample[grep("2t", sample)]<- "S32"
  sample[grep("3t", sample)]<- "S60"
  sample[grep("1w", sample)]<- "D10"
  sample[grep("2w", sample)]<- "D32"
  sample[grep("3w", sample)]<- "D60"
  
  treatment_each <- as.factor(sample) #just conversion

  uds_each.d<-vegdist(uds_median, method="euclidean")

  uds_dist_each <-betadisper(uds_each.d, treatment_each, type="centroid", bias.adjust = TRUE)  #Choose centroid because of euclidean distance
  #Plot
  if(c_LB == "LBonly") boxplot(uds_dist_each)
  else {
    #boxplot(uds_dist_each)
    #ordiplot(uds_dist_each$centroids, type="text", choices=c(1,2), cex=1, xlim=c(-0.6,0.4), ylim=c(-0.6, 0.5))
    #par(new=T)

    if(type == "violin") {
      #boxplot(uds_dist_each)
      ##from ggplot2
      #convert to data.frame (which is necessary for ggplot2)
      
      uds_dist_each2 <- cbind.data.frame(uds_dist_each$distances, as.character(uds_dist_each$group))
      
      treatment<-as.character(uds_dist_each$group)
      #Replace the sample name into the treatment
      treatment[grep("U", treatment)]<- "U"
      treatment[grep("D", treatment)]<- "D"
      treatment[grep("S", treatment)]<- "S"
      
      uds_dist_each2<-cbind.data.frame(uds_dist_each2, treatment)  
      colnames(uds_dist_each2) <- c("distances", "group", "treatment")
      uds_dist_each2$distances <- as.numeric(as.character(uds_dist_each2$distances))
      
      g_box<-ggplot(uds_dist_each2, aes(y=distances, x=group,colour=treatment))
      g_box<-g_box + geom_violin()
      g_box<-g_box + geom_boxplot(fill="black", width=0.1) 
      g_box<-g_box + coord_fixed(ratio = 15/1)
      print(g_box)
    }      
      else if(type == "PCoA") {
        #PCoA 
        test0<-capscale(uds_each.d~1)  
        test<-cmdscale(uds_each.d, eig=FALSE, k = 2)
        #prepare dataframe for each point
        test2<-summary(test)
        #uds_dist_each_test<-data.frame(PCoA1=-test2$sites[,1], PCoA2=-test2$sites[,2])
        uds_dist_each_test<-as.data.frame(test)
        uds_dist_each_test<-cbind.data.frame(sample, uds_dist_each_test)
        treatment<-as.character(sample)
        day<-as.character(sample)
        #Replace the sample name into the treatment & day
        treatment[grep("U", treatment)]<- "U"
        treatment[grep("D", treatment)]<- "D"
        treatment[grep("S", treatment)]<- "S"
        day[grep("10", day)]<-"Day10"
        day[grep("32", day)]<-"Day32"
        day[grep("60", day)]<-"Day60"
        
        uds_dist_each_test<-cbind.data.frame(cbind.data.frame(treatment, day), uds_dist_each_test)
        colnames(uds_dist_each_test)<-c("treatment", "day", "group", "PCoA1", "PCoA2")
        class(uds_dist_each_test$day)
        #prepare dataframe for centroid
        uds_dist_each_centr<-as.data.frame(uds_dist_each$centroids)
        treatment<-as.character(rownames(uds_dist_each_centr))
        day<-as.character(rownames(uds_dist_each_centr))
        #Replace the sample name into the treatment & day
        treatment[grep("U", treatment)]<- "U"
        treatment[grep("D", treatment)]<- "D"
        treatment[grep("S", treatment)]<- "S"
        day[grep("10", day)]<-"Day10"
        day[grep("32", day)]<-"Day32"
        day[grep("60", day)]<-"Day60"
        
        uds_dist_each_centr<-cbind.data.frame(uds_dist_each_centr, data.frame(treatment=treatment, day=day, group=rownames(uds_dist_each_centr)))
        
        g_PCoA <- ggplot(uds_dist_each_test, aes(x=PCoA1,y=PCoA2))
        g_PCoA <- g_PCoA + geom_point(aes(colour=treatment, shape=day), size=3, alpha=0.5)
        g_PCoA <- g_PCoA + xlim(lim_min, lim_max) + ylim(lim_min,lim_max) + coord_fixed(ratio = 1/1)
        #Add the centroid information
        g_PCoA <- g_PCoA + layer(
          data=uds_dist_each_centr, 
          mapping=aes(x=PCoA1, y=PCoA2, colour=treatment, shape=day), 
          params=list(size=5),
          geom="point", 
          stat="identity", 
          position="identity"
        )
        print(g_PCoA)
       
        #Performance of PCoA  http://d.hatena.ne.jp/hoxo_m/20120313/p1
        summary(test0)
    }
  }
  
  #Performance of PCoA  http://d.hatena.ne.jp/hoxo_m/20120313/p1
  #summary(test0)
  
  if(type == "NMDS") {
    uds_each.mds<-metaMDS(uds_median, distance="euclidean")
    rownames(uds_each.mds$points)<-treatment_each
    ordiplot(uds_each.mds$points, type='points', xlim=c(-1.0, 1.0))
    text(uds_each.mds$points,labels=treatment_each, col=unclass(treatment_each),cex=0.8)
  }

}
```
### Figure 2a & 2b
- With excluding leaf beetle and using chao distance
```{r}
#Figure 2a
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="chao", type="violin")
#Figure 2b
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="chao", type="PCoA")
```
### Figure S1
- (a) Using bray-curtis dissimilarity w/o LB
- (b) Using Chao dissimilarity w/ LB
```{r}
####For Figure S1a
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="bray", type="PCoA", lim_min=-0.75, lim_max=0.5) 
####For Figure S1b
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="wLB", meth="chao", type="PCoA", lim_min=-0.75, lim_max=0.5)   
```

### Statistical tests corresponding to the results shown in Fig.2
#### Definitions of functions fot statistical test
```{r}
#///////////////////////Pre-Evaluation of effectiveness of PERMANOVA (or can be used for assessing beta diversity) through PERMDISP for each of 10d, 32d, 60d/////////////////////
disper_uds_centroid <- function(data_name, date, c_LB, meth, adj=NULL, plot="box", dist_out=FALSE) {
    
    if(c_LB == "wLB") {
      uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
      uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3, -77,-78, -79, -80)]
    }  #When leaf beetle is included
    if(c_LB == "woLB") {
      uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
      uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-4,-77,-78,-79, -80)]
    }  #When leaf beetle is NOT included 
  
    if(c_LB == "LBonly") {  #For leaf beetle only but consider the case with other insects only
      uds_data1 <- subset(data_name, abundance2 > 0)
      #uds_data1 <-data_name
      uds_each <- subset(uds_data1, day==date)[, c(4)]  #picking up LB only
    }
    
    sample_each <- subset(uds_data1, day==date)$sample
    #length(uds_each)
    #length(sample_each)
    
    ##In case of leaf beetle only, different calculation
    if(c_LB == "LBonly") {
      sample_ID<-as.factor(sample_each)
      LB_average<-tapply(uds_each, sample_ID, mean)  #calculate the average of short-term repeated measurement
      LB_average2<-LB_average #just copy
      sample <- row.names(LB_average2)
      sample2 <- row.names(LB_average2)
      #Replace the sample name into the treatment
      sample[grep("n", sample)]<- "U" #undamaged
      sample[grep("t", sample)]<- "S" #Signaled
      sample[grep("w", sample)]<- "D" #Damaged
      treatment_each <- as.factor(sample) #just conversion
      #print(treatment_each)
      
      #Replace the sample name into the chamber
      sample2[grep("n", sample2)]<- "C1" #undamaged
      sample2[grep("t", sample2)]<- "C2" #Signaled
      sample2[grep("w", sample2)]<- "C2" #Damaged
      chamber_each <- as.factor(sample2) #just conversion
      
      average_treatment<-tapply(LB_average, treatment_each, mean)   #calculate the average of U, D, S, respectively
      average_chamber<-tapply(LB_average2, chamber_each, mean)   #calculate the average of C1 amd C2
      #print(average_treatment)
      for(i in 1: length(LB_average)) {  #normalizing the original data to exclude the effect of average on the distance to centroid
        if(treatment_each[i] == "D" && average_treatment[1] > 0.0) LB_average[i] = LB_average[i]/average_treatment[1] 
        if(treatment_each[i] == "S" && average_treatment[2] > 0.0) LB_average[i] = LB_average[i]/average_treatment[2]
        if(treatment_each[i] == "U" && average_treatment[3] > 0.0) LB_average[i] = LB_average[i]/average_treatment[3]
      }
      
      for(i in 1: length(LB_average)) {  #normalizing the original data to exclude the effect of average on the distance to centroid
        if(chamber_each[i] == "C1" && average_chamber[1] > 0.0) LB_average2[i] = LB_average2[i]/average_chamber[1] 
        if(chamber_each[i] == "C2" && average_chamber[2] > 0.0) LB_average2[i] = LB_average2[i]/average_chamber[2]
      }
      
    uds_median<-LB_average  #send the data for further analysis, distinguishing D,S,and U
    uds_median2<-LB_average2  #send the data for further analysis, distinguishing C1 and C2
    #print(uds_median)
    if(dist_out == "C1") uds_median <- subset(uds_median, chamber_each=="C1")
    else if(dist_out == "C2") uds_median <- subset(uds_median, chamber_each=="C2")
    #print(uds_median)
    }
    
    else {
      uds_bray_each0.d<-vegdist(uds_each, method=meth)
      
      tp <- "median"
      uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
      uds_median <- data.frame(uds_dist_each0$centroids)  #calculate the centroid (spatial median) for repeated measures
      
      sample <- row.names(uds_median)
      sample2 <- sample #just copy
      #Replace the sample name into the treatment
      sample[grep("n", sample)]<- "U" #undamaged
      sample[grep("t", sample)]<- "S" #Signaled
      sample[grep("w", sample)]<- "D" #Damaged
      treatment_each <- as.factor(sample) #just conversion
      #Replace the sample name into the chamber
      sample2[grep("n", sample2)]<- "C1" #undamaged
      sample2[grep("t", sample2)]<- "C2" #Signaled
      sample2[grep("w", sample2)]<- "C2" #Damaged
      chamber_each <- as.factor(sample2) #just conversion
      
      if(dist_out == "C1") uds_median <- uds_median[chamber_each=="C1",]
      else if(dist_out == "C2") uds_median <- uds_median[chamber_each=="C2",]
      #print(uds_median)
    }
    
    
    uds_each.d<-vegdist(uds_median, method="euclidean")
    if(dist_out!= FALSE) return(uds_each.d)
    
    #treatment_each
    uds_dist_each <-betadisper(uds_each.d, treatment_each, type="centroid", bias.adjust = TRUE)  #Choose centroid because of euclidean distance
    #chamber_each
    uds_dist_each_C <-betadisper(uds_each.d, chamber_each, type="centroid", bias.adjust = TRUE)  
    #Plot
    if(c_LB == "LBonly") boxplot(uds_dist_each)
    else {
      if(plot=="box") {
        boxplot(uds_dist_each_C)
        boxplot(uds_dist_each)
      }
      else{
        plot(uds_dist_each)
        ordilabel(scores(uds_dist_each, "centroids"))
        plot(uds_dist_each_C)
        ordilabel(scores(uds_dist_each_C, "centroids"))
        
      }
    }
    
    #PERDISP & permutation 
    perm_result_C <- permutest(uds_dist_each_C, pairwise=FALSE, permutations=9999)
    perm_result <- permutest(uds_dist_each, pairwise=TRUE, permutations=9999)
    
    cat("Sampling day was\n")
    print(subset(uds_data1, day==date)$day[1])
    cat("Comparison between chambers\n")
    print(perm_result_C)
    cat("Comparison between treatments\n")
    print(perm_result)
    
    #In fact, the adjustment is not necessary because permutest.betadisper does not repeat permutation for each pair
    if(!is.null(adj)) {
      pairwise_adj<-p.adjust(perm_result$pairwise$permuted, method = adj , n = length(perm_result$pairwise$permuted))
      cat("Pairwise comparision by permutation, adjusted by")
      print(as.character(adj))
      print(pairwise_adj[1])
      print(pairwise_adj[2])
      print(pairwise_adj[3])
    }

    
}

#///PERMANOVA for each initial conditions, based on the spatial median for repeated measures
com_div_uds_centroid <- function(data_name, date, c_LB, meth, adj)
{
  if(c_LB == "wLB") {
    uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-77,-78, -79, -80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
    uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79, -80)]
    uds_each_LB<-subset(uds_data1, day==date)[,c(4)]  #extract leaf beetle abundance
  }  #When leaf beele is NOTincluded
  
  sample_each <- subset(uds_data1, day==date)$sample

  uds_each.d<-vegdist(uds_each, method=meth)
  uds_dist_each <-betadisper(uds_each.d, sample_each, type="median", bias.adjust = TRUE)
  uds_median <- data.frame(uds_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  

  sample <- row.names(uds_median)
  sample2 <- row.names(uds_median)
  #Replace the sample name into the treatment
  sample[grep("n", sample)]<- "U"
  sample[grep("t", sample)]<- "S"
  sample[grep("w", sample)]<- "D"
  treatment_each <- as.factor(sample) #just conversion
  
  #Replace the sample name into the chamber
  sample2[grep("n", sample2)]<- "C1" #undamaged
  sample2[grep("t", sample2)]<- "C2" #Signaled
  sample2[grep("w", sample2)]<- "C2" #Damaged
  chamber_each <- as.factor(sample2) #just conversion
      
  
  adonis_result_C<-adonis(uds_median ~ chamber_each, method="euclidean", permutations=9999)
  adonis_result_T<-adonis(uds_median ~ treatment_each, method="euclidean", permutations=9999)

  print(adonis_result_C) #Chamber effect
  print(adonis_result_T) #Treatment effect
  
  #For pairwise comparison (only between D and S)
  cat("The pairwise PERMANOVA (adonis) analysis:\n")
  uds_dataDS0<-subset(data_name,  treatment != "N") #We still need the symbols N, W, and T because it is used in the raw data files
  uds_dataSU0<-subset(data_name,  treatment != "W")
  uds_dataUD0<-subset(data_name,  treatment != "T")
  if(c_LB == "wLB") {
    uds_dataDS1 <-subset(uds_dataDS0, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    DS_each<-subset(uds_dataDS1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
    uds_dataSU1 <-subset(uds_dataSU0, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    SU_each<-subset(uds_dataSU1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
    uds_dataUD1 <-subset(uds_dataUD0, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    UD_each<-subset(uds_dataUD1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_dataDS1 <-subset(uds_dataDS0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    DS_each<-subset(uds_dataDS1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79, -80)]
    uds_dataSU1 <-subset(uds_dataSU0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    SU_each<-subset(uds_dataSU1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79,-80)]
    uds_dataUD1 <-subset(uds_dataUD0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    UD_each<-subset(uds_dataUD1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79,-80)]
  }  #When leaf beetle is NO Tincluded
  
  sample_DS <- subset(uds_dataDS1, day==date)$sample
  sample_SU <- subset(uds_dataSU1, day==date)$sample
  sample_UD <- subset(uds_dataUD1, day==date)$sample
  
  DS_bray_each.d<-vegdist(DS_each, method=meth)
  SU_bray_each.d<-vegdist(SU_each, method=meth)
  UD_bray_each.d<-vegdist(UD_each, method=meth)
  
  DS_dist_each <-betadisper(DS_bray_each.d, sample_DS, type="median", bias.adjust = TRUE)
  SU_dist_each <-betadisper(SU_bray_each.d, sample_SU, type="median", bias.adjust = TRUE)
  UD_dist_each <-betadisper(UD_bray_each.d, sample_UD, type="median", bias.adjust = TRUE)
  
  DS_median <- data.frame(DS_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  SU_median <- data.frame(SU_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  UD_median <- data.frame(UD_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  
  DS_sample <- row.names(DS_median)
  SU_sample <- row.names(SU_median)
  UD_sample <- row.names(UD_median)
  
  #Replace the sample name into the treatment
  DS_sample[grep("n", DS_sample)]<- "U"
  DS_sample[grep("t", DS_sample)]<- "S"
  DS_sample[grep("w", DS_sample)]<- "D"
  treatment_DS <- as.factor(DS_sample) #just conversion
  SU_sample[grep("n", SU_sample)]<- "U"
  SU_sample[grep("t", SU_sample)]<- "S"
  SU_sample[grep("w", SU_sample)]<- "D"
  treatment_SU <- as.factor(SU_sample) #just conversion
  UD_sample[grep("n", UD_sample)]<- "U"
  UD_sample[grep("t", UD_sample)]<- "S"
  UD_sample[grep("w", UD_sample)]<- "D"
  treatment_UD <- as.factor(UD_sample) #just conversion

  adonis_resultDS<-adonis(DS_median ~ treatment_DS, method="euclidean", permutations=9999)
  adonis_resultSU<-adonis(SU_median ~ treatment_SU, method="euclidean", permutations=9999)
  adonis_resultUD<-adonis(UD_median ~ treatment_UD, method="euclidean", permutations=9999)
  
  p_list <- c(adonis_resultDS$aov$Pr[1], adonis_resultSU$aov$Pr[1], adonis_resultUD$aov$Pr[1])
  
  cat("The comparision between D and S is:\n")
  print(adonis_resultDS)
  
  pairwise_adj<-p.adjust(p_list, method = adj , n = length(p_list))
  
  #cat("Pairwise comparision by permutation, adjusted by")
  #print(adj)
  #cat("The adjusted  P-value for S vs D was\n")
  #print(pairwise_adj[1])
  #cat("The adjusted  P-value for U vs S was\n")
  #print(pairwise_adj[2])
  #cat("The adjusted  P-value for D vs U was\n")
  #print(pairwise_adj[3])
  
}

#///PERMANOVA for each date, based on the spatial median for repeated measures
com_div_DAY_centroid <- function(data_name, treat, c_LB, meth, adj=NULL) {
  
  #data_name <-uds_data_all_day2
  #treat <- "T"
  #meth="chao"
  #date <- "10d","32d", or "60d", or any
  if(c_LB == "wLB") {
    uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3, -77,-78, -79,-80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
    uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
  }  #When leaf beele is NOTincluded
  if(c_LB == "LBonly") {
    uds_data1 <- subset(data_name, abundance2 > 0)
    uds_each <- subset(uds_data1, treatment==treat)[, c(4)]  #picking up LB only
  }
  
  sample_each <- subset(uds_data1, treatment==treat)$sample
  day_each <- subset(uds_data1, treatment==treat)$day

  uds_bray_each0.d<-vegdist(uds_each, method=meth)
  if(c_LB == "LBonly") tp <- "centroid"
  else tp <- "median"
  uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
  uds_median <- data.frame(uds_dist_each0$centroids)  #calculate the centroid (spatial median) for repeasted measures

  sample <- row.names(uds_median)
  
  #Replace the sample name into the date
  sample[grep("^1", sample)]<- "day10"
  sample[grep("^2", sample)]<- "day32"
  sample[grep("^3", sample)]<- "day60"
  
  day_each <- as.factor(sample) #just conversion
  
  adonis_result<-adonis(uds_median ~ day_each, method="euclidean", permutations=9999)
  
  print(adonis_result)
  
  #For pairwise comarison (with Bonferroni adjustment)
  cat("The pairwise PERMANOVA (adonis) analysis:\n")
  uds_dataDay1032_0<-subset(data_name,  day != "60d")
  uds_dataDay1060_0<-subset(data_name,  day != "32d")
  uds_dataDay3260_0<-subset(data_name,  day != "10d")
  
  if(c_LB == "wLB") {
    uds_dataDay1032_1 <-subset(uds_dataDay1032_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day1032_each<-subset(uds_dataDay1032_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
    uds_dataDay1060_1 <-subset(uds_dataDay1060_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day1060_each<-subset(uds_dataDay1060_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
    uds_dataDay3260_1 <-subset(uds_dataDay3260_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day3260_each<-subset(uds_dataDay3260_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_dataDay1032_1 <-subset(uds_dataDay1032_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day1032_each<-subset(uds_dataDay1032_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
    uds_dataDay1060_1 <-subset(uds_dataDay1060_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day1060_each<-subset(uds_dataDay1060_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
    uds_dataDay3260_1 <-subset(uds_dataDay3260_0, abundance2 > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    Day3260_each<-subset(uds_dataDay3260_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
  }  #When leaf beele is NOTincluded
  
  sample_1032 <- subset(uds_dataDay1032_1, treatment==treat)$sample
  sample_1060 <- subset(uds_dataDay1060_1, treatment==treat)$sample
  sample_3260 <- subset(uds_dataDay3260_1, treatment==treat)$sample

  Day1032_each.d<-vegdist(Day1032_each, method=meth)
  Day1060_each.d<-vegdist(Day1060_each, method=meth)
  Day3260_each.d<-vegdist(Day3260_each, method=meth)
  
  Day1032_dist_each <-betadisper(Day1032_each.d, sample_1032, type="median", bias.adjust = TRUE)
  Day1060_dist_each <-betadisper(Day1060_each.d, sample_1060, type="median", bias.adjust = TRUE)
  Day3260_dist_each <-betadisper(Day3260_each.d, sample_3260, type="median", bias.adjust = TRUE)

  Day1032_median <- data.frame(Day1032_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  Day1060_median <- data.frame(Day1060_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  Day3260_median <- data.frame(Day3260_dist_each$centroids)  #calculate the centroid (spatial median) for repeasted measures
  
  Day1032_sample <- row.names(Day1032_median)
  Day1060_sample <- row.names(Day1060_median)
  Day3260_sample <- row.names(Day3260_median)
  
  #Replace the sample name into the day
  Day1032_sample[grep("^1", Day1032_sample)]<- "day10"
  Day1032_sample[grep("^2", Day1032_sample)]<- "day32"
  day_Day1032 <- as.factor(Day1032_sample) #just conversion
  Day1060_sample[grep("^1", Day1060_sample)]<- "day10"
  Day1060_sample[grep("^3", Day1060_sample)]<- "day60"
  day_Day1060 <- as.factor(Day1060_sample) #just conversion
  Day3260_sample[grep("^2", Day3260_sample)]<- "day32"
  Day3260_sample[grep("^3", Day3260_sample)]<- "day60"
  day_Day3260 <- as.factor(Day3260_sample) #just conversion

  adonis_resultDay1032<-adonis(Day1032_median ~ day_Day1032, method="euclidian", permutations=9999)
  adonis_resultDay1060<-adonis(Day1060_median ~ day_Day1060, method="euclidian", permutations=9999)
  adonis_resultDay3260<-adonis(Day3260_median ~ day_Day3260, method="euclidian", permutations=9999)
  
  p_list <- c(adonis_resultDay1032$aov$Pr[1], adonis_resultDay1060$aov$Pr[1], adonis_resultDay3260$aov$Pr[1])
  
  
  pairwise_adj<-p.adjust(p_list, method = adj , n = length(p_list))
  
  cat("Pairwise comparision by permutation, adjusted by")
  print(adj)
  cat("The adjusted  P-value for day10 vs day32 was\n")
  print(pairwise_adj[1])
  cat("The adjusted  P-value for day10 vs day60 was\n")
  print(pairwise_adj[2])
  cat("The adjusted  P-value for day32 vs day60 was\n")
  print(pairwise_adj[3])
}


#///PERMDISP or each date, based on the spatial median for repeated measures
disper_DAY_centroid <- function(data_name, treat, c_LB, meth, adj=NULL) {
  
  #resolving inconsistency of the symbols to represent the treatments (and chamber)
  if(treat == "U") treat <- "N"
  else if(treat == "E") treat <- "T"
  else if(treat == "D") treat <- "W"
  else if(treat == "S") treat <- "T"
  
  if(c_LB == "wLB") {
    uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3, -77,-78, -79,-80)]
  }  #When leaf beele is included
  if(c_LB == "woLB") {
    uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
    uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
  }  #When leaf beele is NOTincluded
  if(c_LB == "LBonly") {
    uds_data1 <- subset(data_name, abundance2 > 0)
    uds_each <- subset(uds_data1, treatment==treat)[, c(4)]  #picking up LB only
  }
  
  sample_each <- subset(uds_data1, treatment==treat)$sample
  day_each <- subset(uds_data1, treatment==treat)$day
  
  
  uds_bray_each0.d<-vegdist(uds_each, method=meth)
  if(c_LB == "LBonly") tp <- "centroid"
  else tp <- "median"
  uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
  uds_median <- data.frame(uds_dist_each0$centroids)  #calculate the centroid (spatial median) for repeasted measures
  
  sample <- row.names(uds_median)
  
  #Replace the sample name into the date
  sample[grep("^1", sample)]<- "day10"
  sample[grep("^2", sample)]<- "day32"
  sample[grep("^3", sample)]<- "day60"
  
  day_each <- as.factor(sample) #just conversion
  
  #treatment_each
  uds_each.d<-vegdist(uds_median, method="euclidean")
  uds_dist_each <-betadisper(uds_each.d, day_each, type="centroid", bias.adjust = TRUE)  #Choose centroid becasue of euclidean distance
  #Plot
  if(c_LB == "LBonly") boxplot(uds_dist_each)
  else {
    boxplot(uds_dist_each)
    ordilabel(scores(uds_dist_each, "centroids"))
  }
  #PERDISP & ANOVA & MULTIPLE COMPARISON
  #anova_dist <- anova(uds_dist_each)
  #TukeyHSD(uds_dist_each, which = "group", ordered = FALSE, conf.levels = 0.95)
  
  #PERDIPST & permutation 
  perm_result <- permutest(uds_dist_each, pairwise=TRUE, permutations=9999)
  #perm_result
  
  #print (anova_dist$Pr[1]) 
  cat("Target treatment was\n")
  print(treat)
 
  print(uds_dist_each)
  print(perm_result)
  #In fact, the adjustment is not necessary because permutest.betadisper does not repeat permutation for each pair
  if(!is.null(adj)) {
    pairwise_adj<-p.adjust(perm_result$pairwise$permuted, method = adj, n = length(perm_result$pairwise$permuted))
    cat("Pairwise comparision by permutation, adjusted by")
    print(as.character(adj))
    print(pairwise_adj[1])
    print(pairwise_adj[2])
    print(pairwise_adj[3])
  }
  
  #TukeyHSD(uds_dist_each)
  
}


```

#### Statistical tests for Fig.2
##### For PERMDISP: comparison between dates (Fig.2A)\
```{r}
#Temporal variability in beta diversity for manuscript
disper_DAY_centroid(data_name=uds_data_all_day2, treat="U", c_LB="woLB", meth="chao")
disper_DAY_centroid(data_name=uds_data_all_day2, treat="S", c_LB="woLB", meth="chao")
disper_DAY_centroid(data_name=uds_data_all_day2, treat="D", c_LB="woLB", meth="chao")

```
#####For PERMDISP: comparison between chambers and treatments (Fig.2a)\
```{r}
#Results from 2-3 times measurement in 10d, 32d, and 60d without leaf beetle through spatial median
#With the recommended measure (Chao index) (for Fig.2b)
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao")
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao")
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao")
```

##### For PERMANOVA (Fig.2b)\
```{r}
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao", adj="holm")
```

##### PERMANOVA (for Supplementary Material)
```{r}
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="N", c_LB="woLB", meth="chao", adj="holm")
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="T", c_LB="woLB", meth="chao", adj="holm")
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="W", c_LB="woLB", meth="chao", adj="holm")
```


#### Summary for the statistical tests on Figure 2
- Dispersion with time
  - Dispersion of U (Chamber 1) decreased with time (day10=day32, but day10 > day60)
  - No differences in dispersion bet. dates in S(E) and D treatments
- Dispersion between chambers and treatment\
  - Dispersion of U (Chamber 1) was smaller than chamber 2 (D+S) at 60d (P=0.0075).
  - No differences in dispersion bet. S(E) and D treatments at 60d
- PERMANOVA (with PERMDISP) results tell that:
  - the centroid community composition was not different between U vs (D + S) (no chamber effects).
  - the difference between D and E at 32d (F = 2.4132, df=2, P = 0.0036)
- Temporal dynamics of centroid within each treatment
  - For all U, D, and E, all pairwise PERMANOVA resulted in P < 0.001
  - But due to heterogeneous dispersion day10 day32 > day60 for U, we cannot conclude the centroid of U at day 60 is difference from day 10 not day 32.
  
#### Statistical tests for Figure S1
##### For PERMDISP: comparision between chambers and treatments with bray-curtis measure
```{r}
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="bray")
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="bray")
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="bray")
```
##### For PERMANOVA: comparision between chambers and treatments with bray-curtis measure
```{r}
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="bray", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="bray", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="bray", adj="holm")
```

##### For PERMDISP: comparision between chambers and treatments with Chao measure but with leaf beetle
```{r}
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="wLB", meth="chao")
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="wLB", meth="chao")
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="wLB", meth="chao")
```
##### For PERMANOVA: comparision between chambers and treatments with bray-curtis measure
```{r}
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="wLB", meth="chao", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="wLB", meth="chao", adj="holm")
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="wLB", meth="chao", adj="holm")
```

### Temporal divergence on each willow tree (Decided not to present)
#### Preparation of functions to calculate/visualize temporal divergence
```{r}
turnover_comparison_short <- function(data_name, date, c_LB, meth, output = NULL)
{
 
 #data_name <- uds_data_all_day2
 #date <- "10d"
 #meth<-"chao"
 #date <- "10d","32d", or "60d", or any
 if(c_LB == "wLB") {
   uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
   uds_data2 <-subset(uds_data1, day==date)
   uds_each<-uds_data2[,c(-1,-2,-3, -77,-78, -79,-80)]
 }  #When leaf beetle is included
 if(c_LB == "woLB") {
   uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
   uds_data2 <-subset(uds_data1, day==date)
   uds_each<-uds_data2[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
 }  #When leaf beetle is NOT included
 
 sample_each <- subset(uds_data1, day==date)$sample
 treatment_each <- subset(uds_data1, day==date)$treatment
 chamber_each <- subset(uds_data1, day==date)$chamber
 #View(chamber_each)
 time_each <- subset(uds_data1, day==date)$time
 
 #prepare dummy vectors
 sample_n <- NULL
 distance_n <-NULL
 treatment_n <-NULL
 chamber_n <-NULL
 timei_n <-NULL
 timej_n <-NULL
 
 #calculate distance between 1P, 2P, and 3P
 for(i in 1:(length(sample_each) - 1)) {
   for(j in (i + 1): length(sample_each)) {
     if(sample_each[i] == sample_each[j]) {
       sample_n <- append(sample_n, as.character(sample_each[i]))
       distance_n <- append(distance_n, as.numeric(vegdist(rbind(uds_each[i,], uds_each[j,])), method=meth))
       treatment_n <-append(treatment_n, as.character(treatment_each[i]))
       chamber_n <- append(chamber_n, as.character(chamber_each[i]))
       timei_n <-append(timei_n, as.character(time_each[i]))
       timej_n <-append(timej_n, as.character(time_each[j]))
     }
    }
 }
 turnover<- data.frame(sample=sample_n, distance=distance_n, treatment=treatment_n, chamber=chamber_n, timei=timei_n, timej=timej_n)

 turnover2 <-subset(turnover, !(timei == "1P" & timej=="3P"))  #exclude the combination of 1P & 3P (redundant)
 
 #turnoverN <- subset(turnover2, treatment=="N")  #separate data by treatment
 turnoverT <- subset(turnover2, treatment=="T")  #separate data by treatment
 turnoverW <- subset(turnover2, treatment=="W")  #separate data by treatment
 turnoverC1 <- subset(turnover2, chamber=="C1") #separate data by chamber
 turnoverC2 <- subset(turnover2, chamber=="C2") #separate data by chamber
 
 #calculate average of 1P vs 2P & 2P vs 3P if any
 #average_turnoverN<-by(as.numeric(turnoverN$distance), turnoverN$sample, mean)
 average_turnoverT<-by(as.numeric(turnoverT$distance), turnoverT$sample, mean)
 average_turnoverW<-by(as.numeric(turnoverW$distance), turnoverW$sample, mean)
 average_turnoverC1<-by(as.numeric(turnoverC1$distance), turnoverC1$sample, mean)
 average_turnoverC2<-by(as.numeric(turnoverC2$distance), turnoverC2$sample, mean)

 #average_N<-subset(average_turnoverN, average_turnoverN>=0)
 average_T<-subset(average_turnoverT, average_turnoverT>=0)
 average_W<-subset(average_turnoverW, average_turnoverW>=0)
 average_C1<-subset(average_turnoverC1, average_turnoverC1>=0)
 average_C2<-subset(average_turnoverC2, average_turnoverC2>=0)
 
 length(average_T)
 #trN<-NULL
 trT<-NULL
 trW<-NULL
 trC1<-NULL
 trC2<-NULL
 #for(i in 1:length(average_N)) trN <-append(trN, "N")
 for(i in 1:length(average_T)) trT <-append(trT, "E")
 for(i in 1:length(average_W)) trW <-append(trW, "D")
 for(i in 1:length(average_C1)) trC1 <-append(trC1, "C1")
 for(i in 1:length(average_C2)) trC2 <-append(trC2, "C2")
 
 #av_N<-data.frame(treatment=trN, distance=average_N)
 av_T<-data.frame(treatment=trT, distance=average_T)
 av_W<-data.frame(treatment=trW, distance=average_W)
 av_C1<-data.frame(chamber=trC1, distance=average_C1)
 av_C2<-data.frame(chamber=trC2, distance=average_C2)
 
 #summary dataframe
 turnover_short <- rbind.data.frame(av_T, av_W)
 turnover_shortC <- rbind.data.frame(av_C1, av_C2)
 
 #print(turnover_shortC$distance)
 #View(turnover_shortC)
 #View(turnover_short)
 #print(distance_n)
 
 g_box<-ggplot(turnover_shortC, aes(y=distance, x=chamber))
 g_box<-g_box + geom_violin()
 g_box<-g_box + geom_boxplot(width=0.1) 
 g_box<-g_box + coord_fixed(ratio = 5/1) + coord_cartesian(ylim = c(0, 1))
 print(g_box)
 g_box2<-ggplot(turnover_short, aes(y=distance, x=treatment))
 g_box2<-g_box2 + geom_violin()
 g_box2<-g_box2 + geom_boxplot(width=0.1) 
 g_box2<-g_box2 + coord_fixed(ratio = 5/1) + coord_cartesian(ylim = c(0, 1))
 print(g_box2)
 
 #plotmeans(turnover_shortC$distance~turnover_shortC$chamber)
 #plotmeans(turnover_short$distance~turnover_short$treatment)

 print(bartlett.test(turnover_shortC$distance~turnover_shortC$chamber))
 #One-way ANOVA with unequal-variance
 print(oneway.test(turnover_shortC$distance~turnover_shortC$chamber,var.equal=FALSE)) 
 
 print(bartlett.test(turnover_short$distance~turnover_short$treatment))
 #One-way ANOVA with unequal-variance
 print(oneway.test(turnover_short$distance~turnover_short$treatment,var.equal=FALSE))
 #Pairwise comparison with Welch's t test
 #print(pairwise.t.test(turnover_short$distance, turnover_short$treatment, p.adj="holm"))
 
 #if(output==TRUE) return(turnover_short)
 
}  

turnover_comparison_long <- function(data_name, from_date, c_LB, meth="chao", output = NULL)
{
  
  if(c_LB == "wLB") {
    uds_data1 <-subset(data_name, abundance > 0)  #delete the sampling data with 0 total abundance when leaf beetle
    uds_each<-uds_data1[,c(-1,-2,-3, -77,-78,-79,-80)]
  }  #When leaf beetle is included
  if(c_LB == "woLB") {
    uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
    uds_each<-uds_data1[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
  }  #When leaf beetle is NOT included
  
  sample_each <- uds_data1$sample
  chamber_each <- uds_data1$chamber
  treatment_each <- uds_data1$treatment
  time_each <- uds_data1$time
  
  uds_time.d<-vegdist(uds_each, method=meth)
  uds_dist_time <-betadisper(uds_time.d, sample_each, type="median", bias.adjust = TRUE)
  
  uds_day.median <- data.frame(uds_dist_time$centroids)  #calculate the centroid (spatial median) for repeated measures
 
  #Replace the sample name into the treatment & day 
  sample0 <- row.names(uds_day.median)
  sample1 <- sample0
  sample2 <- sample0
  sample3 <- sample0
 
  sample0<-gsub("1n", "n", sample0)
  sample0<-gsub("2n", "n", sample0)
  sample0<-gsub("3n", "n", sample0)
  sample0<-gsub("1t", "t", sample0)
  sample0<-gsub("2t", "t", sample0)
  sample0<-gsub("3t", "t", sample0)
  sample0<-gsub("1w", "w", sample0)
  sample0<-gsub("2w", "w", sample0)
  sample0<-gsub("3w", "w", sample0)
  
  sample1[grep("n", sample1)]<- "C1"
  sample1[grep("t", sample1)]<- "C2"
  sample1[grep("w", sample1)]<- "C2"
  
  sample2[grep("n", sample2)]<- "U"
  sample2[grep("t", sample2)]<- "E"
  sample2[grep("w", sample2)]<- "D"
  
  sample3[grep("1n", sample3)] <-"10d"
  sample3[grep("2n", sample3)] <-"32d"
  sample3[grep("3n", sample3)] <-"60d"
  sample3[grep("1t", sample3)] <-"10d"
  sample3[grep("2t", sample3)] <-"32d"
  sample3[grep("3t", sample3)] <-"60d"
  sample3[grep("1w", sample3)] <-"10d"
  sample3[grep("2w", sample3)] <-"32d"
  sample3[grep("3w", sample3)] <-"60d"
  
  sample_each <- as.factor(sample0)
  chamber_each <- as.factor(sample1)
  treatment_each <- as.factor(sample2) #just conversion
  day_each <- as.factor(sample3) #just conversion
  
  uds_day.median2<-cbind(chamber=chamber_each, uds_day.median)
  uds_day.median2<-cbind(treatment=treatment_each, uds_day.median2)
  uds_day.median2<-cbind(day=day_each, uds_day.median2)
  uds_day.median2<-cbind(sample_each, uds_day.median2)
  
  
  uds_day.median2[1,c(-1,-2,-3)] 
  uds_day.median2$sample_each[2]
  #View(uds_day.median2)
  
  #prepare dummy vectors
  sample_n <- NULL
  distance_n <-NULL
  chamber_n <-NULL
  treatment_n <-NULL
  dayi_n <-NULL
  dayj_n <-NULL
  
  #calculate distance between 10d, 32d, and 60d
  for(i in 1:(length(sample_each) - 1)) {
    for(j in (i + 1): length(sample_each)) {
      if(sample_each[i] == sample_each[j]) {
        sample_n <- append(sample_n, as.character(sample_each[i]))
        distance_n <- append(distance_n, as.numeric(vegdist(rbind.data.frame(uds_day.median2[i,c(-1,-2,-3,-4)], uds_day.median2[j,c(-1,-2,-3,-4)]), method="euclidean")))
        chamber_n <-append(chamber_n, as.character(chamber_each[i]))
        treatment_n <-append(treatment_n, as.character(treatment_each[i]))
        dayi_n <-append(dayi_n, as.character(day_each[i]))
        dayj_n <-append(dayj_n, as.character(day_each[j]))
      }
    }
  }
  turnover<- data.frame(sample=sample_n, distance=distance_n, chamber=chamber_n, treatment=treatment_n, dayi=dayi_n, dayj=dayj_n)
  
  if(from_date == "10d") turnover_long <-subset(turnover, (dayi == "10d" & dayj=="32d"))  
  #turnover_10d60d <-subset(turnover, (dayi == "10d" & dayj=="60d"))  #not used excluding the combination of 10d & 60d (redundanct)
  if(from_date == "32d") turnover_long <-subset(turnover, (dayi == "32d" & dayj=="60d"))  
 
  turnover_long_tr <- subset(turnover_long, chamber == "C2")
  
  View(turnover_long)
  View(turnover_long_tr)
  
  g_box<-ggplot(turnover_long, aes(y=distance, x=chamber))
  g_box<-g_box + geom_violin()
  g_box<-g_box + geom_boxplot(width=0.1) 
  g_box<-g_box + coord_fixed(ratio = 5/1)
  print(g_box)
  g_box2<-ggplot(turnover_long_tr, aes(y=distance, x=treatment))
  g_box2<-g_box2 + geom_violin()
  g_box2<-g_box2 + geom_boxplot(width=0.1) 
  g_box2<-g_box2 + coord_fixed(ratio = 5/1)
  print(g_box2)
  #plotmeans(turnover_long$distance~turnover_long$chamber)
  #plotmeans(turnover_long_tr$distance~turnover_long_tr$treatment)
  
  print(bartlett.test(turnover_long$distance~turnover_long$chamber))
  print(bartlett.test(turnover_long_tr$distance~turnover_long_tr$treatment))
  #One-way ANOVA with unequal-variance
  print(oneway.test(turnover_long$distance~turnover_long$chamber,var.equal=FALSE))
  print(oneway.test(turnover_long_tr$distance~turnover_long_tr$treatment,var.equal=FALSE))
  
}
```
####Calculate/Visualize within day turnover (which is not used in the manuscript)
```{r}
#(No difference at all @ 10 day between time)  ####For Fig. S1
#turnover_comparison_short(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao")  #error due to lack of replication within day
turnover_comparison_short(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao")
turnover_comparison_short(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao")
```

####Calculate/Visualize between dates turnover (not used)
```{r}

#(No difference at all @ 10 day between time) ####For Fig. S1
turnover_comparison_long(data_name=uds_data_all_day2, from_date="10d", c_LB="woLB", meth="chao")
turnover_comparison_long(data_name=uds_data_all_day2, from_date="32d", c_LB="woLB", meth="chao")
```
#### Summary of the species turnover analysis (Decided not to use)
- No differences between Chambers
- No differences between E and D

### Relationship between leaf beetle population size and community composition
#### Preparation for Mantel test
```{r}
LB_diff_C1_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C12")
LB_diff_C1_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C12")
LB_diff_C1_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C12")

comm_diff_C1_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C1")
comm_diff_C2_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C2")
comm_diff_C12_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C12")
comm_diff_C1_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C1")
comm_diff_C2_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C2")
comm_diff_C12_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C12")
comm_diff_C1_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C1")
comm_diff_C2_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C2")
comm_diff_C12_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C12")

```
#### Mantel test and simple plots
```{r}
plot(LB_diff_C1_d10.d, comm_diff_C1_d10.d)
mantel(LB_diff_C1_d10.d, comm_diff_C1_d10.d, method="pearson", permutations=9999)
plot(LB_diff_C2_d10.d, comm_diff_C2_d10.d)
mantel(LB_diff_C2_d10.d, comm_diff_C2_d10.d, method="pearson", permutations=9999)
plot(LB_diff_C12_d10.d, comm_diff_C12_d10.d)
mantel(LB_diff_C12_d10.d, comm_diff_C12_d10.d, method="pearson", permutations=9999)

plot(LB_diff_C1_d32.d, comm_diff_C1_d32.d)
mantel(LB_diff_C1_d32.d, comm_diff_C1_d32.d, method="pearson", permutations=9999)
plot(LB_diff_C2_d32.d, comm_diff_C2_d32.d)
mantel(LB_diff_C2_d32.d, comm_diff_C2_d32.d, method="pearson", permutations=9999)
plot(LB_diff_C12_d32.d, comm_diff_C12_d32.d)
mantel(LB_diff_C12_d32.d, comm_diff_C12_d32.d, method="pearson", permutations=9999)

plot(LB_diff_C1_d60.d, comm_diff_C1_d60.d)
#mantel(LB_diff_C1_d60.d, comm_diff_C1_d60.d, method="pearson", permutations=9999)
plot(LB_diff_C2_d60.d, comm_diff_C2_d60.d)
mantel(LB_diff_C2_d60.d, comm_diff_C2_d60.d, method="pearson", permutations=9999)
plot(LB_diff_C12_d60.d, comm_diff_C12_d60.d)
mantel(LB_diff_C12_d60.d, comm_diff_C12_d60.d, method="pearson", permutations=9999)

```
#### ggplot for mantel test
Only for day 32
```{r}
library(reshape)
#C1 only
m_LB <- as.matrix(LB_diff_C1_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C1_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)

g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C1 only")
g_points<-g_points+geom_point(size=2,alpha=0.3) 
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
summary(result)
result$coefficients[1]
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=2,colour="blue") # slope
print(g_points)

#C2 only
m_LB <- as.matrix(LB_diff_C2_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C2_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)

g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C2 only")
g_points<-g_points+geom_point(size=2,alpha=0.3) 
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
summary(result)
result$coefficients[1]
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=1,colour="blue") # slope
print(g_points)

#C1 + C2 
m_LB <- as.matrix(LB_diff_C12_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C12_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)

g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C1 + C2")
g_points<-g_points+geom_point(size=2,alpha=0.3) 
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
result$coefficients[1]
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=1,colour="blue") # slope
print(g_points)

```


#### Summary of mantel test
- Positive relationship between LB population size and community composition at day 32 for the whole data
- Positive relationship between LB population size and community composition at day 32 for C2 chamber
- No relationship between LB population size and community composition at day 32 for C1 chamber


## Species richness analysis (Fig.3 & Fig.S3)
### Preparation of functions with iNEXT package
```{r}

alpha_hill_number <- function(data_name, date, q)
{
  #data_name <-uds_data_all_day2
  #date <- "32d"
  #Including leaf_beetle
  
  uds_data1 <-subset(data_name, abundance2 >= 0) #OR, delete the sampling when species other than leaf beetle does not exist
  uds_data2 <-subset(uds_data1, day==date)
  
  #if excluding leaf beetle
  #uds_each<-uds_data2[,c(-1,-2,-3,-4,-77,-78,-79)]
  #OR inclcding leaf beetle
  uds_each<-uds_data2[,c(-1,-2,-3,-77,-78,-79,-80)]
  
  
  #Hill number q = 0, 1, 2
  uds_hill<-renyi(uds_each, scales=c(0,1,2), hill=TRUE)
  uds_hill2<-cbind(uds_hill, uds_data2[,c(77,78,79,80)])
  colnames(uds_hill2) <-c("H0", "H1", "H2", "day", "treatment", "time", "chamber")
  
  #convert NA to 0
  uds_hill2$H0[is.na(uds_hill2$H0)] <- 0 
  
  
  uds_hill3 <- subset(uds_hill2,chamber=="C2")
  uds_hill3$treatment[uds_hill3$treatment=="T"] <- "E"
  uds_hill3$treatment[uds_hill3$treatment=="W"] <- "D"
  #View(uds_hill2)
  if(q == 0) {
     g_box<-ggplot(uds_hill2, aes(y=H0, x=chamber))
     g_box<-g_box + geom_violin()
     g_box<-g_box + geom_boxplot(width=0.1) 
     g_box<-g_box + coord_fixed(ratio = 1/2)
     print(g_box)
     g_box2<-ggplot(uds_hill3, aes(y=H0, x=treatment))
     g_box2<-g_box2 + geom_violin()
     g_box2<-g_box2 + geom_boxplot(width=0.1) 
     g_box2<-g_box2 + coord_fixed(ratio = 1/2)
     print(g_box2)
     
     plotmeans(uds_hill2$H0~uds_hill2$chamber) 
     
     View(uds_hill2)
    
    #homogeneity of variance
    print(bartlett.test(uds_hill2$H0~uds_hill2$chamber))
    print(bartlett.test(uds_hill3$H0~uds_hill3$treatment))
    #One-way ANOVA with unequal-variance between chambers
    print(oneway.test(uds_hill2$H0~uds_hill2$chamber,var.equal=FALSE))
    print(oneway.test(uds_hill3$H0~uds_hill3$treatment,var.equal=FALSE))
    #Pairwise comparison with Welch's t test
    #print(pairwise.t.test(uds_hill2$H0, uds_hill2$treatment, p.adj="holm"))
  }
  if(q == 1) {
    g_box<-ggplot(uds_hill2, aes(y=H1, x=chamber))
     g_box<-g_box + geom_violin()
     g_box<-g_box + geom_boxplot(width=0.1) 
     g_box<-g_box + coord_fixed(ratio = 1/2)
     print(g_box)
     g_box2<-ggplot(uds_hill3, aes(y=H1, x=treatment))
     g_box2<-g_box2 + geom_violin()
     g_box2<-g_box2 + geom_boxplot(width=0.1) 
     g_box2<-g_box2 + coord_fixed(ratio = 1/2)
     print(g_box2)
     
     plotmeans(uds_hill2$H1~uds_hill2$chamber) 
     
     View(uds_hill2)
    
    #homogeneity of variance
    print(bartlett.test(uds_hill2$H1~uds_hill2$chamber))
    print(bartlett.test(uds_hill3$H1~uds_hill3$treatment))
    #One-way ANOVA with unequal-variance between chambers
    print(oneway.test(uds_hill2$H1~uds_hill2$chamber,var.equal=FALSE))
    print(oneway.test(uds_hill3$H1~uds_hill3$treatment,var.equal=FALSE))
    
  }
  if(q == 2) {
    g_box<-ggplot(uds_hill2, aes(y=H2, x=chamber))
     g_box<-g_box + geom_violin()
     g_box<-g_box + geom_boxplot(width=0.1) 
     g_box<-g_box + coord_fixed(ratio = 1/2)
     print(g_box)
     g_box2<-ggplot(uds_hill3, aes(y=H2, x=treatment))
     g_box2<-g_box2 + geom_violin()
     g_box2<-g_box2 + geom_boxplot(width=0.1) 
     g_box2<-g_box2 + coord_fixed(ratio = 1/2)
     print(g_box2)
     
     plotmeans(uds_hill2$H1~uds_hill2$chamber) 
     
     View(uds_hill2)
    
    #homogeneity of variance
    print(bartlett.test(uds_hill2$H2~uds_hill2$chamber))
    print(bartlett.test(uds_hill3$H2~uds_hill3$treatment))
    #One-way ANOVA with unequal-variance between chambers
    print(oneway.test(uds_hill2$H2~uds_hill2$chamber,var.equal=FALSE))
    print(oneway.test(uds_hill3$H2~uds_hill3$treatment,var.equal=FALSE))
    
  }
  
    
}

```
### Temporal dynamics of species richness
#### Observed species richness per plant individuals (Fig.S2)
```{r}
#all samplings are treated as independent, used in the manuscript 
####For Figure S2
alpha_hill_number(uds_data_all_day2, date="10d", q=0)
alpha_hill_number(uds_data_all_day2, date="32d", q=0)
alpha_hill_number(uds_data_all_day2, date="60d", q=0)

```
#### Summary of observerd richness
- There were significant differences bet.chambers (C1 < C2) at day 32 and day 60 but the difference is minor 
- There were no differences bet. D and E at any days

### Species richness estimates analyzed by iNEXT (Fig.3)
#### Loading data for iNEXT analyses
```{r}
incidence_whole <- read.csv("monthly-for-R_iNEXT.csv", header=T)  #for incidence data calculated by Excel
```
#### Estimating the coverage-based standardized richness (Fig.3a) & asymptotic richness (Fig.3b)
For standardized richness (Results for Fig.3a)\
```{r}
####This is code for graphics in the manuscript (Fig.3)
stat_test<-DataInfo(incidence_whole, datatype = "incidence")  ##This is necessary to find the minimum value of the coverage
estimateD(incidence_whole, datatype="incidence_freq", base="coverage", level=min(stat_test$SC))
```
For asymptotic richness (Results for Fig.3b)\
```{r}
#When Short-term repeated measurement is assumed to be independent
ChaoRichness(incidence_whole, datatype="incidence")  #It is OK to use dataframe, this is from iNEXT
```

#### Loading data for graphics
The above iNEXT results were manually summarized in excel and csv.\
```{r}
diversity_summary <- read.csv("iNEXT_summary.csv",header=T)
class(diversity_summary)
View(diversity_summary)
diversity_summary2<-diversity_summary[1:20,]
class(diversity_summary2$day)
```

#### Figure 3a
```{r}
#colours <- c("#F8766D", "#EE55B7", "#00BA38", "#619CFF", "#E76BF3")
g_div_a <- ggplot(data=diversity_summary2, aes(x=ID,y=SR))
g_div_a <- g_div_a + geom_errorbar(aes(ymax = Hconf, ymin = Lconf, colour=treatment), width = 0.0,size=5, alpha=0.5) #+ scale_colour_manual("Group", values = colours)
g_div_a <- g_div_a + geom_point(aes(colour=treatment), size=3) #+ scale_colour_manual("Group", values = colours)
g_div_a <- g_div_a +  geom_point(aes(y=Observed), size=3) + coord_fixed(ratio = 1/5) #+ scale_colour_manual("Group", values = colours) 
print(g_div_a)
```

#### Figure 3b
```{r}
#Plot Asymptotic richness (Chao index)(Fig3b)
diversity_summary3<-diversity_summary[16:20,]
diversity_summary3$ID <- c("1U", "2ED", "3E", "4D", "5UED")
#colours <- c("#F8766D", "#00BA38", "#619CFF", "#E76BF3", "#00B6EB", "#A58AFF", "#FB61D7")
g_div_b <- ggplot(data=diversity_summary3, aes(x=ID,y=AsyEST))
g_div_b <- g_div_b + geom_errorbar(aes(ymax = Hconf_asy, ymin = Lconf_asy, colour=treatment), width = 0.0,size=5, alpha=1.0) #+ scale_colour_manual("Group", values = colours)
g_div_b <- g_div_b + geom_point(size=5) 
print(g_div_b)
```

#### Summary for richness analysis
- Observed richness (Fig.3a)
  - the whole richness at each day was greater than any of single treatment or single chamber
  - the richness in the whole period was greater than any of single date for every treatment and chamber
- Standardized richness (Fig.3a)

- Asymptotic richness (Fig.3b)
  - Chamber 2 (E&D) had a greater richness than Chamber 1 (U)
  - Chamber 2 had a greater richness than D
  - The whole richness (UED) was greater than Chamber 1 and D
  
#### Species richness partitioning (Table S2, Fig. S3, and Fig.3C)
```{r}
partitioning_summary <- read.csv("richness_partitioning.csv",header=T)

g <- ggplot(partitioning_summary, aes(x = partition, y = richness, fill = fraction))
g <- g + geom_bar(stat = "identity")
g <- g + geom_text(aes(label = richness), position = position_stack(vjust = 0.5))
plot(g)

```


