###################################################
### ‘CONSTRAINED' FUNCTIONAL SIZE-BASED NETWORK ###
###################################################

# Load libraries
library(cluster) # to find groups in data
library(plyr) # to use the function "rename" and "revalue"
library (bipartite) # for bipartite networks


# With only numerical variables, we can perform the cluster analysis using the function "agnes". It computes agglomerative hierarchical clustering of the dataset.

# Load data
Traits <-read.csv("data.csv")
str(Traits) #explore the structure of the data


## (1) PLANTS
## (2) BEES
## (3) NETWORK


### (1) PLANTS ###

# Select the variable of interest (Nectar Holder Depth - NHD)
Plants <- subset(Traits, select = c(Interaction, Plant, Predicted.NHD))

## (1.1) Prepare dataset Plants for AGNES:
# Matrix or dataframe should be numeric. Our continuous variable "Perdicted.NHD" is standardized when applying the function "agnes"

# Remove the two first columns, which are not variables, they are the ID of the interaction and the plant species:
Plants_cluster <- Plants[-c(1,2)]

## (1.2) Function "AGNES: Agglomerative Nesting (hierarchical clustering)"

Plants_AGNES <- agnes(Plants_cluster, diss = FALSE,
                      metric = "euclidean", stand = TRUE, method = "average",
                      keep.diss = FALSE,
                      keep.data = TRUE,
                      trace.lev = 0)

# Visualize the hierarchical tree and set the number of K clusters
TREE_Plants_AGNES <- pltree(Plants_AGNES) 

Plants.AGNES_K <- cutree(as.hclust(Plants_AGNES), h = 0.3) # "cuttree" is a function that cuts a tree into several groups either by specifying the desired number(s) of groups or the cut height(s). By setting "h=0.3", we obtain 10 K (as in the species matrix).


# Cbind cluster assignation to the database "Plants":
Plants_K_agnes <- data.frame(cbind(Plants,Plants.AGNES_K))

# Create a .csv table with results:
write.table(Plants_K_agnes, "Plants.agnes_NHD.K.csv", sep=",", dec=".")


### (2) BEES ###

# Select the variable of interest from the original dataset (Intertegular Distance - IT)
Bees <- subset(Traits, select = c(Interaction, Bee, IT.distance))

## (2.1) Prepare dataset Plants for AGNES
# Matrix or dataframe should be numeric. Our continuous (numeric) variable "IT.distance" is standardized when applying the function "agnes"

# Remove the two first columns, which are not variables, they are the ID of the interaction and bee species
Bees_cluster <- Bees[-c(1,2)]

## (2.2) Function "AGNES: Agglomerative Nesting (hierarchical clustering)"

Bees_AGNES <- agnes(Bees_cluster, diss = FALSE,
                    metric = "euclidean", stand = TRUE, method = "average",
                    keep.diss = FALSE,
                    keep.data = TRUE,
                    trace.lev = 0)

# Visualize the hierarchical tree and set the number of K clusters
TREE_Bees_AGNES <- pltree(Bees_AGNES) 

Bees.AGNES_K <- cutree(as.hclust(Bees_AGNES), h = 0.15) # By setting "h=0.15", we obtain 28 K (as the species matrix).

# Cbind cluster assignation to the database "Bees":
Bees_K_agnes <- data.frame(cbind(Bees,Bees.AGNES_K))

# Create a .csv table with results:
write.table(Bees_K_agnes, "Bees.agnes_IT.K.csv", sep=",", dec=".") 

### (3) NETWORK ###

## (3.1) Interaction Matrix
# Merge datasets
Plants.Bees <- merge(Plants_K_agnes, Bees_K_agnes, by = "Interaction")

# Encode the variables with cluster assignation as factor:
Plants.Bees$Plants.AGNES_K <- factor(Plants.Bees$Plants.AGNES_K)
Plants.Bees$Bees.AGNES_K <- factor(Plants.Bees$Bees.AGNES_K)

# Rename values in "Plants.AGNES_K". Order from the smallest (F01) to bigger size (F10):
Plants.Bees$Plants.AGNES_K <- revalue(Plants.Bees$Plants.AGNES_K, c("1"="F10", "2"="F03","3"="F05", "4"="F01", "5"="F04", "6"="F08", "7"="F06", "8"="F02", "9"="F07", "10"="F09"))

Plants.Bees$Plants.AGNES_K <- factor(Plants.Bees$Plants.AGNES_K,levels=c("F01","F02","F03","F04","F05","F06","F07","F08","F09", "F10"),ordered=TRUE)

