--- title: "Analysis of a mutualistic multipartite ecological network with GREMLINS" author: "Sophie Donnet, Pierre Barbillon" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 4 bibliography: biblio.bib link-citations: yes vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{Analysis of a mutualistic multipartite ecological network with GREMLINS} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", # fig.path = "man/figures/README-", out.width = "100%" ) ``` Note that the `sbm` package is much more easy to use (but implements the same inference methods and algorithm). The same dataset is presented in a vignette. ## The dataset We apply our methodology to an ecological mutualistic multipartite network. The dataset --compiled and conducted by @Dattilo at Centro de Investigaciones Costeras La Mancha (CICOLMA), located on the central coast of the Gulf of Mexico, Veracruz, Mexico-- involves three general types of plant-animal mutualistic interaction: pollination, seed dispersal by frugivorous birds, and protective mutualisms between ants and plants with extrafloral nectaries. The dataset --which is one of the largest compiled so far with respect to species richness, number of interactions and sampling effort-- includes 4 functional groups, namely plants, pollinator species (refered as floral visitors), ant species and frugivorous bird species. Three binary bipartite networks have been collected representing interactions between 1/ plants and florals visitor, 2/ plants and ants, and 3/ plants and seed dispersal birds, resulting into three bipartite networks. The FG are of respective sizes: $n_1 = 141$ plant species, $n_2 = 173$ pollinator species (refered as ), $n_3 = 46$ frugivorous bird species and $n_4 = 30$ ant species. The 3 networks contain $753$ observed interactions of which $55\%$ are plant-pollinator interactions, $17\%$ are plant-birds interactions and $28\%$ are plant-ant interactions. ```{r loading dataset, eval=TRUE} library(GREMLINS) data(MPEcoNetwork, package = "GREMLINS") names(MPEcoNetwork) ``` As required by GREMLINS, our the global network has to be encoded in separate matrices for each network (in our case the $3$ incidence matrices) So, here, our 3 networks are provided in 3 incidence matrices, the plants being in rows. *Note that the order of the individuals within the functional groups must be the same in all the matrices*. ## Formatting the data We format the data to be able to use our R package GREMLINS i.e. we transform the matrices into an list containing *the matrix*, *its type* : `inc` for incidence matrix, `adj` for adjacency symetric, and `diradj` for non symetric (oriented) adjacency matrix, the name of functional group in row and the name of functional group in column. The three matrices are gathered in a list. To do so, we use de the function `defineNetwork`. ```{r transform dataset, eval=TRUE} PlantFlovis = defineNetwork(MPEcoNetwork$Inc_plant_flovis,"inc","Plants","Flovis") PlantAnt = defineNetwork(MPEcoNetwork$Inc_plant_ant,"inc","Plants","Ants") PlantBird = defineNetwork(MPEcoNetwork$Inc_plant_bird,"inc","Plants","Birds") list_net <- list(PlantFlovis,PlantAnt,PlantBird) names(PlantFlovis) ``` If one wants to keep a track of the names of the species, they should be used as rownames and colnames in the matrices. ```{r example of dataset, eval=TRUE} PlantFlovis$mat[1:2,1:2] ``` ## Inference The model selection and the estimation are performed with the function `multipartiteBM`. ```{r MBM, echo = TRUE, eval = FALSE} RES_MBM = multipartiteBM( list_Net = list(PlantFlovis, PlantAnt, PlantBird), namesFG = c('Plants','Flovis','Ants','Birds'), v_distrib = c('bernoulli','bernoulli','bernoulli'), initBM = TRUE, keep = TRUE, nbCores = 2) ``` ```{r MBM load, echo = FALSE, eval = TRUE} load(file='res_EcologicalNetwork.Rda') ``` RES_MBM contains the estimated parameters of the models we run through during the search of the better numbers of blocks. If one sets `keep = FALSE` in the `multipartiteBM` function then we only save the best model. RES_MBM constains de dataset and the results. ```{r MBM what} names(RES_MBM) ``` The better model has the following numbers of blocks ```{r MBM v_K } RES_MBM$fittedModel[[1]]$paramEstim$v_K ``` To see the parameters estimated for the better model we use the following command `RES_MBM$fittedModel[[1]]$paramEstim$***` ```{r MBM param } RES_MBM$fittedModel[[1]]$paramEstim$list_pi$Plants RES_MBM$fittedModel[[1]]$paramEstim$list_theta$PlantsFlovis ``` The clustering supplied by the better model are in `RES_MBM$fittedModel[[1]]$paramEstim$Z$***`. ```{r MBM Z } table(RES_MBM$fittedModel[[1]]$paramEstim$Z$Plants) table(RES_MBM$fittedModel[[1]]$paramEstim$Z$Ants) ``` ## Plots Please use the plot functions included in the R package `sbm`. ## References