First commit

First upload of all files to the repo.

.gitignore

@@ -0,0 +1,2 @@
+
+.DS_Store

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))

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)
+}
+
+

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*.

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)
+}
+
+
+
+