# Rename values in variable "Bees.AGNES_K". Order from the smallest (B01) to bigger size (B28)
Plants.Bees$Bees.AGNES_K <- revalue(Plants.Bees$Bees.AGNES_K, c("1"="B16", "2"="B18","3"="B13", "4"="B05", "5"="B12", "6"="B20", "7"="B10", "8"="B21", "9"="B09", "10"="B14", "11"="B19", "12"="B17","13"="B15", "14"="B03", "15"="B11", "16"="B04", "17"="B23", "18"="B06", "19"="B07", "20"="B08", "21"="B27", "22"="B28","23"="B22", "24"="B01", "25"="B02", "26"="B24", "27"="B26", "28"="B25")) 

Plants.Bees$Bees.AGNES_K <- factor(Plants.Bees$Bees.AGNES_K,levels=c("B01","B02","B03","B04","B05","B06","B07","B08","B09","B10","B11","B12","B13","B14","B15","B16","B17","B18","B19","B20","B21","B22","B23","B24","B25","B26","B27","B28"),ordered=TRUE)

# Create a table with the frequency of interactions
Plants.Bees_K_AGNES <- with(Plants.Bees, table(Plants.AGNES_K,Bees.AGNES_K)) 
# Convert table to matrix
Plants.Bees_K_AGNES <-as.data.frame.matrix(Plants.Bees_K_AGNES) 
# Create a .csv with results
write.table(Plants.Bees_K_AGNES, "NHD.IT_agnes.K.csv", sep=",", dec=".")

## (3.2) Calculate network parameters

# Network level:
NHD_IT_K.agnes_prop <- networklevel(Plants.Bees_K_AGNES, index=c("weighted connectance", "weighted NODF", "interaction evenness", "H2"))

write.table(NHD_IT_K.agnes_prop, "NHD_IT_K.agnes_NetworkParameters.csv", sep=",", dec=".")

# Node level:

NHD_IT_K.agnes_node_prop <- specieslevel(Plants.Bees_K_AGNES, index="ALLBUTD", level="both")

NHD_IT_K.agnesHL_node_prop <- as.data.frame(NHD_IT_K.agnes_node_prop$`higher level`)

NHD_IT_K.agnesLL_node_prop <- as.data.frame(NHD_IT_K.agnes_node_prop$`lower level`)

write.table(NHD_IT_K.agnesHL_node_prop, "NHD_IT_K.agnes_HLParameters.csv", sep=",", dec=".")
write.table(NHD_IT_K.agnesLL_node_prop, "NHD_IT_K.agnes_LLParameters.csv", sep=",", dec=".")

## (3.3) Visualize network

visweb(Plants.Bees_K_AGNES, type="as is", text="interaction")

plotweb(Plants.Bees_K_AGNES, method ="normal", bor.col.interaction ="gray60", text.rot=90, labsize = 0.8, y.lim=c(-1,2.5))



######################################################
### ‘UNCONSTRAINED' FUNCTIONAL SIZE-BASED NETWORK ###
#####################################################

# Load libraries
library(NbClust) # to find the optimal numbers of groups in data set
library(plyr) # to use the function "rename" and "revalue"
library (bipartite) # for bipartite networks

# NbClust provides 30 indexes for determining the optimal number of clusters in a data set and offers the best clustering scheme from different results to the user.

# Load data
Traits <-read.csv("data.csv")
str(Traits) #explore the structure of the data


### (1) PLANTS
### (2) BEES
### (3) NETWORK


### (1) PLANTS

# Select the variable of interest (Nectar Holder Depth - NHD)
Plants <- subset(Traits, select = c(Interaction, Plant, Predicted.NHD))

## (1.1) Prepare dataset Plants for NbClust:

# Matrix or dataframe should be numeric
# Remove the two first columns, which are the ID of the interaction and the plant species
Plants_cluster <- Plants[-c(1,2)]

## (1.2) Apply the function "NbClust":
Plants_NbClust <- NbClust(Plants_cluster, diss = NULL, distance = "euclidean", min.nc=2, max.nc=68, method = "average", index = "all")

# If the matrix is a dissimilarity matrix, diss = (name of the diss matrix). 
# min.nc = minimum number of clusters
# max.nc = maximum number of clusters (it should be a sub-multiple or multiple of the number of rows)
# Method = average <- uses the average pairwise distance between all pairs of points in the different clusters as the measure of distance. The clustering for average linkage usually lies somewhere between the solutions given by single and complete linkage.
# index = to specify the index wanted, or "all" in this case

# View result as matrix:
Plants.All.index <- as.matrix(Plants_NbClust$All.index)
Plants.Best.nc <- as.matrix(Plants_NbClust$Best.nc)

