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


<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
#  fig.path = "man/figures/README-",
  out.width = "100%"
)
```


<!-- badges: start -->
<!-- badges: end -->


 
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