diff --git a/.Rhistory b/.Rhistory
new file mode 100644
index 0000000..3dbfddd
--- /dev/null
+++ b/.Rhistory
@@ -0,0 +1,512 @@
+### Plot
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T)
+tiff("figures/morph.tiff", width = 24, height = 10, units = "cm", res = 300)
+ggplot(ggmorph, aes(x=module, y=measure, fill=variable)) +
+ylab("Measure (mm)")+ xlab("Modules")+ ylim(0, 150) +
+scale_fill_manual(values=c("slateblue", "goldenrod1", "coral1"))+
+geom_boxplot(width=0.5, color="black", position = position_dodge(width=0.5)) +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", size =20),
+axis.title.x = element_text(color="black", size =20),
+axis.text= element_text(color="black", size=19), legend.position = "none")
+dev.off()
+################################################################################
+##### SAMPLING COMPLETENESS ANALYSIS
+################################################################################
+# Loading interaction data for Chao1 estimator
+sampbat<- read.xls("data/sampbat.xlsx", h=T)
+estimateR(sampbat, index =c("chao"))
+str(sampbat)
+samphawk <- read.xls("data/samphawk.xlsx", h=T)
+estimateR(samphawk, index =c("chao"))
+str(samphawk)
+# Loading interaction data for rarefaction curve drawing
+sampling_bats <- read.xls("data/sampling_bats.xlsx", h=T)
+curve_bat<- specaccum(sampling_bats, method="rarefaction")
+sampling_hawkmoths <- read.xls("data/sampling_hawkmoths.xlsx", h=T)
+curve_hawk<- specaccum(sampling_hawkmoths, method="rarefaction")
+# Plot curves
+tiff("figures/sampling.tiff", width = 20, height = 20, units = "cm", res = 600)
+par(mfrow=c(1,2))
+plot(curve_hawk, ci.type = "poly", xvar = "individuals", ci.lty=0, ylab = NA,
+xlab = NA,
+ci.col=rgb(0.7, 0, 0.2, 0.3), ylim=c(0,30))
+abline(h=22.2, lty=1, col=rgb(0.7, 0, 0.2, 0.3), lwd=2.5)
+abline(h=(22.2+0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+abline(h=(22.2-0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+plot(curve_bat, ci.type = "poly", xvar = "individuals", ci.lty=0,
+ci.col=rgb(0, 0, 0.5, 0.3), ylim=c(0,30), ylab=NA, xlab=NA)
+abline(h=14, lty=1, col=rgb(0, 0, 0.5, 0.3), lwd=2.5)
+abline(h=(14-2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+abline(h=(14+2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+dev.off()
+### Plot
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T)
+### Plot
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T)
+ggplot(ggmorph, aes(x=module, y=measure, fill=variable)) +
+ylab("Measure (mm)")+ xlab("Modules")+ ylim(0, 150) +
+scale_fill_manual(values=c("slateblue", "goldenrod1", "coral1"))+
+geom_boxplot(width=0.5, color="black", position = position_dodge(width=0.5)) +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", size =20),
+axis.title.x = element_text(color="black", size =20),
+axis.text= element_text(color="black", size=19), legend.position = "none")
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T)
+tiff("figures/morph.tiff", width = 24, height = 10, units = "cm", res = 300)
+ggplot(ggmorph, aes(x=module, y=measure, fill=variable)) +
+ylab("Measure (mm)")+ xlab("Modules")+ ylim(0, 150) +
+scale_fill_manual(values=c("slateblue", "goldenrod1", "coral1"))+
+geom_boxplot(width=0.5, color="black", position = position_dodge(width=0.5)) +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", size =20),
+axis.title.x = element_text(color="black", size =20),
+axis.text= element_text(color="black", size=19), legend.position = "none")
+dev.off()
+################################################################################
+################################################################################
+#Supplement to the paper Queiroz et al. (2020, Biotropica).
+#Ecological Synthesis Lab (SintECO): https://marcomellolab.wordpress.com.
+#Authors: Joel A. Queiroz, Ugo M. Diniz, Diego P. Vázquez, Zelma M. Quirino,
+#Francisco A.R. Santos, Marco A.R. Mello, Isabel C. Machado.
+#See README for further info:
+#https://github.com/marmello77/queiroz_et_al_2020/blob/main/README.md
+################################################################################
+################################################################################
+################################################################################
+##### SET THE STAGE
+################################################################################
+#Set the working directory
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+#Delete all previous objects
+rm(list= ls())
+#Load the required packages
+library(igraph)
+library(bipartite)
+library(Rmisc)
+library(vegan)
+library(gdata)
+library(ggplot2)
+#Load some custom-made functions
+source("RestNullModel.R")
+source("PosteriorProb.R")
+source("MyTriangle.R")
+################################################################################
+##### PROCESS THE NETWORK
+################################################################################
+#Import the network
+data <- as.matrix(read.delim("data/network.txt", row.names=1))
+#Inspect the network
+class(data)
+data
+dim(data)
+min(data)
+max(data)
+#Plot the matrix
+visweb(data)
+#Convert the network to igraph format
+data2 <- graph_from_incidence_matrix(data, directed = F, weighted = TRUE)
+#Inspect object
+class(data2)
+data2
+E(data2)
+V(data2)$name
+#Inform which nodes represent which taxonomic groups
+V(data2)$set[1:nrow(data)] = c("Moths",	"Moths", "Bats",	"Bats",	"Moths",	"Moths",
+"Moths",	"Moths",	"Moths",	"Bats",	"Bats",	"Moths",
+"Moths",	"Moths",	"Moths",	"Moths",	"Moths",
+"Moths",	"Moths")
+V(data2)$set[(nrow(data)+1):(nrow(data)+ncol(data))] = "Plants"
+################################################################################
+##### DRAW THE NETWORK
+################################################################################
+#Set layout
+lay1 <- layout_nicely(data2)
+#Set edge curvatures
+curves1 = curve_multiple(data2)
+#Set edge mode and width
+E(data2)$arrow.mode = 0
+E(data2)$width = E(data2)$weight/5+1
+#Calculate Louvain modularity (resolution = 1.0, similar to DIRT_LPA+)
+data2.lou = cluster_louvain(data2)
+#Import "diamond" vertex shape
+source("MyTriangle.R")
+#Set vertex shapes
+V(data2)$shape = V(data2)$set
+V(data2)$shape = gsub("Bats","diamond",V(data2)$shape)
+V(data2)$shape = gsub("Moths","square",V(data2)$shape)
+V(data2)$shape = gsub("Plants","circle",V(data2)$shape)
+##Set node and cloud colors by modularity
+colrs <- rainbow(length(data2.lou), alpha = 1.0, s = 1, v = 0.8)
+V(data2)$color <- colrs[data2.lou$membership]
+clouds = rainbow(length(data2.lou), alpha = 0.1)
+#Plot and export the graph
+tiff(filename= "figures/network.tif", res= 300, height= 3000, width= 3100)
+par(mfrow=c(1,1),mar=c(1,1,1,5))
+plot(data2.lou,
+data2,
+col = V(data2)$color,
+mark.border="lightgrey",
+mark.col=clouds,
+vertex.size=7.5,
+vertex.label=V(data2)$name,
+vertex.label.color="white",
+vertex.label.cex=.3,
+edge.color = adjustcolor("grey", alpha.f = .5),
+edge.curved=0.3,
+edge.width = 3,
+layout=lay1)
+legend(x = 0.9,y = 1.0, legend = c("Bats", "Moths", "Plants"),
+pch = c(18,15,19),  title="Taxon",
+text.col = "gray20", title.col = "black",
+box.lwd = 0, cex = 2, col=c("grey", "grey", "grey"))
+par(mfrow=c(1,1))
+dev.off()
+################################################################################
+##### NETWORK LEVEL ANALYSIS
+################################################################################
+#Set seed
+set.seed(14)
+#Set the number of permutations to be used in all null model analyses
+permutations <- 1000
+#Generate randomized matrices
+nulls <- nullmodel(data, N=permutations, method="vaznull")
+##### MODULARITY
+#Calculate metric for the original network
+Mod <- computeModules(data, method = "Beckett")
+Mod@likelihood
+#Extract module membership
+Part <- bipartite::module2constraints(Mod)
+row.Part <- Part[1:nrow(data)]
+col.Part <- Part[(nrow(data)+1):(nrow(data)+ncol(data))]
+#Calculate metric for the randomized networks
+nullmod <- sapply(nulls, computeModules, method = "Beckett")
+modnull <- sapply(nullmod, function(x) x@likelihood)
+(Mod@likelihood - mean(modnull))/sd(modnull) # Z value
+Mod.sig <- sum(modnull>(Mod@likelihood)) / length(modnull) # p value
+Mod.sig
+#Plot the observed value against the distribution of randomized values
+png("specialization.png", width = 3000, height = 3000, res = 300)
+plot(density(modnull), main="Observed vs. randomized",
+xlim=c(min((Mod@likelihood), min(modnull)),
+max((Mod@likelihood), max(modnull))))
+abline(v=Mod@likelihood, col="red", lwd=2, xlab="")
+dev.off()
+Mod@likelihood #observed
+mean(modnull) #randomized mean
+sd(modnull) #randomized SD
+(Mod@likelihood - mean(modnull))/sd(modnull) # Z-value
+sum(modnull>(Mod@likelihood)) / length(modnull) #randomized > observed
+sum(modnull<(Mod@likelihood)) / length(modnull) #randomized < observed
+##### SPECIALIZATION
+#Calculate metric for the original network
+Spec <- networklevel(data, index="H2")
+class(Spec)
+#Calculate metric for the randomized networks
+randomized.Spec <- unlist(sapply(nulls, networklevel, index="H2"))
+(Spec - mean(randomized.Spec))/sd(randomized.Spec) # Z value
+Spec.sig <- sum(randomized.Spec>Spec)/length(randomized.Spec) # p value
+Spec.sig
+#Plot the observed value against the distribution of randomized values
+png("specialization.png", width = 3000, height = 3000, res = 300)
+plot(density(randomized.Spec), main="Observed vs. randomized",
+xlim=c(min((Spec), min(randomized.Spec)),
+max((Spec), max(randomized.Spec))))
+abline(v=Spec, col="red", lwd=2, xlab="")
+dev.off()
+Spec #observed
+mean(randomized.Spec) #randomized mean
+sd(randomized.Spec) #randomized SD
+(Spec - mean(randomized.Spec))/sd(randomized.Spec) # Z-value
+sum(randomized.Spec>(Spec)) / length(randomized.Spec) #randomized > observed
+sum(randomized.Spec<(Spec)) / length(randomized.Spec) #randomized < observed
+##### NESTEDNESS
+#Calculate metric for the original network
+Nest <- networklevel(data, index="weighted NODF")
+#Calculate metric for the randomized networks
+randomized.Nest <- unlist(sapply(nulls, networklevel, index="weighted NODF"))
+(Nest - mean(randomized.Nest))/sd(randomized.Nest) # Z value
+Nest.sig <- sum(randomized.Nest>Nest)/length(randomized.Nest) # p value
+Nest.sig
+#Plot the observed value against the distribution of randomized values
+png("nestedness.png", width = 3000, height = 3000, res = 300)
+plot(density(randomized.Nest), main="Observed vs. randomized",
+xlim=c(min((Nest), min(randomized.Nest)),
+max((Nest), max(randomized.Nest))))
+abline(v=Nest, col="red", lwd=2, xlab="")
+dev.off()
+Nest #observed
+mean(randomized.Nest) #randomized mean
+sd(randomized.Nest) #randomized SD
+(Nest - mean(randomized.Nest))/sd(randomized.Nest) # Z-value
+sum(randomized.Nest>(Nest)) / length(randomized.Nest) #randomized > observed
+sum(randomized.Nest<(Nest)) / length(randomized.Nest) #randomized < observed
+##### COMPOUND TOPOLOGY
+#Calculate the desired nestedness metric (here WNODA) for the original network.
+obs.com <- unlist(bipartite::nest.smdm(x = data,
+constraints = Part, #Input the modular structured recovered from step 2
+weighted = T, #By considering the edge weights, you are choosing WNODA
+decreasing = "abund"))
+#Check the scores
+obs.com
+#Calculate constrained interaction probabilities considering the network's modular structure
+Pij <- PosteriorProb(M = data,
+R.partitions = row.Part, C.partitions = col.Part, #Input the modular structured recovered from step 2
+Prior.Pij = "degreeprob", #Choose the null model
+Conditional.level = "modules") #Choose the kind of constraints
+#Check what those probabilities look like
+Pij
+dim(Pij)
+#Generate randomized networks with the null model of your choice, considering the interaction probabilities calculated before.
+nulls.com <- RestNullModel(M = data,
+Pij.Prob = Pij, #Recover the probabilities calculated in the previous command
+Numbernulls = permutations, #This step may take long, so start experimenting with low values
+Print.null = T,
+allow.degeneration = F, #Choose whether you allow orphan rows and columns to be removed or not
+return.nonrm.species = F,
+connectance = T, byarea = T,
+R.partitions = row.Part,
+C.partitions = col.Part)
+#Calculate the nestedness within and between modules
+rest.nest <- nest.smdm(data, constraints = Part,
+weighted = T,
+decreasing = "abund",
+sort = T)
+unlist(rest.nest)
+null.com <- sapply(nulls.com,
+function(x) bipartite::nest.smdm(x = x,
+constraints = Part,
+weighted = T,
+decreasing = "abund"))
+WNODA.null.com <- unlist(null.com[3,])
+WNODAsm.null.com <- unlist(null.com[8,])
+WNODAdm.null.com <- unlist(null.com[9,])
+#Plot the observed nestedness value against the distribution of randomized values
+par(mfrow = c(1,3))
+plot(density(WNODA.null.com), xlim=c(min(obs.com[3], min(WNODA.null.com)),
+max(obs.com[3], max(WNODA.null.com))),
+main="observed vs. randomized", xlab = "WNODA matrix")
+abline(v=obs.com[3], col="red", lwd=2)
+plot(density(WNODAsm.null.com), xlim=c(min(obs.com[8], min(WNODAsm.null.com)),
+max(obs.com[8], max(WNODAsm.null.com))),
+main="observed vs. randomized", xlab = "WNODAsm matrix")
+abline(v=obs.com[8], col="red", lwd=2)
+plot(density(WNODAdm.null.com), xlim=c(min(obs.com[9], min(WNODAdm.null.com)),
+max(obs.com[9], max(WNODAdm.null.com))),
+main="observed vs. randomized", xlab = "WNODAdm matrix")
+abline(v=obs.com[9], col="red", lwd=2)
+par(mfrow = c(1,1))
+#Estimate the p-values
+#Nestedness in th entire network
+praw.WNODA <- sum(WNODA.null.com>obs.com[3]) / length(WNODA.null.com)
+p.WNODA <- ifelse(praw.WNODA > 0.5, 1- praw.WNODA, praw.WNODA)    # P-value
+p.WNODA
+#Nestedness within the modules
+praw.WNODAsm <- sum(WNODAsm.null.com>obs.com[8]) / length(WNODAsm.null.com)
+p.WNODAsm <- ifelse(praw.WNODAsm > 0.5, 1- praw.WNODAsm, praw.WNODAsm)    # P-value
+p.WNODAsm
+#Nestedness between the modules
+praw.WNODAdm <- sum(WNODAdm.null.com>obs.com[9]) / length(WNODAdm.null.com)
+p.WNODAdm <- ifelse(praw.WNODAdm > 0.5, 1- praw.WNODAdm, praw.WNODAdm)    # P-value
+p.WNODAdm
+##### PLOT THE COMPOUND TOPOLOGY
+par(mfrow = c(1,1))
+#Sort the matrix in a way that facilitates visualizing the compound topology
+data.comp <- bipartite::sortmatrix(matrix = data, topology = "compound",
+sort_by = "weights",
+row_partitions = row.Part,
+col_partitions = col.Part)
+#Assign colors for the modules
+modcol <- rainbow((length(unique(Part))), alpha=1, s = 1, v = 1)
+#Plot the matrix
+png("figures/compound.png", width = 3000, height = 3000, res = 300)
+plotmatrix(data.comp$matrix,
+row_partitions = data.comp$row_partitions,
+col_partitions = data.comp$col_partitions,
+border = T,
+binary = F,
+modules_colors = modcol,
+within_color = modcol,
+between_color = "lightgrey")
+dev.off()
+##### EXPORT A SUMMARY OF THE TOPOLOGICAL RESULTS
+sink(file = "results/results_topology.txt")
+paste("Topological analysis of the nocturnal pollination network")
+paste("Queiroz et al. 2020, Biotropica")
+cat("\n")
+paste("The network has", nrow(data), "rows and", ncol(data), "columns.")
+cat("\n")
+paste("The network's specialization (H2) is", round(Spec, 2),",", "P =", round(Spec.sig, 2))
+cat("\n")
+paste("The network's modularity (DIRT_LPA+) is", round(Mod@likelihood, 2), ",", "P =", round(Mod.sig, 2), ",", "and it contains", length(unique(Part)),"modules.")
+cat("\n")
+paste("The network's nestedness (WNODF) is", round(Nest/100, 2),",", "P =", round(Nest.sig, 2))
+cat("\n")
+paste("The network shows the following scores of nestedness (WNODA):")
+cat("\n")
+paste("Entire network =", round(rest.nest$WNODAmatrix/100, 2), ",", "P =", round(p.WNODA, 2))
+cat("\n")
+paste("Between the modules =", round(rest.nest$WNODA_DM_matrix/100, 2), ",", "P =", round(p.WNODAdm, 2))
+cat("\n")
+paste("Within the modules =", round(rest.nest$WNODA_SM_matrix/100, 2), ",", "P =", round(p.WNODAsm, 2))
+cat("\n")
+sink(file = NULL, )
+################################################################################
+##### SPECIES LEVEL ANALYSIS (CENTRALITY)
+################################################################################
+# Specialization (d')
+d <- specieslevel(data,index="d")
+dplants <- d$`higher level`
+write.csv(dplants, "results/dplants.csv") # writing a separate csv file
+# Betweenness centrality (BC)
+BC <- specieslevel(data, index="betweenness")
+BCplants <- BC$higher
+write.csv(BCplants, "results/BCplants.csv") # writing a separate csv file
+# Normalized degree (nk)
+ND <-ND(data, normalised=T)
+NDplants <- ND$higher
+write.csv(NDplants, "results/NDplants.csv") # writing a separate csv file
+##### Plot centrality metrics by syndrome
+#Import data
+plants <- read.xls("data/plants.xlsx", h=T) # reading compiled spreadsheet with species & metrics classified by guild
+# Change reference level for GLMs
+ord <-  ordered(plants$Guild, levels = c("sphin", "chiro", "other"))
+# BC
+tiff("figures/bc.tiff", width = 12, height = 15, units = "cm", res = 300)
+ggplot(plants, aes(x=ord, y=bc, fill=Guild)) +
+ylab("BC")+ xlab("")+ ylim(0, 0.15) +
+scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+geom_boxplot(width=0.5, color="black") +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", face ="italic", size =23),
+axis.text= element_text(color="black", size=19),
+legend.position = "none")
+dev.off()
+# d'
+tiff("figures/d.tiff", width = 12, height = 15, units = "cm", res = 300)
+ggplot(plants, aes(x=ord, y=d, fill=Guild)) +
+ylab("d'")+ xlab("")+ ylim(0, 0.8) +
+scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+geom_boxplot(width=0.5, color="black") +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", face ="italic", size =23),
+axis.text= element_text(color="black", size=19),
+legend.position = "none")
+dev.off()
+# ND
+tiff("figures/nk.tiff", width = 12, height = 15, units = "cm", res = 300)
+ggplot(plants, aes(x=ord, y=nk, fill=Guild)) +
+ylab("nk")+ xlab("")+ ylim(0, 1.1) +
+scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+geom_boxplot(width=0.5, color="black") +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", face ="italic", size =23),
+axis.text= element_text(color="black", size=19),
+legend.position = "none")
+dev.off()
+##### Run GLMs to compare centrality scores
+table(plants$Guild)
+plants$Guild <- relevel(plants$Guild, ref="chiro") #changing reference level for GLMs
+plants$Guild <- relevel(plants$Guild, ref="sphin")
+plants$Guild <- relevel(plants$Guild, ref="other")
+# d'
+glmd <- glm(plants$d ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmd)
+glm_d <- anova(glmd, test = "Chisq")
+glm_d
+## BC
+glmbc <- glm(plants$bc ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmbc)
+glm_bc <- anova(glmbc, test = "Chisq")
+glm_bc
+## nk
+glmnk <- glm(plants$nk ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmnk)
+glm_k <- anova(glmnk, test = "Chisq")
+glm_k
+#Export the results
+sink(file = "results/results_centrality.txt")
+paste("Comparison of centrality by syndrome")
+paste("Queiroz et al. 2020, Biotropica")
+cat("\n")
+paste(capture.output(glmd))
+paste(capture.output(glm_d))
+cat("\n")
+paste(capture.output(glmbc))
+paste(capture.output(glm_bc))
+cat("\n")
+paste(capture.output(glmnk))
+paste(capture.output(glm_k))
+sink(file = NULL, )
+################################################################################
+##### MORPHOMETRIC ANALYSIS
+################################################################################
+##### Import morphology data
+morph_plants<-read.xls("data/morph_pla.xlsx", h=T)
+morph_pol <- read.xls("data/morph_pol.xlsx", h=T)
+str(morph_plants)
+str(morph_pol)
+# Change reference level for GLMs
+morph_plants$module <- relevel(morph_plants$module, ref="bat")
+morph_pol$module <- relevel(morph_pol$module, ref="hawk1")
+##### Run GLMs to compare modules
+# Pollinator tongues
+glm_pol <- glm(morph_pol$length_pol~morph_pol$module, family=gaussian())
+summary(glm_pol)
+anova(glm_pol, test = "Chisq")
+# Floral width (w) and length (l)
+glm_pla_l <- glm(morph_plants$length_pla~morph_plants$module, family=gaussian())
+glm_pla_w <- glm(morph_plants$width_pla~morph_plants$module, family=gaussian())
+summary(glm_pla_l)
+anova(glm_pla_l, test = "Chisq")
+summary(glm_pla_w)
+anova(glm_pla_w, test = "Chisq")
+### Plot
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T)
+tiff("figures/morph.tiff", width = 24, height = 10, units = "cm", res = 300)
+ggplot(ggmorph, aes(x=module, y=measure, fill=variable)) +
+ylab("Measure (mm)")+ xlab("Modules")+ ylim(0, 150) +
+scale_fill_manual(values=c("slateblue", "goldenrod1", "coral1"))+
+geom_boxplot(width=0.5, color="black", position = position_dodge(width=0.5)) +
+theme_classic() +
+theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+axis.title.y = element_text(color="black", size =20),
+axis.title.x = element_text(color="black", size =20),
+axis.text= element_text(color="black", size=19), legend.position = "none")
+dev.off()
+################################################################################
+##### SAMPLING COMPLETENESS ANALYSIS
+################################################################################
+# Loading interaction data for Chao1 estimator
+sampbat<- read.xls("data/sampbat.xlsx", h=T)
+estimateR(sampbat, index =c("chao"))
+str(sampbat)
+samphawk <- read.xls("data/samphawk.xlsx", h=T)
+estimateR(samphawk, index =c("chao"))
+str(samphawk)
+# Loading interaction data for rarefaction curve drawing
+sampling_bats <- read.xls("data/sampling_bats.xlsx", h=T)
+curve_bat<- specaccum(sampling_bats, method="rarefaction")
+sampling_hawkmoths <- read.xls("data/sampling_hawkmoths.xlsx", h=T)
+curve_hawk<- specaccum(sampling_hawkmoths, method="rarefaction")
+# Plot curves
+tiff("figures/sampling.tiff", width = 20, height = 20, units = "cm", res = 600)
+par(mfrow=c(1,2))
+plot(curve_hawk, ci.type = "poly", xvar = "individuals", ci.lty=0, ylab = NA,
+xlab = NA,
+ci.col=rgb(0.7, 0, 0.2, 0.3), ylim=c(0,30))
+abline(h=22.2, lty=1, col=rgb(0.7, 0, 0.2, 0.3), lwd=2.5)
+abline(h=(22.2+0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+abline(h=(22.2-0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+plot(curve_bat, ci.type = "poly", xvar = "individuals", ci.lty=0,
+ci.col=rgb(0, 0, 0.5, 0.3), ylim=c(0,30), ylab=NA, xlab=NA)
+abline(h=14, lty=1, col=rgb(0, 0, 0.5, 0.3), lwd=2.5)
+abline(h=(14-2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+abline(h=(14+2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+dev.off()
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..9bea433
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+
+.DS_Store
diff --git a/MyTriangle.R b/MyTriangle.R
new file mode 100644
index 0000000..26b85b9
--- /dev/null
+++ b/MyTriangle.R
@@ -0,0 +1,21 @@
+MyTriangle <- function(coords, v=NULL, params) {
+  vertex.color <- params("vertex", "color")
+  if (length(vertex.color) != 1 && !is.null(v)) {
+    vertex.color <- vertex.color[v]
+  }
+  vertex.frame.color <- params("vertex", "frame.color")
+  if (length(vertex.frame.color) != 1 && !is.null(v)) {
+    vertex.frame.color <- vertex.frame.color[v]
+  }
+  vertex.size <- 1/200 * params("vertex", "size")
+  if (length(vertex.size) != 1 && !is.null(v)) {
+    vertex.size <- vertex.size[v]
+  }
+  
+  symbols(x=coords[,1], y=coords[,2], bg=vertex.color,
+          stars=cbind(vertex.size, vertex.size, vertex.size, vertex.size),
+          add=TRUE, inches=FALSE)
+}
+add_shape("diamond", clip=shapes("circle")$clip,
+          plot=MyTriangle, parameters=list(vertex.frame.color="white",
+                                           vertex.frame.width=1))
\ No newline at end of file
diff --git a/PosteriorProb.R b/PosteriorProb.R
new file mode 100644
index 0000000..2a8d8d2
--- /dev/null
+++ b/PosteriorProb.R
@@ -0,0 +1,101 @@
+#### Ecological Synthesis Lab (SintECO)
+
+#### Authors: Gabriel M. Felix, Rafael B. P. Pinheiro, and Marco A. R. Mello.
+
+#### See README for further info.
+
+#### Compute probabilities of interaction in a network with a modular structure.
+
+
+PosteriorProb <- function(M, R.partitions, C.partitions, Prior.Pij, Conditional.level){
+  
+  # Test of assumptions
+  if (!is.matrix(M)){stop("M is not a matrix")}
+  if (0 %in% rowSums(M) | 0 %in% colSums(M)) {stop("M is degenerated. There are rows and/or columns without interactions in the matrix. Remove them before proceding")}
+  if (!is.numeric(R.partitions) | !is.numeric(C.partitions)) {stop("Partitions are not numeric")}
+  if (length(R.partitions) != nrow(M) | length(C.partitions) != ncol(M)) {stop("Partitions and matrix dimensions have different sizes")}
+  if (!(Conditional.level %in% c("matrix","modules","areas"))) {stop("Conditional.level should be 'matrix','modules' or 'areas'")}
+  if (Prior.Pij != "degreeprob" & Prior.Pij != "equiprobable" & Prior.Pij != "degreeprob.byarea") {stop("Pij.probs should be 'equiprobable' or 'degreeprob' or 'degreeprob.byarea")}
+  
+  # M dimensions
+  r <- dim(M)[1] # Number of rows
+  c <- dim(M)[2] # Number of columns
+  array()
+  
+  # Making an array with r rows, c columns, and 3 slices. This array represents the modular structure.
+  # The first slice informs if a given cell M(rc) is within (1) or outside (0) a module.
+  # The second slice informs to which module the species in the row (r) of a given cell M(rc) belongs.
+  # The third slice informs to which module the species in the column (c) of a given cell M(rc) belongs .
+  
+  Matrix.mod <- array(0, dim = c(r, c, 3))
+  for (rr in 1:r){
+    for (cc in 1:c){
+      Matrix.mod[rr,cc,1] <- ifelse(R.partitions[rr] == C.partitions[cc], 1,0)
+      Matrix.mod[rr,cc,2] <- R.partitions[rr]
+      Matrix.mod[rr,cc,3] <- C.partitions[cc]
+    }
+  }
+  
+  # Defining a priori Pij probabilities.
+  if (Prior.Pij == "equiprobable"){
+    Pi <- rep(1 / r, times = r) 
+    Pj <- rep(1 / c, times = c)
+    Prior.Pij.species <- tcrossprod(Pi, Pj)
+  }else if (Prior.Pij == "degreeprob"){
+    Pi <- rowSums(M) / sum(rowSums(M)) 
+    Pj <- colSums(M) / sum(colSums(M))
+    Prior.Pij.species <- tcrossprod(Pi, Pj)
+  }else if(Prior.Pij == "degreeprob.byarea"){
+    Prior.Pij.species <- M
+    RMod <- sort(unique(R.partitions))
+    CMod <- sort(unique(C.partitions))
+    for (rr in RMod){
+      for (cc in CMod){
+        M.rr.cc <- matrix(M[R.partitions == rr,C.partitions == cc], sum(1*(R.partitions == rr)), sum(1*(C.partitions == cc)))
+        Pi.rr.cc <- rowSums(M.rr.cc) / sum(rowSums(M.rr.cc)) 
+        Pj.rr.cc <- colSums(M.rr.cc) / sum(colSums(M.rr.cc))
+        Prior.Pij.species[R.partitions == rr, C.partitions == cc] <- tcrossprod(Pi.rr.cc, Pj.rr.cc)
+      }
+    }
+  }
+
+  # Defining conditional probabilities by area based on species degrees and connectance by area.
+  
+  if (Conditional.level == "matrix"){
+    Post.Pij <- Prior.Pij.species
+  }else {
+    Prior.Pij.area <- matrix(NA,r,c)
+    Cond.Pij.area <- matrix(NA,r,c)
+    if (Conditional.level == "modules"){
+      WMod.prior <- sum(Prior.Pij.species[Matrix.mod[,,1] == 1])
+      OMod.prior <- sum(Prior.Pij.species[Matrix.mod[,,1] == 0])
+      Prior.Pij.area[Matrix.mod[,,1] == 1] <- WMod.prior
+      Prior.Pij.area[Matrix.mod[,,1] == 0] <- OMod.prior
+    
+      WMod.cond <- sum(M[Matrix.mod[,,1] == 1]) / sum(M)
+      OMod.cond <- sum(M[Matrix.mod[,,1] == 0]) / sum(M)
+      Cond.Pij.area[Matrix.mod[,,1] == 1] <- WMod.cond
+      Cond.Pij.area[Matrix.mod[,,1] == 0] <- OMod.cond
+    }else if (Conditional.level == "areas"){
+      RMod <- sort(unique(R.partitions))
+      CMod <- sort(unique(C.partitions))
+      for (rr in RMod){
+        for (cc in CMod){
+          WArea.prior <- sum(Prior.Pij.species[Matrix.mod[,,2] == rr & Matrix.mod[,,3] == cc])
+          Prior.Pij.area[Matrix.mod[,,2] == rr & Matrix.mod[,,3] == cc] <- WArea.prior
+        
+          WArea.cond <- sum(M[Matrix.mod[,,2] == rr & Matrix.mod[,,3] == cc]) / sum(M)
+          Cond.Pij.area[Matrix.mod[,,2] == rr & Matrix.mod[,,3] == cc] <- WArea.cond
+        }
+      } 
+    }
+  
+    # Adjusting the prior Pij prob by conditional probabilities. 
+  
+    Post.Pij <- Prior.Pij.species * (Cond.Pij.area / Prior.Pij.area)
+  }
+  
+  return(Post.Pij = Post.Pij)
+}
+
+
diff --git a/README.md b/README.md
index 9b3baf0..35ad7ca 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,114 @@
 # queiroz_et_al_2020
-Supplements to the paper Queiroz et al. (2020, Biotropica)
+
+Supplement to the paper Queiroz et al. (2020, Biotropica).
+
+[Ecological Synthesis Lab](https://marcomellolab.wordpress.com) (SintECO).
+
+Authors: Joel A. Queiroz, Ugo M. Diniz, Diego P. Vázquez, Zelma M. Quirino, Francisco A.R. Santos, Marco A.R. Mello, Isabel C. Machado.
+
+E-mail: imachado@ufpe.br. 
+
+Published on October 2nd, 2020 (English version).
+
+Run in R version 4.0.2 (2020-06-22) -- "Taking Off Again".
+
+Disclaimer: You may freely use the software and data provided here for commercial or non-commercial purposes at your own risk. We assume no responsibility or liability for the use of this material, convey no license or title under any patent, copyright, or mask work right to the product. We reserve the right to make changes in the material without notification. We also make no representation or warranty that such application will be suitable for the specified use without further testing or modification. If this material helps you produce any academic work (paper, book, chapter, monograph, dissertation, report or similar), please acknowledge the authors and cite the source.
+
+
+## Functionality and origin
+
+The data, scripts, and functions provided here aim at making or study fully reproducible. You will find code to reproduce both the analysis and the figures, as well as the main supplementary material. We've also added some bonus material related to data visualization.
+
+
+## List of folders and files
+
+1. data (folder)
+
+    a. morph_pla.xlsx -> morphometric data of plants.
+    
+    b. morph_pol.xlsx - > morphometric data of pollinators.
+    
+    c. network.txt -> the nocturnal pollination network composed of plants, bats, and hawkmoths.
+    
+    d. plants.xlsx -> guild data of plants.
+    
+    e. sampbat.xlsx -> capture data of bats.
+    
+    f. samphawk.xlsx -> capture data of hawkmoths.
+    
+    g. sampling_bats.xlsx -> sampling data of bats at the individual level.
+    
+    h. sampling_hawkmoths.xlsx -> sampling data of hawkmoths at the individual level.
+    
+
+2. figures (folder)
+
+    a. bc.tiff -> boxplot of betweenness centrality by pollination syndrome.
+    
+    b. compound.png -> matrix of the nocturnal pollination network.
+
+    c. d.tiff -> boxplot of specialization (d') by pollination syndrome.
+    
+    d. morph.tiff -> boxplot comparing morphometric data between network modules.
+
+    e. network.tif -> graph of the nocturnal pollination network.
+
+    f. nk.tiff -> boxplot of normalized degree by pollination syndrome.
+
+    g. sampling.tiff -> rarefaction curves representing interaction sampling completeness.
+
+
+3. results (folder)
+
+    a. dplants.csv -> scores of specialization (d') for the plants.
+
+    b. BCplants.csv -> scores of betweenness centrality for the plants.
+
+    c. NDplants.csv -> scores of normalized degree for the plants.
+
+    d. results_centrality.txt -> results of the GLMs comparing centraliy by pollination syndrome.
+    
+    e. results_topology.txt -> results of the topological analyses.
+
+
+4. analysis.R -> main script to run the analyses and plot the graphs.
+
+5. MyTriangle.R -> additional function to create the diamond shape used in "network.tif".
+
+6. PosteriorProb.R -> additional function to calculate restricted probabilities to be used in the compound topology analysis. See the original [repo](https://github.com/gabrielmfelix/Restricted-Null-Model).
+
+7. RestNullModel.R -> additional function to run the compound topology analysis. See the original [repo](https://github.com/gabrielmfelix/Restricted-Null-Model).
+
+
+## Instructions
+
+1. Run the script "analysis.R" and follow the steps provided in it;
+
+2. You can also run only some sections of the script to reproduce specific analyses, but beware of dependences between sections.
+
+3. The analyses based on Monte Carlo simulations can be quite memory-consuming. Check your computer's power before setting the number of permutations. If you want to analyze some data for real using our script, consider that you'll need at least 1,000 permutations. If you want only to explore some data, 10 will do it.
+
+4. If have any questions, corrections, or suggestions, please feel free to open an *issue* or make a *pull request*.
+
+5. Have fun!
+
+
+## Original sources
+
+The codes used in this repo have borrowed solutions from other repos developed by or lab:
+
+1. [network-significance](https://github.com/marmello77/network-significance)
+
+2. [Restricted-Null-Model](https://github.com/gabrielmfelix/Restricted-Null-Model)
+
+3. [ihsmodel](https://github.com/pinheirorbp/ihsmodel)
+
+
+## Acknowledgments
+
+We thank our colleagues, who helped us at different stages of this project. We also thank our sponsors, who funded this study through grants, fellowships, and scholarships. Last, but not least, we thank the [Stack Overflow Community](https://stackoverflow.com), where we solve most of our coding dilemmas.
+
+
+## Reference
+
+Queiroz JA, Diniz UM, Vázquez DP, Quirino ZM, Santos FAR, Mello MAR, Machado IC. 2020. Bats and hawkmoths form mixed modules with flowering plants in a nocturnal interaction network. Biotropica, *accepted*.
diff --git a/RestNullModel.R b/RestNullModel.R
new file mode 100644
index 0000000..3fb9d43
--- /dev/null
+++ b/RestNullModel.R
@@ -0,0 +1,150 @@
+#### Ecological Synthesis Lab (SintECO)
+
+#### Authors: Gabriel M. Felix, Rafael B. P. Pinheiro, and Marco A. R. Mello.
+
+#### See README for further info.
+
+#### Vaznull algorithm of bipartite modified to run the restricted null model
+
+
+RestNullModel <- function(M, Pij.Prob, Numbernulls, Print.null = F, allow.degeneration = F, 
+                              return.nonrm.species = T, connectance = T, byarea = F, R.partitions = NULL, C.partitions = NULL){
+  
+  ### Test of assumptions
+  
+  if (!is.matrix(M)){stop("M is not a matrix")}
+  if (0 %in% rowSums(M) | 0 %in% colSums(M)) {stop("M is degenerated")}
+  
+  if (!is.matrix(Pij.Prob)){stop("Pij is not a matrix")}
+  if (T %in% c(Pij.Prob < 0)){stop("Pij must contain only numbers >= 0")}
+  
+  if (nrow(M) != nrow(Pij.Prob) | ncol(M) != ncol(Pij.Prob)){stop("Dimensions of M and Pij.Prob should be identical")}
+  
+  if (byarea == T){
+    if(is.null(C.partitions) | is.null(R.partitions)){stop("Partitions missing")}
+    if (length(unique(c(length(R.partitions),nrow(M),nrow(Pij.Prob)))) != 1){stop("The number of elements of R.partition should be the same as the number of rows of M and Pij.prob")}
+    if (length(unique(c(length(C.partitions),ncol(M),ncol(Pij.Prob)))) != 1){stop("The number of elements of C.partition should be the same as the number of column of M and Pij.prob")}
+     if(!identical(sort(unique(R.partitions)), sort(unique(C.partitions)))){stop("The number and labels of modules in R.partition and C.partition must be the same")}
+  }  
+  
+  if (Numbernulls <= 0 | !is.numeric(Numbernulls)) {stop("Numbernulls should be a number > 0")}
+  if (!is.logical(connectance)){stop("connectance should be logical (T or F)")}
+  if (!is.logical(allow.degeneration)){stop("allow.degeneration should be logical (T or F)")}
+  if (!is.logical(return.nonrm.species)){stop("return.nonrm.species should be logical (T or F)")}
+  if (!is.logical(byarea)){stop("byarea should be logical (T or F)")}
+  
+  ### M dimensions
+  r <- dim(M)[1] # Number of rows
+  c <- dim(M)[2] # Number of collums
+  
+  ### Constructing a array with r rows, c columns and 2 slices. This array represents the matrix area structure
+  if (byarea == T){
+    Matrix.area <- array(0, dim = c(r, c, 2))
+    for (rr in 1:r){
+      for (cc in 1:c){
+        Matrix.area[rr,cc,1] <- R.partitions[rr]
+        Matrix.area[rr,cc,2] <- C.partitions[cc]
+      }
+    }
+    
+  }else if (byarea == F){
+    ## Assigning all rows and columns to the same partition in order to run the code bellow
+    Matrix.area <- array(1, dim = c(r, c, 2))
+    R.partitions <- rep(1, nrow(M))
+    C.partitions <- rep(1, ncol(M))
+    
+  }
+  
+  ### Null model simulation
+  
+  NullMatrices <- list() # list where the null matrices will be storage
+  length(NullMatrices) <- Numbernulls #assigning the number of null matrices to be saved in NullMatrices 
+  
+  ## Drawing interaction in each null matrix 
+  for (nn in 1:Numbernulls){
+      
+    R.part <- sort(unique(as.vector(Matrix.area[,,1])))
+    C.part <- sort(unique(as.vector(Matrix.area[,,2])))
+    finalmat <- matrix(NA, r, c)
+    
+    for (R.p in R.part){
+      for (C.p in C.part){
+          
+        M.a <- as.matrix(M[R.partitions == R.p, C.partitions == C.p])
+        Pij.a <- Pij.Prob[R.partitions == R.p, C.partitions == C.p]
+          
+        r.a <- dim(M.a)[1]
+        c.a <- dim(M.a)[2]
+        
+        P.a <- P1.a <- Pij.a
+        finalmat.a <- matrix(0, r.a, c.a)
+        
+        if(allow.degeneration == F & R.p == C.p){
+          
+          ## Ensuring that the dimensions of the null matrix will be the same of the original matrix
+          
+          D.int.finalmat.a <- 0 # The number of rows + columns occupied of the null matrix 
+          while (D.int.finalmat.a < sum(dim(M.a))) { # While the dimensions of the null matrix was smaller then the original matrix, keep going
+            sel <- sample(1:length(M.a), 1, prob = P.a) # Sample an cell of M.a with probability P.a
+            selc <- floor((sel - 1)/(dim(M.a)[1])) + 1 # Recovering column and 
+            selr <- ((sel - 1)%%dim(M.a)[1]) + 1 # row of the cell sampled
+            if (sum(finalmat.a[, selc]) == 0 | sum(finalmat.a[selr,]) == 0) { # Checking if row or column of the sampled cell is empty
+              finalmat.a[sel] <- 1 
+              P.a[sel] <- 0
+            }
+          D.int.finalmat.a <- sum(rowSums(finalmat.a) > 0) + sum(colSums(finalmat.a) > 0) # Setting the new number of dimensions occupied
+          }
+          # When the number of occupied dimensions of the null matrix was the same as the original matrix, continue
+        }
+        
+        conn.remain <- sum(M.a > 0) - sum(finalmat.a > 0) # The number of cells remaining to be occupied to mantain the original connectance
+          
+          if (conn.remain > 0) {
+              if(connectance == T){
+                if (length(which(finalmat.a == 0)) == 1) {
+                  add <- which(finalmat.a == 0)
+                } else {
+                add <- sample(which(finalmat.a == 0), conn.remain, 
+                              prob = P1.a[finalmat.a == 0], replace = F)
+                }
+              }else {
+                add <- sample(1:length(finalmat.a), conn.remain, 
+                              prob = P1.a, replace = T)
+              }
+            for (add1 in add){
+              finalmat.a[add1] <- finalmat.a[add1] + 1
+            }
+          }
+          
+          ### Checking if there are still interactions to be drawn. If applicable, draw.
+          int.remain <- (sum(M.a) - sum(finalmat.a))
+          if (int.remain > 0) {
+            if (length(which(finalmat.a > 0)) == 1) {
+              add <- rep(which(finalmat.a > 0),int.remain)
+            }else{
+              add <- sample(which(finalmat.a > 0), int.remain, prob = P1.a[which(finalmat.a >0)], replace = T)
+            }
+            finalmat.a[as.numeric(names(table(add)))] <- finalmat.a[as.numeric(names(table(add)))] + (table(add)) 
+          }
+      
+        finalmat[R.partitions == R.p, C.partitions == C.p] <- finalmat.a
+      }
+    }
+    
+    # Saving outputs
+    R2keep <- which(rowSums(finalmat) != 0)
+    C2keep <- which(colSums(finalmat) != 0)
+    finalmat2 <- finalmat[R2keep,C2keep]
+    if (return.nonrm.species == T){
+      NullMatrices[[nn]] = list(NullMatrix = finalmat2, R.Kept = R2keep, C.Kept = C2keep)
+    }else if(return.nonrm.species == F){
+      NullMatrices[[nn]] = finalmat2
+    }
+    if (Print.null == T){print(nn)}
+  }
+  return(NullMatrices = NullMatrices)
+}
+
+
+
+
diff --git a/analysis.R b/analysis.R
new file mode 100644
index 0000000..f7dcb1d
--- /dev/null
+++ b/analysis.R
@@ -0,0 +1,592 @@
+################################################################################
+################################################################################
+#Supplement to the paper Queiroz et al. (2020, Biotropica).
+
+#Ecological Synthesis Lab (SintECO): https://marcomellolab.wordpress.com.
+
+#Authors: Joel A. Queiroz, Ugo M. Diniz, Diego P. Vázquez, Zelma M. Quirino,
+#Francisco A.R. Santos, Marco A.R. Mello, Isabel C. Machado.
+
+#See README for further info:
+#https://github.com/marmello77/queiroz_et_al_2020/blob/main/README.md
+################################################################################
+################################################################################
+
+
+
+################################################################################
+##### SET THE STAGE
+################################################################################
+
+
+#Set the working directory
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+
+#Delete all previous objects
+rm(list= ls())
+
+#Load the required packages
+library(igraph)
+library(bipartite)
+library(Rmisc)
+library(vegan)
+library(gdata)
+library(ggplot2)
+
+#Load some custom-made functions
+source("RestNullModel.R")
+source("PosteriorProb.R")
+source("MyTriangle.R")
+
+
+
+################################################################################
+##### PROCESS THE NETWORK
+################################################################################
+
+
+#Import the network
+data <- as.matrix(read.delim("data/network.txt", row.names=1))
+
+#Inspect the network
+class(data)
+data
+dim(data)
+min(data)
+max(data)
+
+#Plot the matrix
+visweb(data)
+
+#Convert the network to igraph format
+data2 <- graph_from_incidence_matrix(data, directed = F, weighted = TRUE) 
+
+#Inspect object
+class(data2)
+data2
+E(data2)
+V(data2)$name
+
+#Inform which nodes represent which taxonomic groups
+V(data2)$set[1:nrow(data)] = c("Moths",	"Moths", "Bats",	"Bats",	"Moths",	"Moths",
+                       "Moths",	"Moths",	"Moths",	"Bats",	"Bats",	"Moths",
+                       "Moths",	"Moths",	"Moths",	"Moths",	"Moths",
+                       "Moths",	"Moths")
+V(data2)$set[(nrow(data)+1):(nrow(data)+ncol(data))] = "Plants"
+
+
+
+################################################################################
+##### DRAW THE NETWORK
+################################################################################
+
+
+#Set layout
+lay1 <- layout_nicely(data2)
+
+#Set edge curvatures
+curves1 = curve_multiple(data2)
+
+#Set edge mode and width
+E(data2)$arrow.mode = 0
+E(data2)$width = E(data2)$weight/5+1
+
+#Calculate Louvain modularity (resolution = 1.0, similar to DIRT_LPA+)
+data2.lou = cluster_louvain(data2)
+
+#Import "diamond" vertex shape
+source("MyTriangle.R")
+
+#Set vertex shapes
+V(data2)$shape = V(data2)$set
+V(data2)$shape = gsub("Bats","diamond",V(data2)$shape)
+V(data2)$shape = gsub("Moths","square",V(data2)$shape)
+V(data2)$shape = gsub("Plants","circle",V(data2)$shape)
+
+##Set node and cloud colors by modularity
+colrs <- rainbow(length(data2.lou), alpha = 1.0, s = 1, v = 0.8)
+V(data2)$color <- colrs[data2.lou$membership]
+clouds = rainbow(length(data2.lou), alpha = 0.1)
+
+#Plot and export the graph
+tiff(filename= "figures/network.tif", res= 300, height= 3000, width= 3100)
+par(mfrow=c(1,1),mar=c(1,1,1,5))
+plot(data2.lou,
+     data2,
+     col = V(data2)$color,
+     mark.border="lightgrey", 
+     mark.col=clouds,
+     vertex.size=7.5,
+     vertex.label=V(data2)$name,
+     vertex.label.color="white",
+     vertex.label.cex=.3,
+     edge.color = adjustcolor("grey", alpha.f = .5),
+     edge.curved=0.3,
+     edge.width = 3,
+     layout=lay1)
+legend(x = 0.9,y = 1.0, legend = c("Bats", "Moths", "Plants"),
+       pch = c(18,15,19),  title="Taxon",
+       text.col = "gray20", title.col = "black",
+       box.lwd = 0, cex = 2, col=c("grey", "grey", "grey"))
+par(mfrow=c(1,1))
+dev.off()
+
+
+
+################################################################################
+##### NETWORK LEVEL ANALYSIS
+################################################################################
+
+
+#Set seed
+set.seed(14)
+
+#Set the number of permutations to be used in all null model analyses
+permutations <- 1000
+
+#Generate randomized matrices
+nulls <- nullmodel(data, N=permutations, method="vaznull")
+
+
+##### MODULARITY
+
+#Calculate metric for the original network
+Mod <- computeModules(data, method = "Beckett")
+Mod@likelihood
+
+#Extract module membership
+Part <- bipartite::module2constraints(Mod)
+row.Part <- Part[1:nrow(data)]
+col.Part <- Part[(nrow(data)+1):(nrow(data)+ncol(data))]
+
+#Calculate metric for the randomized networks
+nullmod <- sapply(nulls, computeModules, method = "Beckett")
+modnull <- sapply(nullmod, function(x) x@likelihood)
+(Mod@likelihood - mean(modnull))/sd(modnull) # Z value
+Mod.sig <- sum(modnull>(Mod@likelihood)) / length(modnull) # p value
+Mod.sig
+
+#Plot the observed value against the distribution of randomized values
+png("specialization.png", width = 3000, height = 3000, res = 300)
+plot(density(modnull), main="Observed vs. randomized",
+     xlim=c(min((Mod@likelihood), min(modnull)), 
+            max((Mod@likelihood), max(modnull))))
+abline(v=Mod@likelihood, col="red", lwd=2, xlab="")
+dev.off()
+
+
+Mod@likelihood #observed
+mean(modnull) #randomized mean
+sd(modnull) #randomized SD
+(Mod@likelihood - mean(modnull))/sd(modnull) # Z-value
+sum(modnull>(Mod@likelihood)) / length(modnull) #randomized > observed
+sum(modnull<(Mod@likelihood)) / length(modnull) #randomized < observed
+
+
+##### SPECIALIZATION 
+
+#Calculate metric for the original network
+Spec <- networklevel(data, index="H2")
+class(Spec)
+
+#Calculate metric for the randomized networks
+randomized.Spec <- unlist(sapply(nulls, networklevel, index="H2"))
+(Spec - mean(randomized.Spec))/sd(randomized.Spec) # Z value
+Spec.sig <- sum(randomized.Spec>Spec)/length(randomized.Spec) # p value
+Spec.sig
+
+#Plot the observed value against the distribution of randomized values
+png("specialization.png", width = 3000, height = 3000, res = 300)
+plot(density(randomized.Spec), main="Observed vs. randomized",
+     xlim=c(min((Spec), min(randomized.Spec)), 
+            max((Spec), max(randomized.Spec))))
+abline(v=Spec, col="red", lwd=2, xlab="")
+dev.off()
+
+Spec #observed
+mean(randomized.Spec) #randomized mean
+sd(randomized.Spec) #randomized SD
+(Spec - mean(randomized.Spec))/sd(randomized.Spec) # Z-value
+sum(randomized.Spec>(Spec)) / length(randomized.Spec) #randomized > observed
+sum(randomized.Spec<(Spec)) / length(randomized.Spec) #randomized < observed
+
+
+##### NESTEDNESS
+
+#Calculate metric for the original network
+Nest <- networklevel(data, index="weighted NODF")
+
+#Calculate metric for the randomized networks
+randomized.Nest <- unlist(sapply(nulls, networklevel, index="weighted NODF"))
+(Nest - mean(randomized.Nest))/sd(randomized.Nest) # Z value
+Nest.sig <- sum(randomized.Nest>Nest)/length(randomized.Nest) # p value
+Nest.sig
+
+#Plot the observed value against the distribution of randomized values
+png("nestedness.png", width = 3000, height = 3000, res = 300)
+plot(density(randomized.Nest), main="Observed vs. randomized",
+     xlim=c(min((Nest), min(randomized.Nest)), 
+            max((Nest), max(randomized.Nest))))
+abline(v=Nest, col="red", lwd=2, xlab="")
+dev.off()
+
+Nest #observed
+mean(randomized.Nest) #randomized mean
+sd(randomized.Nest) #randomized SD
+(Nest - mean(randomized.Nest))/sd(randomized.Nest) # Z-value
+sum(randomized.Nest>(Nest)) / length(randomized.Nest) #randomized > observed
+sum(randomized.Nest<(Nest)) / length(randomized.Nest) #randomized < observed
+
+
+##### COMPOUND TOPOLOGY 
+
+#Calculate the desired nestedness metric (here WNODA) for the original network.
+obs.com <- unlist(bipartite::nest.smdm(x = data, 
+                                       constraints = Part, #Input the modular structured recovered from step 2
+                                       weighted = T, #By considering the edge weights, you are choosing WNODA
+                                       decreasing = "abund"))
+
+#Check the scores
+obs.com
+
+#Calculate constrained interaction probabilities considering the network's modular structure
+Pij <- PosteriorProb(M = data, 
+                     R.partitions = row.Part, C.partitions = col.Part, #Input the modular structured recovered from step 2
+                     Prior.Pij = "degreeprob", #Choose the null model
+                     Conditional.level = "modules") #Choose the kind of constraints
+
+#Check what those probabilities look like
+Pij
+dim(Pij)
+
+#Generate randomized networks with the null model of your choice, considering the interaction probabilities calculated before. 
+nulls.com <- RestNullModel(M = data, 
+                           Pij.Prob = Pij, #Recover the probabilities calculated in the previous command
+                           Numbernulls = permutations, #This step may take long, so start experimenting with low values
+                           Print.null = T, 
+                           allow.degeneration = F, #Choose whether you allow orphan rows and columns to be removed or not
+                           return.nonrm.species = F, 
+                           connectance = T, byarea = T, 
+                           R.partitions = row.Part, 
+                           C.partitions = col.Part)
+
+#Calculate the nestedness within and between modules
+rest.nest <- nest.smdm(data, constraints = Part, 
+                       weighted = T, 
+                       decreasing = "abund", 
+                       sort = T)
+
+unlist(rest.nest)
+
+null.com <- sapply(nulls.com, 
+                   function(x) bipartite::nest.smdm(x = x,
+                                                    constraints = Part, 
+                                                    weighted = T, 
+                                                    decreasing = "abund"))
+WNODA.null.com <- unlist(null.com[3,])
+WNODAsm.null.com <- unlist(null.com[8,])
+WNODAdm.null.com <- unlist(null.com[9,])
+
+#Plot the observed nestedness value against the distribution of randomized values
+par(mfrow = c(1,3))
+
+plot(density(WNODA.null.com), xlim=c(min(obs.com[3], min(WNODA.null.com)),
+                                     max(obs.com[3], max(WNODA.null.com))), 
+     main="observed vs. randomized", xlab = "WNODA matrix")
+abline(v=obs.com[3], col="red", lwd=2)
+
+plot(density(WNODAsm.null.com), xlim=c(min(obs.com[8], min(WNODAsm.null.com)),
+                                       max(obs.com[8], max(WNODAsm.null.com))), 
+     main="observed vs. randomized", xlab = "WNODAsm matrix")
+abline(v=obs.com[8], col="red", lwd=2)    
+
+plot(density(WNODAdm.null.com), xlim=c(min(obs.com[9], min(WNODAdm.null.com)),
+                                       max(obs.com[9], max(WNODAdm.null.com))), 
+     main="observed vs. randomized", xlab = "WNODAdm matrix")
+abline(v=obs.com[9], col="red", lwd=2)    
+
+par(mfrow = c(1,1))
+
+#Estimate the p-values
+
+#Nestedness in th entire network
+praw.WNODA <- sum(WNODA.null.com>obs.com[3]) / length(WNODA.null.com)
+p.WNODA <- ifelse(praw.WNODA > 0.5, 1- praw.WNODA, praw.WNODA)    # P-value
+p.WNODA
+
+#Nestedness within the modules
+praw.WNODAsm <- sum(WNODAsm.null.com>obs.com[8]) / length(WNODAsm.null.com)
+p.WNODAsm <- ifelse(praw.WNODAsm > 0.5, 1- praw.WNODAsm, praw.WNODAsm)    # P-value
+p.WNODAsm
+
+#Nestedness between the modules
+praw.WNODAdm <- sum(WNODAdm.null.com>obs.com[9]) / length(WNODAdm.null.com)
+p.WNODAdm <- ifelse(praw.WNODAdm > 0.5, 1- praw.WNODAdm, praw.WNODAdm)    # P-value
+p.WNODAdm
+
+
+##### PLOT THE COMPOUND TOPOLOGY
+
+par(mfrow = c(1,1))
+
+#Sort the matrix in a way that facilitates visualizing the compound topology
+data.comp <- bipartite::sortmatrix(matrix = data, topology = "compound", 
+                                   sort_by = "weights", 
+                                   row_partitions = row.Part, 
+                                   col_partitions = col.Part)
+
+#Assign colors for the modules
+modcol <- rainbow((length(unique(Part))), alpha=1, s = 1, v = 1)
+
+#Plot the matrix
+png("figures/compound.png", width = 3000, height = 3000, res = 300)
+plotmatrix(data.comp$matrix, 
+           row_partitions = data.comp$row_partitions, 
+           col_partitions = data.comp$col_partitions, 
+           border = T,
+           binary = F,
+           modules_colors = modcol,
+           within_color = modcol, 
+           between_color = "lightgrey")
+dev.off()
+
+
+##### EXPORT A SUMMARY OF THE TOPOLOGICAL RESULTS
+
+sink(file = "results/results_topology.txt")
+
+paste("Topological analysis of the nocturnal pollination network")
+paste("Queiroz et al. 2020, Biotropica")
+cat("\n")
+paste("The network has", nrow(data), "rows and", ncol(data), "columns.")
+cat("\n")
+paste("The network's specialization (H2) is", round(Spec, 2),",", "P =", round(Spec.sig, 2))
+cat("\n")
+paste("The network's modularity (DIRT_LPA+) is", round(Mod@likelihood, 2), ",", "P =", round(Mod.sig, 2), ",", "and it contains", length(unique(Part)),"modules.")
+cat("\n")
+paste("The network's nestedness (WNODF) is", round(Nest/100, 2),",", "P =", round(Nest.sig, 2))
+cat("\n")
+paste("The network shows the following scores of nestedness (WNODA):")
+cat("\n")
+paste("Entire network =", round(rest.nest$WNODAmatrix/100, 2), ",", "P =", round(p.WNODA, 2))
+cat("\n")
+paste("Between the modules =", round(rest.nest$WNODA_DM_matrix/100, 2), ",", "P =", round(p.WNODAdm, 2))
+cat("\n")
+paste("Within the modules =", round(rest.nest$WNODA_SM_matrix/100, 2), ",", "P =", round(p.WNODAsm, 2))
+cat("\n")
+
+sink(file = NULL, )
+
+
+
+################################################################################
+##### SPECIES LEVEL ANALYSIS (CENTRALITY)
+################################################################################
+
+
+# Specialization (d') 
+d <- specieslevel(data,index="d")
+dplants <- d$`higher level`
+write.csv(dplants, "results/dplants.csv") # writing a separate csv file   
+
+# Betweenness centrality (BC)
+BC <- specieslevel(data, index="betweenness")
+BCplants <- BC$higher
+write.csv(BCplants, "results/BCplants.csv") # writing a separate csv file
+
+# Normalized degree (nk)
+ND <-ND(data, normalised=T)
+NDplants <- ND$higher
+write.csv(NDplants, "results/NDplants.csv") # writing a separate csv file
+
+
+##### Plot centrality metrics by syndrome
+
+#Import data
+plants <- read.xls("data/plants.xlsx", h=T) # reading compiled spreadsheet with species & metrics classified by guild
+
+# Change reference level for GLMs
+ord <-  ordered(plants$Guild, levels = c("sphin", "chiro", "other"))
+
+# BC
+
+tiff("figures/bc.tiff", width = 12, height = 15, units = "cm", res = 300)
+ggplot(plants, aes(x=ord, y=bc, fill=Guild)) + 
+  ylab("BC")+ xlab("")+ ylim(0, 0.15) +
+  scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+  geom_boxplot(width=0.5, color="black") +
+  theme_classic() +
+  theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+        axis.title.y = element_text(color="black", face ="italic", size =23),
+        axis.text= element_text(color="black", size=19),
+        legend.position = "none")
+dev.off()
+
+# d'
+
+tiff("figures/d.tiff", width = 12, height = 15, units = "cm", res = 300)
+
+ggplot(plants, aes(x=ord, y=d, fill=Guild)) + 
+  ylab("d'")+ xlab("")+ ylim(0, 0.8) +
+  scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+  geom_boxplot(width=0.5, color="black") +
+  theme_classic() +
+  theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+        axis.title.y = element_text(color="black", face ="italic", size =23),
+        axis.text= element_text(color="black", size=19),
+        legend.position = "none")
+
+dev.off()
+
+# ND
+
+tiff("figures/nk.tiff", width = 12, height = 15, units = "cm", res = 300)
+
+ggplot(plants, aes(x=ord, y=nk, fill=Guild)) + 
+  ylab("nk")+ xlab("")+ ylim(0, 1.1) +
+  scale_fill_manual(values=c("darkolivegreen1", "sandybrown", "orchid1"))+
+  geom_boxplot(width=0.5, color="black") +
+  theme_classic() +
+  theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+        axis.title.y = element_text(color="black", face ="italic", size =23),
+        axis.text= element_text(color="black", size=19),
+        legend.position = "none")
+
+dev.off()
+
+
+##### Run GLMs to compare centrality scores
+
+table(plants$Guild)
+plants$Guild <- relevel(plants$Guild, ref="chiro") #changing reference level for GLMs
+plants$Guild <- relevel(plants$Guild, ref="sphin")
+plants$Guild <- relevel(plants$Guild, ref="other")
+
+# d'
+
+glmd <- glm(plants$d ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmd)
+glm_d <- anova(glmd, test = "Chisq")
+glm_d
+
+## BC
+
+glmbc <- glm(plants$bc ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmbc)
+glm_bc <- anova(glmbc, test = "Chisq")
+glm_bc
+
+## nk
+
+glmnk <- glm(plants$nk ~ plants$Guild, family=quasibinomial("logit"))
+summary(glmnk)
+glm_k <- anova(glmnk, test = "Chisq")
+glm_k
+
+#Export the results
+sink(file = "results/results_centrality.txt")
+
+paste("Comparison of centrality by syndrome")
+paste("Queiroz et al. 2020, Biotropica")
+cat("\n")
+paste(capture.output(glmd))
+paste(capture.output(glm_d))
+cat("\n")
+paste(capture.output(glmbc))
+paste(capture.output(glm_bc))
+cat("\n")
+paste(capture.output(glmnk))
+paste(capture.output(glm_k))
+
+sink(file = NULL, )
+
+
+
+################################################################################
+##### MORPHOMETRIC ANALYSIS
+################################################################################
+
+
+##### Import morphology data
+morph_plants<-read.xls("data/morph_pla.xlsx", h=T)
+morph_pol <- read.xls("data/morph_pol.xlsx", h=T)
+
+str(morph_plants)
+str(morph_pol)
+
+# Change reference level for GLMs
+morph_plants$module <- relevel(morph_plants$module, ref="bat")
+morph_pol$module <- relevel(morph_pol$module, ref="hawk1")
+
+
+##### Run GLMs to compare modules
+
+# Pollinator tongues
+glm_pol <- glm(morph_pol$length_pol~morph_pol$module, family=gaussian())
+summary(glm_pol)
+anova(glm_pol, test = "Chisq")
+
+# Floral width (w) and length (l)
+glm_pla_l <- glm(morph_plants$length_pla~morph_plants$module, family=gaussian())
+glm_pla_w <- glm(morph_plants$width_pla~morph_plants$module, family=gaussian())
+summary(glm_pla_l)
+anova(glm_pla_l, test = "Chisq")
+summary(glm_pla_w)
+anova(glm_pla_w, test = "Chisq")
+
+### Plot
+ggmorph<-read.xls("data/morph_graph.xlsx",h=T) 
+
+tiff("figures/morph.tiff", width = 24, height = 10, units = "cm", res = 300)
+ggplot(ggmorph, aes(x=module, y=measure, fill=variable)) +
+  ylab("Measure (mm)")+ xlab("Modules")+ ylim(0, 150) +
+  scale_fill_manual(values=c("slateblue", "goldenrod1", "coral1"))+
+  geom_boxplot(width=0.5, color="black", position = position_dodge(width=0.5)) +
+  theme_classic() +
+  theme(panel.border = element_rect(colour = "black", fill=NA, size=.5) ,
+        axis.title.y = element_text(color="black", size =20),
+        axis.title.x = element_text(color="black", size =20),
+        axis.text= element_text(color="black", size=19), legend.position = "none")
+dev.off()
+
+
+
+################################################################################
+##### SAMPLING COMPLETENESS ANALYSIS
+################################################################################
+
+
+# Loading interaction data for Chao1 estimator
+sampbat<- read.xls("data/sampbat.xlsx", h=T)
+estimateR(sampbat, index =c("chao")) 
+str(sampbat)
+
+samphawk <- read.xls("data/samphawk.xlsx", h=T)
+estimateR(samphawk, index =c("chao"))
+str(samphawk)
+
+# Loading interaction data for rarefaction curve drawing
+sampling_bats <- read.xls("data/sampling_bats.xlsx", h=T)
+curve_bat<- specaccum(sampling_bats, method="rarefaction")
+
+sampling_hawkmoths <- read.xls("data/sampling_hawkmoths.xlsx", h=T)
+curve_hawk<- specaccum(sampling_hawkmoths, method="rarefaction")
+
+# Plot curves
+tiff("figures/sampling.tiff", width = 20, height = 20, units = "cm", res = 600)
+par(mfrow=c(1,2))
+
+plot(curve_hawk, ci.type = "poly", xvar = "individuals", ci.lty=0, ylab = NA,
+     xlab = NA,
+     ci.col=rgb(0.7, 0, 0.2, 0.3), ylim=c(0,30))
+abline(h=22.2, lty=1, col=rgb(0.7, 0, 0.2, 0.3), lwd=2.5)
+abline(h=(22.2+0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+abline(h=(22.2-0.6195203), lty=3, col=rgb(0.7, 0, 0.2, 0.3), lwd=2)
+
+plot(curve_bat, ci.type = "poly", xvar = "individuals", ci.lty=0,
+     ci.col=rgb(0, 0, 0.5, 0.3), ylim=c(0,30), ylab=NA, xlab=NA)
+abline(h=14, lty=1, col=rgb(0, 0, 0.5, 0.3), lwd=2.5)
+abline(h=(14-2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+abline(h=(14+2.283481), lty=3, col=rgb(0, 0, 0.5, 0.3), lwd=2)
+
+dev.off()
\ No newline at end of file
diff --git a/data/morph_graph.xlsx b/data/morph_graph.xlsx
new file mode 100644
index 0000000..1653f30
Binary files /dev/null and b/data/morph_graph.xlsx differ
diff --git a/data/morph_pla.xlsx b/data/morph_pla.xlsx
new file mode 100644
index 0000000..ed32290
Binary files /dev/null and b/data/morph_pla.xlsx differ
diff --git a/data/morph_pol.xlsx b/data/morph_pol.xlsx
new file mode 100644
index 0000000..62f1dba
Binary files /dev/null and b/data/morph_pol.xlsx differ
diff --git a/data/network.txt b/data/network.txt
new file mode 100644
index 0000000..f699a0b
--- /dev/null
+++ b/data/network.txt
@@ -0,0 +1,20 @@
+	Pilo_goun	Bauh_chei	Ceib_glaz	Ench_spec	Ipom_marc	Guet_ange	Toco_form	Ambu_cear	Heli_baru	Crot_sp	Cere_jama	Pseu_marg	Indet1	Hipp_sp	Pipt_stip	Pilo_chry	Aspi_pyri	Comb_sp	Anad_colu	Cord_rig	Allo_sp	Caes_sp	Schi_bras	Indet2
+Call_gris	13	3	0	0	2	26	0	8	0	11	0	0	1	0	4	0	13	6	1	0	3	0	2	1
+Erin_ello	20	18	1	9	0	28	3	3	5	11	0	0	0	0	1	0	15	0	0	1	0	2	0	0
+Lonc_mord	14	9	4	10	18	0	0	0	16	0	0	14	1	2	1	23	0	1	0	0	0	0	0	0
+Glos_sori	23	10	4	22	8	0	0	0	34	4	0	24	0	24	2	27	0	0	0	0	0	0	0	0
+Isog_alla	9	6	0	6	0	4	1	13	1	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0
+Agri_cing	7	2	0	0	4	0	2	0	3	0	5	0	2	2	0	0	0	0	0	0	0	0	0	0
+Eumo_viti	8	6	0	0	0	3	4	0	0	1	0	0	1	0	0	0	0	2	0	0	0	0	0	0
+Eumo_anal	2	0	0	0	0	1	1	0	1	3	3	0	0	0	0	0	0	0	1	0	0	0	0	0
+Mand_rust	6	2	0	0	1	1	3	0	0	0	4	0	0	0	0	0	0	0	0	0	0	0	0	0
+Phyl_disc	0	5	5	8	0	0	0	0	0	0	0	9	0	6	0	0	0	0	0	0	0	0	0	0
+Arti_plan	3	0	2	0	0	0	0	0	0	0	0	2	0	0	0	2	0	0	0	0	0	0	0	0
+Mand_sext	2	1	0	0	0	0	4	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0	0
+Xylo_ters	5	0	3	7	0	0	0	8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
+Prot_stri	6	0	4	8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
+Erin_obsc	2	0	2	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
+Aell_sp	0	0	0	0	0	5	0	15	0	0	0	0	0	0	0	0	10	0	0	0	0	0	0	0
+Cocy_anta	2	1	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
+Call_pars	0	0	0	0	0	0	0	1	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
+Pseu_tetr	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
diff --git a/data/plants.xlsx b/data/plants.xlsx
new file mode 100644
index 0000000..9595ab6
Binary files /dev/null and b/data/plants.xlsx differ
diff --git a/data/sampbat.xlsx b/data/sampbat.xlsx
new file mode 100644
index 0000000..5139e3f
Binary files /dev/null and b/data/sampbat.xlsx differ
diff --git a/data/samphawk.xlsx b/data/samphawk.xlsx
new file mode 100644
index 0000000..43b8098
Binary files /dev/null and b/data/samphawk.xlsx differ
diff --git a/data/sampling_bats.xlsx b/data/sampling_bats.xlsx
new file mode 100644
index 0000000..152df2e
Binary files /dev/null and b/data/sampling_bats.xlsx differ
diff --git a/data/sampling_hawkmoths.xlsx b/data/sampling_hawkmoths.xlsx
new file mode 100644
index 0000000..35d0518
Binary files /dev/null and b/data/sampling_hawkmoths.xlsx differ
diff --git a/figures/bc.tiff b/figures/bc.tiff
new file mode 100644
index 0000000..7a297b0
Binary files /dev/null and b/figures/bc.tiff differ
diff --git a/figures/compound.png b/figures/compound.png
new file mode 100644
index 0000000..9a10ca7
Binary files /dev/null and b/figures/compound.png differ
diff --git a/figures/d.tiff b/figures/d.tiff
new file mode 100644
index 0000000..c5efa96
Binary files /dev/null and b/figures/d.tiff differ
diff --git a/figures/morph.tiff b/figures/morph.tiff
new file mode 100644
index 0000000..095109a
Binary files /dev/null and b/figures/morph.tiff differ
diff --git a/figures/network.tif b/figures/network.tif
new file mode 100644
index 0000000..cece752
Binary files /dev/null and b/figures/network.tif differ
diff --git a/figures/nk.tiff b/figures/nk.tiff
new file mode 100644
index 0000000..8d4fea3
Binary files /dev/null and b/figures/nk.tiff differ
diff --git a/figures/sampling.tiff b/figures/sampling.tiff
new file mode 100644
index 0000000..608fa57
Binary files /dev/null and b/figures/sampling.tiff differ
diff --git a/nestedness.png b/nestedness.png
new file mode 100644
index 0000000..acdad17
Binary files /dev/null and b/nestedness.png differ
diff --git a/results/BCplants.csv b/results/BCplants.csv
new file mode 100644
index 0000000..f7e4c67
--- /dev/null
+++ b/results/BCplants.csv
@@ -0,0 +1,25 @@
+"","betweenness","weighted.betweenness"
+"Pilo_goun",0.112508120964003,0.716112531969309
+"Bauh_chei",0.112508120964003,0.0051150895140665
+"Ceib_glaz",0.035386427298192,0
+"Ench_spec",0.0478641456582633,0.00255754475703325
+"Ipom_marc",0.0728574121956475,0
+"Guet_ange",0.0574119916031681,0.258312020460358
+"Toco_form",0.0313822615293204,0
+"Ambu_cear",0.0574119916031681,0
+"Heli_baru",0.0586358902535373,0.0179028132992327
+"Crot_sp",0.112508120964003,0
+"Cere_jama",0.00505279034690799,0
+"Pseu_marg",0,0
+"Indet1",0.0728574121956475,0
+"Hipp_sp",0.0101381461675579,0
+"Pipt_stip",0.0945156967215791,0
+"Pilo_chry",0,0
+"Aspi_pyri",0.0430960379489791,0
+"Comb_sp",0.0585414585414585,0
+"Anad_colu",0.0173239750445633,0
+"Cord_rig",0,0
+"Allo_sp",0,0
+"Caes_sp",0,0
+"Schi_bras",0,0
+"Indet2",0,0
diff --git a/results/NDplants.csv b/results/NDplants.csv
new file mode 100644
index 0000000..7f801a0
--- /dev/null
+++ b/results/NDplants.csv
@@ -0,0 +1,25 @@
+"","x"
+"Pilo_goun",0.842105263157895
+"Bauh_chei",0.578947368421053
+"Ceib_glaz",0.421052631578947
+"Ench_spec",0.368421052631579
+"Ipom_marc",0.368421052631579
+"Guet_ange",0.368421052631579
+"Toco_form",0.368421052631579
+"Ambu_cear",0.315789473684211
+"Heli_baru",0.315789473684211
+"Crot_sp",0.315789473684211
+"Cere_jama",0.263157894736842
+"Pseu_marg",0.210526315789474
+"Indet1",0.210526315789474
+"Hipp_sp",0.210526315789474
+"Pipt_stip",0.210526315789474
+"Pilo_chry",0.157894736842105
+"Aspi_pyri",0.157894736842105
+"Comb_sp",0.157894736842105
+"Anad_colu",0.105263157894737
+"Cord_rig",0.0526315789473684
+"Allo_sp",0.0526315789473684
+"Caes_sp",0.0526315789473684
+"Schi_bras",0.0526315789473684
+"Indet2",0.0526315789473684
diff --git a/results/dplants.csv b/results/dplants.csv
new file mode 100644
index 0000000..fcf0a73
--- /dev/null
+++ b/results/dplants.csv
@@ -0,0 +1,25 @@
+"","d"
+"Pilo_goun",0.0870684996573851
+"Bauh_chei",0.109869600599462
+"Ceib_glaz",0.285666228720193
+"Ench_spec",0.208445059351663
+"Ipom_marc",0.265958761343998
+"Guet_ange",0.366920370363653
+"Toco_form",0.407610228547915
+"Ambu_cear",0.500337014187968
+"Heli_baru",0.228839087598398
+"Crot_sp",0.253786440481666
+"Cere_jama",0.540421123018025
+"Pseu_marg",0.306313815816232
+"Indet1",0.232476097455208
+"Hipp_sp",0.308948643771313
+"Pipt_stip",0.0763551585328956
+"Pilo_chry",0.342617618479856
+"Aspi_pyri",0.400626621068713
+"Comb_sp",0.291223767218696
+"Anad_colu",0.294794896990733
+"Cord_rig",0.0849024182417464
+"Allo_sp",0.29669586237059
+"Caes_sp",0.18343307057735
+"Schi_bras",0.227357450604607
+"Indet2",0.126962155226962
diff --git a/results/results_centrality.txt b/results/results_centrality.txt
new file mode 100644
index 0000000..1a6b8f8
--- /dev/null
+++ b/results/results_centrality.txt
@@ -0,0 +1,55 @@
+[1] "Comparison of centrality by syndrome"
+[1] "Queiroz et al. 2020, Biotropica"
+
+ [1] ""                                                                                
+ [2] "Call:  glm(formula = plants$d ~ plants$Guild, family = quasibinomial(\"logit\"))"
+ [3] ""                                                                                
+ [4] "Coefficients:"                                                                   
+ [5] "      (Intercept)  plants$Guildother  plants$Guildsphin  "                       
+ [6] "          -1.1626            -0.1409             0.9344  "                       
+ [7] ""                                                                                
+ [8] "Degrees of Freedom: 21 Total (i.e. Null);  19 Residual"                          
+ [9] "Null Deviance:\t    1.729 "                                                      
+[10] "Residual Deviance: 0.861 \tAIC: NA"                                              
+ [1] "Analysis of Deviance Table"                                     ""                                                              
+ [3] "Model: quasibinomial, link: logit"                              ""                                                              
+ [5] "Response: plants$d"                                             ""                                                              
+ [7] "Terms added sequentially (first to last)"                       ""                                                              
+ [9] ""                                                               "             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    "   
+[11] "NULL                            21    1.72876              "    "plants$Guild  2  0.86774        19    0.86102 2.664e-05 ***"   
+[13] "---"                                                            "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1"
+
+ [1] ""                                                                                 
+ [2] "Call:  glm(formula = plants$bc ~ plants$Guild, family = quasibinomial(\"logit\"))"
+ [3] ""                                                                                 
+ [4] "Coefficients:"                                                                    
+ [5] "      (Intercept)  plants$Guildother  plants$Guildsphin  "                        
+ [6] "          -2.9447            -0.3615            -0.2632  "                        
+ [7] ""                                                                                 
+ [8] "Degrees of Freedom: 21 Total (i.e. Null);  19 Residual"                           
+ [9] "Null Deviance:\t    0.9808 "                                                      
+[10] "Residual Deviance: 0.9569 \tAIC: NA"                                              
+ [1] "Analysis of Deviance Table"                             ""                                                      
+ [3] "Model: quasibinomial, link: logit"                      ""                                                      
+ [5] "Response: plants$bc"                                    ""                                                      
+ [7] "Terms added sequentially (first to last)"               ""                                                      
+ [9] ""                                                       "             Df Deviance Resid. Df Resid. Dev Pr(>Chi)"
+[11] "NULL                            21    0.98081         " "plants$Guild  2 0.023963        19    0.95685   0.7577"
+
+ [1] ""                                                                                 
+ [2] "Call:  glm(formula = plants$nk ~ plants$Guild, family = quasibinomial(\"logit\"))"
+ [3] ""                                                                                 
+ [4] "Coefficients:"                                                                    
+ [5] "      (Intercept)  plants$Guildother  plants$Guildsphin  "                        
+ [6] "          -0.4643            -1.3664            -0.4082  "                        
+ [7] ""                                                                                 
+ [8] "Degrees of Freedom: 21 Total (i.e. Null);  19 Residual"                           
+ [9] "Null Deviance:\t    3.643 "                                                       
+[10] "Residual Deviance: 2.252 \tAIC: NA"                                               
+ [1] "Analysis of Deviance Table"                                     ""                                                              
+ [3] "Model: quasibinomial, link: logit"                              ""                                                              
+ [5] "Response: plants$nk"                                            ""                                                              
+ [7] "Terms added sequentially (first to last)"                       ""                                                              
+ [9] ""                                                               "             Df Deviance Resid. Df Resid. Dev Pr(>Chi)   "     
+[11] "NULL                            21     3.6428            "      "plants$Guild  2   1.3911        19     2.2518 0.002469 **"     
+[13] "---"                                                            "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1"
diff --git a/results/results_topology.txt b/results/results_topology.txt
new file mode 100644
index 0000000..896e5e9
--- /dev/null
+++ b/results/results_topology.txt
@@ -0,0 +1,19 @@
+[1] "Topological analysis of the nocturnal pollination network"
+[1] "Queiroz et al. 2020, Biotropica"
+
+[1] "The network has 19 rows and 24 columns."
+
+[1] "The network's specialization (H2) is 0.34 , P = 0"
+
+[1] "The network's modularity (DIRT_LPA+) is 0.36 , P = 0 , and it contains 4 modules."
+
+[1] "The network's nestedness (WNODF) is 0.34 , P = 1"
+
+[1] "The network shows the following scores of nestedness (WNODA):"
+
+[1] "Entire network = 0.38 , P = 0.32"
+
+[1] "Between the modules = 0.31 , P = 0.38"
+
+[1] "Within the modules = 0.56 , P = 0.31"
+
diff --git a/specialization.png b/specialization.png
new file mode 100644
index 0000000..18bf3fe
Binary files /dev/null and b/specialization.png differ