# cbind best partitions to the plants dataset:
Plants.Best.partition <- data.frame(cbind(Plants,Plants_NbClust$Best.partition))

# Create a .csv table with results
write.table(Plants.Best.partition, "Plants.NbClust_NHD.K.csv", sep=",", dec=".")

### (2) BEES

# Select the variable of interest from the original dataset (Intertegular Distance - IT)
Bees <- subset(Traits, select = c(Interaction, Bee, IT.distance))

## (2.1) Prepare dataset Bees for NbClust:

# Matrix or dataframe should be numeric
# Remove the two first columns, which are the ID of the interaction and bee species
Bees_cluster <- Bees[-c(1,2)]

## (2.2) Apply the function "NbClust":

Bees_NbClust <- NbClust(Bees_cluster, diss = NULL, distance = "euclidean", min.nc=2, max.nc=136, method = "average", index = "all")

# View result as matrix:
Bees.All.index <- as.matrix(Bees_NbClust$All.index)
Bees.Best.nc <- as.matrix(Bees_NbClust$Best.nc)

# cbind best partitions to the bees dataset:
Bees.Best.partition <- data.frame(cbind(Bees,Bees_NbClust$Best.partition))

# Create a .csv table with results
write.table(Bees.Best.partition, "Bees.NbClust_IT.K.csv", sep=",", dec=".")

### (3) NETWORK

## (3.1) Interaction Matrix:

# Merge datasets
Plants.Bees <- merge(Plants.Best.partition, Bees.Best.partition, by = "Interaction")

# Encode the variables with cluster assignation as factor:
Plants.Bees$Plants_NbClust.Best.partition <- factor(Plants.Bees$Plants_NbClust.Best.partition)
Plants.Bees$Bees_NbClust.Best.partition <- factor(Plants.Bees$Bees_NbClust.Best.partition)

# Rename values in variable "Plants_NbClust.Best.partition" and order
Plants.Bees$Plants_NbClust.Best.partition <- revalue(Plants.Bees$Plants_NbClust.Best.partition, c("1"="F05", "2"="F02","3"="F03", "4"="F01", "5"="F04"))

Plants.Bees$Plants_NbClust.Best.partition <- factor(Plants.Bees$Plants_NbClust.Best.partition,levels=c("F01","F02","F03","F04","F05"),ordered=TRUE)

# Rename values in variable "Bees_NbClust.Best.partition" and order
Plants.Bees$Bees_NbClust.Best.partition <- revalue(Plants.Bees$Bees_NbClust.Best.partition, c("1"="B03", "2"="B02","3"="B01", "4"="B04"))

Plants.Bees$Bees_NbClust.Best.partition <- factor(Plants.Bees$Bees_NbClust.Best.partition,levels=c("B01","B02","B03","B04"),ordered=TRUE)

# Create a table with frequency of interactions
Plants.Bees_K_NbClust <- with(Plants.Bees, table(Plants_NbClust.Best.partition,Bees_NbClust.Best.partition))

# Convert table to matrix
Plants.Bees_K_NbClust <-as.data.frame.matrix(Plants.Bees_K_NbClust)

# Create a .csv with results
write.table(Plants.Bees_K_NbClust, "NHD.IT_matrix_NbClust.csv", sep=",", dec=".")

## (3.2) Network parameters:

#Network level
NHD_IT_K.NbClust_prop <- networklevel(Plants.Bees_K_NbClust, index=c("weighted connectance", "weighted NODF", "interaction evenness", "H2"))

write.table(NHD_IT_K.NbClust_prop, "NHD_IT_K.NbClust_NetworkParameters.csv", sep=",", dec=".")

# Node level:
NHD_IT_K.NbClust_node_prop <- specieslevel(Plants.Bees_K_NbClust, index="ALLBUTD", level="both")
NHD_IT_K.NbClustHL_node_prop <- as.data.frame(NHD_IT_K.NbClust_node_prop$`higher level`)
NHD_IT_K.NbClustLL_node_prop <- as.data.frame(NHD_IT_K.NbClust_node_prop$`lower level`)

write.table(NHD_IT_K.NbClustHL_node_prop, "NHD_IT.NbClust_HLParameters.csv", sep=",", dec=".")
write.table(NHD_IT_K.NbClustLL_node_prop, "NHD_IT.NbClust_LLParameters.csv", sep=",", dec=".")

## (3.3) Visualize network
visweb(Plants.Bees_K_NbClust, type="as is", text="interaction")
plotweb(Plants.Bees_K_NbClust, method ="normal", bor.col.interaction ="gray60", text.rot=90, labsize = 0.8, y.lim=c(-1,2.5))