Advanced metrics on networks and null models
Outline
- Advanced metrics
1.1. degrees
1.2. trophic levels
1.3. generalism, network specialization
1.4. modularity, clustering, nestedness
1.5. robustness to species removal
1.6. \(\beta\) diversity of networks
- Models
2.1. de novo generation of networks
2.2. null models
Preamble
The script script_metricsnull_v2.R contains the examples given in the presentation. The source script functions.R contains the functions especially devised for this course.
load("mg1.Rdata")
source('functions.R')We will make use of two real datasets: Food web data from Cruz-Escalona, V. H., Arreguin-Sanchez, F. & Zetina-Rejon, M. (2007) Analysis of the ecosystem structure of Laguna Alvarado, western Gulf of Mexico, by means of a mass balance model. Estuarine, Coastal and Shelf Science, 72, 155-167.
mat_foodweb<-t(as.matrix(mg1[[1]][[433]]))
rownames(mat_foodweb)<-names(mg1[[1]][[433]])
colnames(mat_foodweb)<-names(mg1[[1]][[433]])Plant-pollinator data from Barrett, S. C. H. & Helenurm, K. (1987) The reproductive biology of boreal forest herbs. I. Breeding systems and pollination. Canadian Journal of Botany, 65, 2036-2046.
mat_plantpol<-barrett1987
mat_plantpol_bin<-mat_plantpol
mat_plantpol_bin[mat_plantpol>0]<-1Building graphs
One of the most useful R objects to work on networks is the graph, created through the igraph package. Here is a collection of examples, using functions from different packages for graphical display.
foodweb<-graph_from_adjacency_matrix(mat_foodweb)
foodweb<-largest_component(foodweb)
plot(foodweb,layout=layout_as_food_web3(foodweb))
plotMyMatrix(as_adj(foodweb,sparse=FALSE))
undirected_foodweb<-as.undirected(foodweb)
plotMyMatrix(as_adj(undirected_foodweb,sparse=FALSE))
pollination<-graph_from_biadjacency_matrix(mat_plantpol)
plot(pollination,layout=layout_as_bipartite)
plotMyMatrix(as_adj(pollination,sparse=FALSE))
plotMyMatrix(as_biadjacency_matrix(pollination,sparse=FALSE))
pollination_bin<-graph_from_biadjacency_matrix(mat_plantpol_bin)
plotweb(mat_plantpol)
In the following: - foodweb is a directed binary network - undirected_foodweb is its undirected equivalent (the adjacency matrix has been made symmetric) - pollination is a weighted bipartite network - pollination_bin is its binarized version
Advanced metrics
Built-in metrics
The functions salvaged from FWebs and the bipartite package both yield lots of ‘’built-in automatic outputs’’, accessible through the following commands for instance:
fw.metrics(list(list(mat_foodweb)))## Computing network metrics for food web 1 of 1...## $number_nodes
## [1] 106
##
## $number_links
## [1] 233
##
## $link_density
## [1] 2.198113
##
## $connectance
## [1] 0.02093441
##
## $compartmentalization
## [1] NaN
##
## $maximum_trophic_level
## [1] 3.84058networklevel(mat_plantpol)## connectance web asymmetry links per species number of compartments
## 0.13643791 0.78947368 1.46491228 2.00000000
## compartment diversity cluster coefficient modularity Q nestedness
## 1.09233738 0.08333333 0.56737521 10.23333209
## NODF weighted nestedness weighted NODF spectral radius
## 30.78467181 0.66719667 13.49547377 72.55427745
## interaction strength asymmetry specialisation asymmetry linkage density weighted connectance
## 0.23565934 -0.20755966 9.29902902 0.08157043
## Fisher alpha Shannon diversity interaction evenness Alatalo interaction evenness
## 81.60982933 4.11901613 0.57933698 0.38703061
## H2 number.of.species.HL number.of.species.LL mean.number.of.shared.partners.HL
## 0.57807328 102.00000000 12.00000000 0.56105611
## mean.number.of.shared.partners.LL cluster.coefficient.HL cluster.coefficient.LL weighted.cluster.coefficient.HL
## 2.07575758 0.29803030 0.28659537 0.83564121
## weighted.cluster.coefficient.LL niche.overlap.HL niche.overlap.LL togetherness.HL
## 0.17336473 0.34762949 0.13855758 0.14499265
## togetherness.LL C.score.HL C.score.LL V.ratio.HL
## 0.05157815 0.50388420 0.59112597 1.54171552
## V.ratio.LL discrepancy.HL discrepancy.LL extinction.slope.HL
## 27.81509917 71.00000000 63.00000000 5.84839013
## extinction.slope.LL robustness.HL robustness.LL functional.complementarity.HL
## 1.54269007 0.82379068 0.52452619 382.39912254
## functional.complementarity.LL partner.diversity.HL partner.diversity.LL generality.HL
## 338.58919107 0.68675599 2.17835524 2.31225322
## vulnerability.LL
## 16.28580482The issue is that not all these outputs are equally useful… Let’s dissect them one by one.
Trophic levels
Trophic levels classically correspond to the number of trophic links between a species and any basal species (primary producers) + 1. Basal species always belong to level 1 and consumers of basal species belong to species 2, etc. Complexity arises when species are omnivorous and can feed on multiple trophic levels. The package cheddar has methods to compute all the classic trophic levels. For instance, the first column of foodweb_TL gives trophic levels as computed using the ‘’shortest path’’ method (TL of a species is exactly 1+the size of the shortest path to any basal species) while the third column yields the longest trophic level (the longest path linking the focal species to any basal species).
foodweb_TL<-TrophicLevels(as_Community(foodweb))
plotMyMatrix(as_adj(foodweb,sparse=FALSE),clustering=list("row"=foodweb_TL[,1],"col"=foodweb_TL[,1]))
plotMyMatrix(as_adj(foodweb,sparse=FALSE),clustering=list("row"=foodweb_TL[,3],"col"=foodweb_TL[,3]))
Another possible measure of trophic levels is the top-down measure of MacKay et al., which is based on the use of a normalized Laplacian of the graph. It yields a measure of trophic level that is maximal for apex predators and decreases according to the distance between species and the apex predators (more or less).
tl.1<-trophic_levels(foodweb)
plot(TrophicLevels(as_Community(foodweb,"."))[,1],tl.1[,1],xlab="ShortestTL",ylab="MacKayTL")
plot(TrophicLevels(as_Community(foodweb,"."))[,6],tl.1[,1],xlab="PreyAveragedTL",ylab="MacKayTL")
These plots show the correspondence between the classic TL 1 (shortest path) and 6 (prey-averaged TL) and MacKay’s TL definition.
Network structure
Ecological interaction networks can have non-random structures. In ecology, we tend to focus on two such structures: modularity (on the left) and nestedness (on the right).
Node specialization
In ecology, one nagging question is whether we can reliably tell the level of specialization of the different species, i.e. whether they depend on a single species or can substitute one interactor species for another one.
Blüthgen et al. (2006) proposed a metric (often noted \(d_i\) or \(d^\prime_i\) in its relative version) to work out the level of specialisation of a species based on weighted interactions, using the following equation: \[d_i = \sum_j \frac{w_{ij}}{W_{i\cdot}} ln\left[ \frac{w_{ij}W_{\cdot \cdot}}{W_{i\cdot}W_{\cdot j}}\right]\] where \(w_{ij}\) is the weight of the interaction between species \(i\) and \(j\), \(W_{i\cdot}\) is the sum of all interactions of species \(i\) (row sum), \(W_{\cdot j}\) is the sum of all interactions of species \(j\) (column sum), and \(W_{\cdot \cdot}\) is the sum of all interaction weights between all species.
The following formula yields the relative specialization score, which is always between 0 and 1: \[d^\prime_i = \frac{d_i-d_{min}}{d_{max}-d_{min}}\] where the \(d_{min}\) and \(d_{max}\) values are obtained from theoretical minimum and maximum.
The dfun function in the package bipartite does the computation of both \(d\) and \(d^\prime\), but the computation of \(d_{min}\) is not always accurate. One can check the relationship between specialization and degree with a simple scatter plot.
dfun(mat_plantpol)## $dprime
## Aralia.nudicaulis Chimaphila.umbellata Clintonia.borealis Cornus.canadensis Cypripedium.acaule Linnaea.borealis
## 0.4029868 0.8734263 0.4617572 0.5573480 0.1978600 0.5606655
## Maianthemum.canadense Medeola.virginiana Oxalis.montana Pyrola.secunda Trientalis.borealis Trillium.undulatum
## 0.6348963 1.0000000 0.7139439 0.5462431 0.7462165 0.4919476
##
## $d
## Aralia.nudicaulis Chimaphila.umbellata Clintonia.borealis Cornus.canadensis Cypripedium.acaule Linnaea.borealis
## 0.9987162 1.8296061 1.6615624 0.7119762 1.3358969 1.6530716
## Maianthemum.canadense Medeola.virginiana Oxalis.montana Pyrola.secunda Trientalis.borealis Trillium.undulatum
## 1.2600139 6.3099183 1.4272914 3.0608135 4.1383474 2.8379140
##
## $dmin
## Aralia.nudicaulis Chimaphila.umbellata Clintonia.borealis Cornus.canadensis Cypripedium.acaule Linnaea.borealis
## 0.2416664 0.2299450 0.4001940 0.1076680 0.6219049 0.3182885
## Maianthemum.canadense Medeola.virginiana Oxalis.montana Pyrola.secunda Trientalis.borealis Trillium.undulatum
## 0.2040323 1.6655274 0.2088886 0.8183194 0.9834566 0.8183194
##
## $dmax
## Aralia.nudicaulis Chimaphila.umbellata Clintonia.borealis Cornus.canadensis Cypripedium.acaule Linnaea.borealis
## 2.120264 2.061423 3.131864 1.191924 4.230477 2.699000
## Maianthemum.canadense Medeola.virginiana Oxalis.montana Pyrola.secunda Trientalis.borealis Trillium.undulatum
## 1.867267 6.309918 1.915469 4.923624 5.211306 4.923624plot(degree(pollination)[!V(pollination_bin)$type],dfun(mat_plantpol)$dprime,xlab="degrees",ylab="d'")
This function was coded for bipartite networks, but you need to transpose the bi-adjacency matrix to obtain the values of \(d\) or \(d^\prime\) for the other level.
dfun(t(mat_plantpol))## $dprime
## Acmaeopsoides.rufula Agiotes.stabilis Ancistrocerus.sp. Andrena.melanochroa Andrena.miranda
## 0.00000000 0.21563785 0.07436233 0.13593241 0.04749160
## Andrena.nivalis Andrena.rufosignata Andrena.sigmundi Andrena.thaspii Andrena.wheeleri
## 0.00000000 0.07436233 0.00000000 0.00000000 0.18138730
## Andrena.wilkella Andrena.w.scripta Anthaxia.expansa Anthonomus.sp. Blera.badia
## 0.37904305 0.19882256 0.26144023 0.00000000 0.00000000
## Blera.confusa Blera.nigripes Bombus.perplexus Bombus.ternarius Bombus.terricola
## 0.28139346 0.00000000 0.19070143 0.27712368 0.17343875
## Bombus.vagans Bombylius.major Bombylius.pigmaeus Cantharis.sp. Carposcalis.confusus
## 0.15129125 0.16064599 0.13195455 0.00000000 0.38747019
## Carposcalis.obscura Cateres.pennatus Celastrina.argiolus Chalcosyrphus.curvarius Chalcosyrphus.inarmatus
## 0.78534318 0.18138730 0.07436233 0.00000000 0.00000000
## Chalcosyrphus.memorum Chalcosyrphus.nigra Cheilosia.pallipes Cheilosia.sialia Cheilosia.slossonae
## 0.00000000 0.07436233 0.22572129 0.07436233 0.13195455
## Cheilosia.tristis Crabro.sp. Dalopius.sp. Dialictus.cressoni Dialictus.imitatus
## 0.18138730 0.08699553 0.00000000 0.20792957 0.22937119
## Dialictus.sp. Eclimus.harrisi Ephemerella.sp. Eremomyioides.setosa Eudasyphora.cyanicolor
## 0.16289528 0.30481234 1.00000000 0.13195455 0.00000000
## Eupyses.vestris Eusphalerum.convexum Eusphalerum.phothos Eusphalerum.sp. Evodinus.monticola
## 0.29446614 0.38985945 0.32822619 0.70755994 0.00000000
## Evylaeus.divergens Evylaeus.quebecensis Fannia.ungulata Halictus.confusus Halictus.rubicundus
## 0.00000000 0.14829973 0.38962890 0.00000000 0.09710005
## Hoplitis.producta Hylaeus.basalis Hylaeus.ellipticus Hylaeus.stevensi Judolia.montivagans
## 0.72913325 0.00000000 0.55310629 0.29446614 0.00000000
## Lygus.refidorsus Mallota.posticata Melagyna.lasiophthalma Melanostoma.sp. Nomada.cuneata
## 0.00000000 0.02646228 0.15763329 0.22572129 0.13195455
## Nomada.sayi Orsodacne.atra Orthoneura.pulchella Osmia.lignaria Osmia.proxima
## 0.00000000 0.00000000 0.04749160 0.37904305 0.10709229
## Osmia.pumila Paramyia.nitens Parasyrphus.relictus Passaloecus.sp. Pegohylemyia.fugax
## 0.29446614 0.89940438 0.51915999 0.29446614 0.13195455
## Phalaenophana.pyramusalis Phaonia.serva Pidonia.ruficolis Pieris.napi Poanes.hobomok
## 0.13195455 0.21508697 0.23907334 0.00000000 0.37904305
## Psithyrus.sp. Pyrophaena.rosarum Rhopalum.sp. Sarcophaga.nearctica Sphaerophoria.bifurcata
## 0.00000000 0.29446614 0.37904305 0.13227114 0.30083960
## Sphaerophoria.longipilosa Sphaerophoria.sp. Sphecodes.sp. Syritta.pipiens Syritta.rectus
## 0.00000000 0.19515964 0.07436233 0.00000000 0.00000000
## Temnostoma.alternans Temnostoma.balyras Temnostoma.barberi Thamnotettix.confinis Trachysida.aspera
## 0.00000000 0.13593241 0.07436233 0.13195455 0.18138730
## Tropidia.quadrata Tychius.stephensi Vespula.arenaria Xylota.bigelowi Xylota.hinei
## 0.10189683 0.13593241 0.72913325 0.28139346 0.04749160
## Xylota.sp. Zeraea.americana
## 0.00000000 0.08699553
##
## $d
## Acmaeopsoides.rufula Agiotes.stabilis Ancistrocerus.sp. Andrena.melanochroa Andrena.miranda
## 1.1919245 1.8672670 1.1919245 1.1919245 0.7805245
## Andrena.nivalis Andrena.rufosignata Andrena.sigmundi Andrena.thaspii Andrena.wheeleri
## 0.8364486 1.1919245 1.1919245 1.1919245 2.1202635
## Andrena.wilkella Andrena.w.scripta Anthaxia.expansa Anthonomus.sp. Blera.badia
## 3.1318644 0.9311967 1.5625142 1.1919245 1.1919245
## Blera.confusa Blera.nigripes Bombus.perplexus Bombus.ternarius Bombus.terricola
## 1.1919245 1.1919245 1.0759162 1.8487036 0.6982811
## Bombus.vagans Bombylius.major Bombylius.pigmaeus Cantharis.sp. Carposcalis.confusus
## 0.4146760 0.7984489 1.8672670 1.1919245 1.9571744
## Carposcalis.obscura Cateres.pennatus Celastrina.argiolus Chalcosyrphus.curvarius Chalcosyrphus.inarmatus
## 5.2113060 2.1202635 1.1919245 1.1919245 1.1919245
## Chalcosyrphus.memorum Chalcosyrphus.nigra Cheilosia.pallipes Cheilosia.sialia Cheilosia.slossonae
## 1.1919245 1.1919245 1.9154691 1.1919245 1.8672670
## Cheilosia.tristis Crabro.sp. Dalopius.sp. Dialictus.cressoni Dialictus.imitatus
## 2.1202635 1.2523152 1.1919245 0.9686997 1.9329168
## Dialictus.sp. Eclimus.harrisi Ephemerella.sp. Eremomyioides.setosa Eudasyphora.cyanicolor
## 0.4403134 1.7598962 6.3099183 1.8672670 1.1919245
## Eupyses.vestris Eusphalerum.convexum Eusphalerum.phothos Eusphalerum.sp. Evodinus.monticola
## 2.6990004 1.7178941 1.2370680 1.1798908 1.1919245
## Evylaeus.divergens Evylaeus.quebecensis Fannia.ungulata Halictus.confusus Halictus.rubicundus
## 0.8364486 0.6271943 2.6990004 1.1919245 1.3006181
## Hoplitis.producta Hylaeus.basalis Hylaeus.ellipticus Hylaeus.stevensi Judolia.montivagans
## 4.9236239 1.1919245 2.6990004 2.6990004 1.1919245
## Lygus.refidorsus Mallota.posticata Melagyna.lasiophthalma Melanostoma.sp. Nomada.cuneata
## 1.1919245 0.9629468 1.5899865 1.9154691 1.8672670
## Nomada.sayi Orsodacne.atra Orthoneura.pulchella Osmia.lignaria Osmia.proxima
## 1.1919245 1.1919245 0.7805245 3.1318644 1.0577689
## Osmia.pumila Paramyia.nitens Parasyrphus.relictus Passaloecus.sp. Pegohylemyia.fugax
## 2.6990004 1.9648339 1.2565624 2.6990004 1.8672670
## Phalaenophana.pyramusalis Phaonia.serva Pidonia.ruficolis Pieris.napi Poanes.hobomok
## 1.8672670 1.5601276 1.0586413 1.1919245 3.1318644
## Psithyrus.sp. Pyrophaena.rosarum Rhopalum.sp. Sarcophaga.nearctica Sphaerophoria.bifurcata
## 1.1919245 2.6990004 3.1318644 1.4687473 1.5691864
## Sphaerophoria.longipilosa Sphaerophoria.sp. Sphecodes.sp. Syritta.pipiens Syritta.rectus
## 1.1919245 1.4674317 1.1919245 1.1919245 1.1919245
## Temnostoma.alternans Temnostoma.balyras Temnostoma.barberi Thamnotettix.confinis Trachysida.aspera
## 1.1919245 1.1919245 1.1919245 1.8672670 2.1202635
## Tropidia.quadrata Tychius.stephensi Vespula.arenaria Xylota.bigelowi Xylota.hinei
## 0.8364486 1.1919245 4.9236239 1.1919245 0.7805245
## Xylota.sp. Zeraea.americana
## 1.1919245 1.2523152
##
## $dmin
## Acmaeopsoides.rufula Agiotes.stabilis Ancistrocerus.sp. Andrena.melanochroa Andrena.miranda
## 1.191924466 0.836448563 0.836448563 0.559607915 0.559607915
## Andrena.nivalis Andrena.rufosignata Andrena.sigmundi Andrena.thaspii Andrena.wheeleri
## 0.836448563 0.836448563 1.191924466 1.191924466 1.191924466
## Andrena.wilkella Andrena.w.scripta Anthaxia.expansa Anthonomus.sp. Blera.badia
## 1.191924466 0.112437627 0.372726551 1.191924466 1.191924466
## Blera.confusa Blera.nigripes Bombus.perplexus Bombus.ternarius Bombus.terricola
## 0.089460825 1.191924466 0.221831524 0.559607915 0.058989497
## Bombus.vagans Bombylius.major Bombylius.pigmaeus Cantharis.sp. Carposcalis.confusus
## 0.025770376 0.116027905 1.191924466 1.191924466 0.221831524
## Carposcalis.obscura Cateres.pennatus Celastrina.argiolus Chalcosyrphus.curvarius Chalcosyrphus.inarmatus
## 1.191924466 1.191924466 0.836448563 1.191924466 1.191924466
## Chalcosyrphus.memorum Chalcosyrphus.nigra Cheilosia.pallipes Cheilosia.sialia Cheilosia.slossonae
## 1.191924466 0.836448563 0.836448563 0.836448563 1.191924466
## Cheilosia.tristis Crabro.sp. Dalopius.sp. Dialictus.cressoni Dialictus.imitatus
## 1.191924466 0.836448563 1.191924466 0.112437627 0.836448563
## Dialictus.sp. Eclimus.harrisi Ephemerella.sp. Eremomyioides.setosa Eudasyphora.cyanicolor
## 0.025452320 0.372726551 1.191924466 1.191924466 1.191924466
## Eupyses.vestris Eusphalerum.convexum Eusphalerum.phothos Eusphalerum.sp. Evodinus.monticola
## 1.191924466 0.112437627 0.047925742 0.004890966 1.191924466
## Evylaeus.divergens Evylaeus.quebecensis Fannia.ungulata Halictus.confusus Halictus.rubicundus
## 0.836448563 0.070383992 0.836448563 1.191924466 0.836448563
## Hoplitis.producta Hylaeus.basalis Hylaeus.ellipticus Hylaeus.stevensi Judolia.montivagans
## 1.191924466 1.191924466 0.221831524 1.191924466 1.191924466
## Lygus.refidorsus Mallota.posticata Melagyna.lasiophthalma Melanostoma.sp. Nomada.cuneata
## 1.191924466 0.836448563 0.836448563 0.836448563 1.191924466
## Nomada.sayi Orsodacne.atra Orthoneura.pulchella Osmia.lignaria Osmia.proxima
## 1.191924466 1.191924466 0.559607915 1.191924466 0.559607915
## Osmia.pumila Paramyia.nitens Parasyrphus.relictus Passaloecus.sp. Pegohylemyia.fugax
## 1.191924466 0.016188829 0.024272240 1.191924466 1.191924466
## Phalaenophana.pyramusalis Phaonia.serva Pidonia.ruficolis Pieris.napi Poanes.hobomok
## 1.191924466 0.559607915 0.099097400 1.191924466 1.191924466
## Psithyrus.sp. Pyrophaena.rosarum Rhopalum.sp. Sarcophaga.nearctica Sphaerophoria.bifurcata
## 1.191924466 1.191924466 1.191924466 0.836448563 0.221831524
## Sphaerophoria.longipilosa Sphaerophoria.sp. Sphecodes.sp. Syritta.pipiens Syritta.rectus
## 1.191924466 0.559607915 0.836448563 1.191924466 1.191924466
## Temnostoma.alternans Temnostoma.balyras Temnostoma.barberi Thamnotettix.confinis Trachysida.aspera
## 1.191924466 0.559607915 0.836448563 1.191924466 1.191924466
## Tropidia.quadrata Tychius.stephensi Vespula.arenaria Xylota.bigelowi Xylota.hinei
## 0.372726551 0.559607915 1.191924466 0.089460825 0.559607915
## Xylota.sp. Zeraea.americana
## 1.191924466 0.836448563
##
## $dmax
## Acmaeopsoides.rufula Agiotes.stabilis Ancistrocerus.sp. Andrena.melanochroa Andrena.miranda
## 6.309918 5.616771 5.616771 5.211306 5.211306
## Andrena.nivalis Andrena.rufosignata Andrena.sigmundi Andrena.thaspii Andrena.wheeleri
## 5.616771 5.616771 6.309918 6.309918 6.309918
## Andrena.wilkella Andrena.w.scripta Anthaxia.expansa Anthonomus.sp. Blera.badia
## 6.309918 4.230477 4.923624 6.309918 6.309918
## Blera.confusa Blera.nigripes Bombus.perplexus Bombus.ternarius Bombus.terricola
## 4.007333 6.309918 4.700480 5.211306 3.744969
## Bombus.vagans Bombylius.major Bombylius.pigmaeus Cantharis.sp. Carposcalis.confusus
## 2.596346 4.364008 6.309918 6.309918 4.700480
## Carposcalis.obscura Cateres.pennatus Celastrina.argiolus Chalcosyrphus.curvarius Chalcosyrphus.inarmatus
## 6.309918 6.309918 5.616771 6.309918 6.309918
## Chalcosyrphus.memorum Chalcosyrphus.nigra Cheilosia.pallipes Cheilosia.sialia Cheilosia.slossonae
## 6.309918 5.616771 5.616771 5.616771 6.309918
## Cheilosia.tristis Crabro.sp. Dalopius.sp. Dialictus.cressoni Dialictus.imitatus
## 6.309918 5.616771 6.309918 4.230477 5.616771
## Dialictus.sp. Eclimus.harrisi Ephemerella.sp. Eremomyioides.setosa Eudasyphora.cyanicolor
## 2.572249 4.923624 6.309918 6.309918 6.309918
## Eupyses.vestris Eusphalerum.convexum Eusphalerum.phothos Eusphalerum.sp. Evodinus.monticola
## 6.309918 4.230477 3.670861 1.665527 6.309918
## Evylaeus.divergens Evylaeus.quebecensis Fannia.ungulata Halictus.confusus Halictus.rubicundus
## 5.616771 3.825012 5.616771 6.309918 5.616771
## Hoplitis.producta Hylaeus.basalis Hylaeus.ellipticus Hylaeus.stevensi Judolia.montivagans
## 6.309918 6.309918 4.700480 6.309918 6.309918
## Lygus.refidorsus Mallota.posticata Melagyna.lasiophthalma Melanostoma.sp. Nomada.cuneata
## 6.309918 5.616771 5.616771 5.616771 6.309918
## Nomada.sayi Orsodacne.atra Orthoneura.pulchella Osmia.lignaria Osmia.proxima
## 6.309918 6.309918 5.211306 6.309918 5.211306
## Osmia.pumila Paramyia.nitens Parasyrphus.relictus Passaloecus.sp. Pegohylemyia.fugax
## 6.309918 2.182784 2.397895 6.309918 6.309918
## Phalaenophana.pyramusalis Phaonia.serva Pidonia.ruficolis Pieris.napi Poanes.hobomok
## 6.309918 5.211306 4.112694 6.309918 6.309918
## Psithyrus.sp. Pyrophaena.rosarum Rhopalum.sp. Sarcophaga.nearctica Sphaerophoria.bifurcata
## 6.309918 6.309918 6.309918 5.616771 4.700480
## Sphaerophoria.longipilosa Sphaerophoria.sp. Sphecodes.sp. Syritta.pipiens Syritta.rectus
## 6.309918 5.211306 5.616771 6.309918 6.309918
## Temnostoma.alternans Temnostoma.balyras Temnostoma.barberi Thamnotettix.confinis Trachysida.aspera
## 6.309918 5.211306 5.616771 6.309918 6.309918
## Tropidia.quadrata Tychius.stephensi Vespula.arenaria Xylota.bigelowi Xylota.hinei
## 4.923624 5.211306 6.309918 4.007333 5.211306
## Xylota.sp. Zeraea.americana
## 6.309918 5.616771plot(degree(pollination)[V(pollination_bin)$type],dfun(t(mat_plantpol))$dprime,xlab="degrees",ylab="d'")
In the case of unipartite directed networks / food webs, dfun can also be used.
foodweb_spe.res<-dfun(as_adj(foodweb,sparse=FALSE))$dprime #as resource
foodweb_spe.con<-dfun(t(as_adj(foodweb,sparse=FALSE)))$dprime #as consumer
interm_species_list<-intersect(names(foodweb_spe.res),names(foodweb_spe.con))
plot(foodweb_spe.res[which(names(foodweb_spe.res)%in%interm_species_list)],foodweb_spe.con[which(names(foodweb_spe.con)%in%interm_species_list)],xlab="d' as prey",ylab="d' as consumer")
Network specialization
Blüthgen et al. (2006) also introduced a network-level metric of specialization which is based on Shannon’s entropy. The metric \(H_2\) is simply the entropy of interaction weights and the relative measure \(H_2^\prime\) is its opposite, scaled between 0 and 1.
H2fun(mat_plantpol)## H2 H2min H2max H2uncorr
## 0.5780733 3.4320000 5.0600000 4.1190000Modularity
Modularity, as defined by Newman, is given by the following formula: \[\begin{equation} Q = \frac{1}{A} \sum_{i,j} \left[ a_{ij} - \frac{d_i d_j}{A}\right]\delta_{ij} \end{equation}\] where \(A\) is the total number of connections in the network (i.e. \(A = \sum_i \sum_j a_{ij}\)), \(d_i\) is the degree of species \(i\) (i.e. \(d_i = \sum_j a_{ij}\)) and the \(\delta_{ij}\) are dummy variables indicating whether species/nodes \(i\) and \(j\) are assumed to belong to the same module/community/cluster. In layman’s terms, modularity is a quantity that measures how ‘’intense’’ interactions are within modules, i.e. compares \(a_{ij}\) to its expectation given the degrees of \(i\) and \(j\), but only taking into account edges within modules. The principle of modularity search is to look for the partitioning of nodes into communities/modules that maximizes the associated modularity score. The problem of modularity search is quite complex (NP-complete problem) and thus can be solved with different algorithms that have their pros and cons.
Modularity, in its classic form, only works for undirected networks. Three algorithms that are often used in ecology are the edge-betweenness algorithm (EB), the leading-eigenvector one (LE) and the Louvain algorithm (ML). Another useful one (and also more efficient) is simulated annealing, as implementing in rnetcarto.
foodweb_EB.mod<-cluster_edge_betweenness(undirected_foodweb)
foodweb_LE.mod<-cluster_leading_eigen(undirected_foodweb)
foodweb_ML.mod<-cluster_louvain(undirected_foodweb)
foodweb_SA.mod<-cluster_netcarto(undirected_foodweb)Which algorithm to compute modularity?
Because there are tons of algorithms already developed to find clusters/communities/modules in networks, it is easy to get lost… Yet two papers (at least) have made our life easier by comparing all these algorithms, bot for unipartite networks (Yang et al. 2016) and bipartite networks (Leger et al. 2015).
Plotting modularity
Here is a comparison of modules as obtained through the four algorithms used before.
par(mfrow=c(2,2))
plot(foodweb_EB.mod,foodweb,layout = layout_as_food_web(foodweb))
plot(foodweb_LE.mod,foodweb,layout = layout_as_food_web(foodweb))
plot(foodweb_ML.mod,foodweb,layout = layout_as_food_web(foodweb))
plot(foodweb_SA.mod,foodweb,layout = layout_as_food_web(foodweb))
It is also possible to plot the modules using the ‘’matrix’’ visualization of networks.
plotMyMatrix(as_adj(undirected_foodweb,sparse=FALSE),clustering=list("row"=foodweb_SA.mod$membership,"col"=foodweb_SA.mod$membership))
Comparing modules from the different algorithms
For the sake of visualizing whether the modules obtained with the different algorithms correspond, we can make use of alluvial plots.
make_alluvial_2(foodweb_ML.mod$membership,foodweb_LE.mod$membership,"Louvain","Leading_eigen")
make_alluvial_2(foodweb_ML.mod$membership,foodweb_SA.mod$membership,"Louvain","Simulated annealing")
Modularity on bipartite network
Modularity work as easily on bipartite networks as on unipartite ones. However, the representation of the modules is more useful when done on the bi-adjacency matrix.
pollination_SA.mod<-cluster_netcarto(pollination_bin)
plotMyMatrix(mat_plantpol_bin,clustering=list("row"=pollination_SA.mod$membership[!V(pollination_bin)$type],"col"=pollination_SA.mod$membership[V(pollination_bin)$type]))
Spectral clustering
A convenient alternative algorithm to find clusters within a network is spectral clustering. The idea of spectral clustering is to make use of the Laplacian matrix \(L\) of the graph, given by: \[L = D - A\] where \(D\) is the diagonal matrix of nodes’ degrees and \(A\) is the adjacency matrix.
The Laplacian matrix has interesting properties, the best of all being that the number of zeros among its eigenvalues yield the number of components of the graph (one 0 means there is a single connected component, two zeros means the network is split in two disconnected components, etc.). The eigenvectors associated the zero eigenvalues can inform on the membership of nodes to the different components.
The heuristics of the spectral clustering is to use the eigenvectors associated with the set of small non-zero eigenvalues of \(L\) to deduce the ‘’almost components’’ of the graph. To do so, one chooses a number of no-zero eigenvalues to keep and then the algorithm uses the associated eigenvector to find the groups through a K-means algorithm on the space obtained with these eigenvectors.
The function laplacian_spectral_gap computes the ‘’spectral gap’’ plot and tries to find the optimum number of clusters.
foodweb_SG<-laplacian_spectral_gap(undirected_foodweb)
foodweb_SG$optim_n## [1] 3Once we know how many clusters we are looking for, we can compute the spectral clustering.
foodweb_SC<-spectral_clustering(undirected_foodweb,foodweb_SG$optim_n)
plotMyMatrix(as_adj(undirected_foodweb,sparse=FALSE),clustering=list("row"=foodweb_SC,"col"=foodweb_SC))
modularity(undirected_foodweb,foodweb_SC)## [1] 0.3296616In passing, we note that the modularity of the spectral clustering is not so bad.
Comparison between modules and spectral clusters
With an alluvial plot it is possible to “see” how congruent the modules obtained by the Louvain algorithm are with the clusters identified by the spectral clustering method.
make_alluvial_2(foodweb_SA.mod$membership,foodweb_SC,"Simulated annealing","Spectral clustering")
Nestedness
The idea of nestedness is to measure the tendency for specialists to only interact with a subsample of the interactors of generalists. In other words, a network would be said “non-nested” only when specialists would interact more likely with specialists than with generalists.
Several indices have been developed over the years. Navigating their complexity is the purpose of the synthesis written by Podani and Schmera (2012).
The package bipartite proposes some methods to do this (partly borrowed from the package vegan), in the function nested.
pollination_nestdness<-nested(mat_plantpol, method = "ALL", rescale=TRUE, normalised=TRUE)Another option is to use the various functions directly from vegan:
nestednodf(mat_plantpol)## N columns : 30.72672
## N rows : 35.30757
## NODF : 30.78467
## Matrix fill: 0.1364379nestednodf(mat_plantpol,weighted=TRUE)## N columns : 13.49511
## N rows : 13.52354
## NODF : 13.49547
## Matrix fill: 0.1364379nestednodf(mat_plantpol,weighted=TRUE,wbinary=T)## N columns : 25.53355
## N rows : 25.85196
## NODF : 25.53758
## Matrix fill: 0.1364379nestednodf(mat_plantpol_bin)## N columns : 30.72672
## N rows : 35.30757
## NODF : 30.78467
## Matrix fill: 0.1364379nestednodf(mat_plantpol_bin,weighted=TRUE)## N columns : 0
## N rows : 0
## NODF : 0
## Matrix fill: 0.1364379nestednodf(mat_plantpol_bin,weighted=TRUE,wbinary=T)## N columns : 30.72672
## N rows : 35.30757
## NODF : 30.78467
## Matrix fill: 0.1364379Recent research seems to indicate that the concept of nestedness is not really fruitful for a variety of reasons: - it seems to be very strongly linked to connectance; - the diversity of nestedness indices tend to obfuscate its meaning (results reported show more about the index than about the concept); - these different indices can yield very different results when applied to the same dataset; - in the classic case of power law distributions of degrees, the nestedness scores tend to be completely detrmined by the exponent of the power law; - finally, when putting the common datasets in the right matrix ensemble (i.e. using correct null models to sample similar matrices), reported nestedness values seem to be completely expected.
Robustness to secondary extinctions
The general idea is to assume that when a species is removed from a network, the species interacting with it could go extinct (secondary extinction). Robustness analysis then consists in assessing how many species are lost when R species are removed. The most used indicator of robustness is the R50 which is the percentage of species on needs to remove from a web to lose 50% of its species.
In the paper of J. Dunne and collaborators, they used examples based on a variety of datasets. The basic rule then is that a species go extinct when it loses all its prey items (except itself if it is cannibalistic). One can obtain different scenarios based on how removed species are chosen (intentionality: how strongly the species’ degree changes its probability of being drawn randomly as the next one to be removed).
We can first try a robustness analysis on a virtual dataset, using FWebs (or the functions in crutch_functions.R). We will use the “niche model” to generate a virtual food web.
The plot yields the value of the R50 as a function of intentionality, which is a proxy for how strongly species removals target species based on their degree.
niche<-niche_matrix(0.2,200)
m<-niche$matrix
netw<-graph_from_adjacency_matrix(m,mode="directed")
V(netw)$name<-paste0("x",as.character(1:200))
i_index <- seq(from = 0, to = 0.95, by =0.05)
prob_exp<-exponent.removal(degree(netw,mode="out"), i_index)
netw_iterate_deg<-iterate(fw_to_attack=netw, prob_exp, alpha1=50, iter=10, i_index, plot = TRUE)## Iteration 1## Iteration 2## Iteration 3## Iteration 4## Iteration 5## Iteration 6## Iteration 7## Iteration 8## Iteration 9## Iteration 10
We can also try it on the real food web data:
prob_exp_FW<-exponent.removal(degree(foodweb,mode="out"), i_index)
fw_iterate_deg<-iterate(fw_to_attack=foodweb, prob_exp_FW, alpha1=50, iter=20, i_index, plot = TRUE)## Iteration 1## Iteration 2## Iteration 3## Iteration 4## Iteration 5## Iteration 6## Iteration 7## Iteration 8## Iteration 9## Iteration 10## Iteration 11## Iteration 12## Iteration 13## Iteration 14## Iteration 15## Iteration 16## Iteration 17## Iteration 18## Iteration 19## Iteration 20
prob_exp<-exponent.removal(trophic_levels(netw), i_index)
netw_iterate_TL<-iterate(fw_to_attack=netw, prob_exp, alpha1=50, iter=20, i_index, plot = TRUE)## Iteration 1## Iteration 2## Iteration 3## Iteration 4## Iteration 5## Iteration 6## Iteration 7## Iteration 8## Iteration 9## Iteration 10## Iteration 11## Iteration 12## Iteration 13## Iteration 14## Iteration 15## Iteration 16## Iteration 17## Iteration 18## Iteration 19## Iteration 20
foodweb_TL<-as.data.frame(TrophicLevels(as_Community(foodweb)))
prob_exp_FW<-exponent.removal(foodweb_TL$LongestTL, i_index)
fw_iterate_TL<-iterate(fw_to_attack=foodweb, prob_exp_FW, alpha1=50, iter=20, i_index, plot = TRUE)## Iteration 1## Iteration 2## Iteration 3## Iteration 4## Iteration 5## Iteration 6## Iteration 7## Iteration 8## Iteration 9## Iteration 10## Iteration 11## Iteration 12## Iteration 13## Iteration 14## Iteration 15## Iteration 16## Iteration 17## Iteration 18## Iteration 19## Iteration 20
Beta diversity of networks
Ohlmann et al. developed a measure of the “beta diversity” of nodes and links among networks that share some nodes. For nodes, this amounts to the classic beta diversity of communities. For links, two measures are possible, based on the following formula: \[L_{ij} = \pi_{ij}p_ip_j\] where - \(L_{ij}\) is the relative weight of the link between species \(i\) and \(j\), - \(p_i\) is the relative abundance of species \(i^\), - \(\pi_{ij}\) is the probability of interaction between species \(i\) and \(j\)j. The beta diversity of the \(p\)’s corresponds to the beta diversity of nodes (community diversity), while the beta diversity of either \(L_{ij}\) or \(\pi_{ij}\) provides a measure of link diversity among networks.
The functions are implement in package econetwork. We test it on four networks of the same 60 species.
gList <- c()
for(i in 1:4){
graphLocal <- sample_gnp(60, 0.1, directed=TRUE)
V(graphLocal)$name <- as.character(1:60)
gList <- c(gList, list(graphLocal))
}
names(gList) <- c("A","B","C","D")
disPairwise(gList, type='P')## | | | 0% | |======================= | 17% | |============================================== | 33% | |===================================================================== | 50% | |============================================================================================ | 67% | |=================================================================================================================== | 83% | |==========================================================================================================================================| 100%## A B C
## B 0
## C 0 0
## D 0 0 0disPairwise(gList, type='L')## | | | 0% | |======================= | 17% | |============================================== | 33% | |===================================================================== | 50% | |============================================================================================ | 67% | |=================================================================================================================== | 83% | |==========================================================================================================================================| 100%## A B C
## B 0.8912001
## C 0.8771519 0.9100341
## D 0.8768809 0.8770819 0.8930309disPairwise(gList, type='Pi')## | | | 0% | |======================= | 17% | |============================================== | 33% | |===================================================================== | 50% | |============================================================================================ | 67% | |=================================================================================================================== | 83% | |==========================================================================================================================================| 100%## A B C
## B 0.8912001
## C 0.8771519 0.9100341
## D 0.8768809 0.8770819 0.8930309As expected, the dissimilarity/distance between node sets is 0, but the distance among edge sets is greater than 0.
Network models
Statistics on networks are hard to test because everything depends on everything else. In order to test statistics, we need to construct models. These are generally classified in two types: - null models, which correspond to algorithmic models shuffling the initial network to obtain a series of networks that obey some constraints (e.g. same sequence of degrees of the nodes) - generative models, which provide a network based on a few parameters and some working hypotheses (e.g. all connections are equally likely and their probability is known)
Among generative models, some (but not all) are probabilistic (i.e. can be fitted), which clearly opens the door to the elaboration of statistical tests in the classic sense (testing whether the working hypothesis of the model holds).
The other generative models (e.g. preferential attachment) and all null models can be used to re-simulate the network in order to obtain the p-value corresponding to the occurrence of the statistic value under the hypothesis represented by the generative or null model.
How to make a food web form an undirected unipartite graph?
Before tackling models per se, one important step is to be able to convert any undirected unipartite network into a directed network that can be thought of as a food web. The main difficulties are that randomly orienting the links is likely - to cause cycles (which in principle do not exist in food webs), and - not to have a few species that feed on no one else (the “basal” level).
The make_food_web_from_undirected function tries to cope with this two difficulties. It asks for an undirected graph and a number of basal species. Then it proceeds in two steps:
- it colours the graph to find one color that has at least as many nodes as the basal set, and then it assigns the “basal species” into this colour;
- it orients the graph using a breadth-first search (see the following figure).
The configuration model
One of the most useful and simplest null model is the configuration model. Under this model, edges are randomly reattached to other nodes such that all nodes keep their degrees. This behaves effectively as if all nodes were attached to half-edges which were randomly paired until the resulting network yields no self-link and no double link.
Generate configurations for an undirected unipartite network
Here is an example to generate configurations starting with a random graph (Erdos-Rényi model).
net<-sample_gnp(50,0.2, directed=FALSE)
sample.config.undirected<-lapply(1:100,function(x) sample_degseq(degree(net), method = "vl"))
length(sample.config.undirected)## [1] 100The function config_VL (in the source file) does that automatically from the input network or from its adjacency matrix.
alt.config<-config_VL(net,100)## [1] "unipartite graph detected"plot(alt.config[[42]])
Example of use of unipartite configuration model
Here, we exemplify the use of the configuration model with a simple issue. We want to assess whether the modularity of the largest component of the food web data is significant (in the sense: not expected from the configuration model). In this example, food webs are considered in their undirected form.
one_comp_foodweb<-as.undirected(foodweb)
one_comp_foodweb_LE.mod<-cluster_leading_eigen(one_comp_foodweb)
one_comp_foodweb_LE.mod$mod## [1] 0.3630662foodweb_configs<-config_VL(one_comp_foodweb,1000)## [1] "unipartite graph detected"mods_configs<-sapply(1:1000, function(x) cluster_leading_eigen(foodweb_configs[[x]])$mod)The function p.val (from the source file) computes the p-value based on a collection of networks (here generated with the configuration model) and their modularities. Here, we specify we want a “larger than” test.
p.val(test_val=one_comp_foodweb_LE.mod$mod,test_collection=mods_configs,method="larger",label="modularity")
## [1] 0.000999001one_comp_foodweb_SA.mod<-cluster_netcarto(one_comp_foodweb)
mods_configs_SA<-sapply(1:1000, function(x) cluster_netcarto(foodweb_configs[[x]])$mod)
p.val(test_val=one_comp_foodweb_SA.mod$mod,test_collection=mods_configs_SA,method="larger",label="modularity")
## [1] 0.000999001Randomization of directed networks
To randomize a directed network, we can make use of “sample_degseq” also, but we are not assured of uniform sampling anymore.
net<-sample_gnp(50,0.2, directed =FALSE)
net_directed<-make_food_web_from_undirected(net,5)
sample.config.directed<-lapply(1:100,function(x) sample_degseq(degree(net_directed,mode="out"), degree(net_directed,mode="in"), method = "fast.heur.simple"))We can plot the initial food web and one of its configurations side by side just to check that these are different.
par(mfrow=c(1,2))
plot(net_directed,layout=layout_as_food_web3(net_directed))
plot(sample.config.directed[[1]],layout=layout_as_food_web3(sample.config.directed[[1]]))
Configuration model for bipartite networks
In the case of bipartite networks, Strona et al. proposed the “curveball” algorithm which is efficient and ensures uniformity of the sampling of equivalent bipartite networks.
It is implemented in the package vegan with functions simulate and nullmodel (pay attention: there is also a nullmodel function in bipartite). The result is given as a “cube” of data (same matrix dimensions, but with depth for all simulated networks).
net<-sample_bipartite(40,60,"gnp",0.1)
net<-largest_component(net)
net_mat<-as.matrix(as_biadjacency_matrix(net))
sample.bip.config<-simulate(vegan::nullmodel(net_mat,"curveball"),nsim=1000)
dim(sample.bip.config)## [1] 40 60 1000This is also implented in the config_VL function.
alt.config<-config_VL(net,1000)## [1] "bipartite graph detected"Example of use of bipartite configuration model
Here, we exemplify the use of the bipartite configuration model to ascertain the significance of the nestedness of the bipartite network “net_mat” (created earlier usinf the bipartite ER model).
net_nestedness<-nested(net_mat, method = "NODF2")
nestedness_configs<-sapply(1:1000, function(x) nested(sample.bip.config[,,x], method = "NODF2"))
p.val(net_nestedness,nestedness_configs,method="two-sided",label="nestedness")
## [1] 0.3476523Generative models
Now we turn to generative models, i.e. those models that yield a network based on summary parameters (not models that turn out a network from a network). In the following, we will focus on a few of those, mainly Erdos-Rényi, preferential attachment, models with invariants, models with prescribed degree sequence and the expected degree distribution model.
Food web generation using Erdos-Rényi
The ER algorithm assumes that all edges are equally likely. To construct a food web, one can use the food web construction from an undirected algorithm presented above.
net<-make_food_web_from_undirected(sample_gnp(50,0.2, "gnp",directed =TRUE),basal = 10)## [1] "not enough nodes in the largest color to make basal species, "
## [1] "the number of basal species is set to"
## [1] 9V(net)$name<-sapply(1:50,function(x) paste0("sp_",x))
plot(net,layout=layout_as_food_web2(net))
It is possible to test whether the distribution of degrees follows a binomial (the expectation), but it looks significant a bit too easily…
ks.test(degree(net),"pbinom",length(V(net))-1,mean(degree(net))/(length(V(net))-1))## Warning in ks.test.default(degree(net), "pbinom", length(V(net)) - 1, mean(degree(net))/(length(V(net)) - : ties should not be present for the
## one-sample Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: degree(net)
## D = 0.17473, p-value = 0.09443
## alternative hypothesis: two-sidedFood web generation using preferential attachment
Undirected networks with PA can be obtained using the function sample_pa. The argument power controls how fat the tail of the distribution is. As above, the undirected network is generated and then turned into a food web.
net<-make_food_web_from_undirected(sample_pa(500, power = 1, directed = FALSE, out.pref=TRUE),basal = 50)
plot(net,layout=layout_as_food_web2(net))
It is also possible to test whether the distribution of degrees follows a power law distribution (“ppldis”).
net_dist<-displ$new(degree(net))
net_dist$setXmin(estimate_xmin(displ$new(degree(net)))$xmin)
net_pars<-estimate_pars(net_dist)
ks.test(degree(as.undirected(net)),"ppldis",xmin = estimate_xmin(displ$new(degree(net)))$xmin, alpha = net_pars$pars)## Warning in ks.test.default(degree(as.undirected(net)), "ppldis", xmin = estimate_xmin(displ$new(degree(net)))$xmin, : ties should not be present
## for the one-sample Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: degree(as.undirected(net))
## D = 0.596, p-value < 2.2e-16
## alternative hypothesis: two-sidedFood web models with invariants: the cascade model
In the 80-90’s, one important goal of interaction network research was to find invariant properties, i.e. statistics or metrics that did not vary among food webs. With this goal in mind, two models were proposed.
First, the cascade model of Cohen & Briand proposed that predators only eat a rather fixed number of prey species. This corresponds to predators eating prey that are are smaller than them with probability \(c/S\) (\(c\) being a fixed parameter), so that connectance decreases as \(1/S\).
net<-cascade_matrix(3,50)
sum(net)/(dim(net)[1]*(dim(net)[1]-1))## [1] 0.03183673graph_cascade<-graph_from_adjacency_matrix(net)
V(graph_cascade)$name<-sapply(1:50,function(x) paste0("sp_",x))
plot(graph_cascade,layout=layout_as_food_web3(graph_cascade))
net<-cascade_matrix(3,200)
sum(net)/(dim(net)[1]*(dim(net)[1]-1))## [1] 0.007035176Food web models with invariants: the niche model
The second model that was developed to find invariants is the niche model of Williams and Martinez. In this model, all species are supposed to be defined by their size and to eat all species with size contained within an interval centred around a size lower than the predator’s. With the right parameter values, one can use this model to generate a food web that has a given connectance, whatever the number of species in the food web.
net<-niche_matrix(0.1,50)
sum(net$matrix)/(dim(net$matrix)[1]*(dim(net$matrix)[1]-1))## [1] 0.1081633graph_niche<-graph_from_adjacency_matrix(net$matrix)
V(graph_niche)$name<-sapply(1:50,function(x) paste0("sp_",x))
plot(graph_niche,layout=layout_as_food_web3(graph_niche))## Warning in min(vec[!is.infinite(vec)]): aucun argument trouvé pour min ; Inf est renvoyé
## Warning in min(vec[!is.infinite(vec)]): aucun argument trouvé pour min ; Inf est renvoyé
## Warning in min(vec[!is.infinite(vec)]): aucun argument trouvé pour min ; Inf est renvoyé
## Warning in min(vec[!is.infinite(vec)]): aucun argument trouvé pour min ; Inf est renvoyé
## Warning in min(vec[!is.infinite(vec)]): aucun argument trouvé pour min ; Inf est renvoyé## Warning in max(vec[!is.infinite(vec)]): aucun argument pour max ; -Inf est renvoyé
## Warning in max(vec[!is.infinite(vec)]): aucun argument pour max ; -Inf est renvoyé
## Warning in max(vec[!is.infinite(vec)]): aucun argument pour max ; -Inf est renvoyé
## Warning in max(vec[!is.infinite(vec)]): aucun argument pour max ; -Inf est renvoyé
## Warning in max(vec[!is.infinite(vec)]): aucun argument pour max ; -Inf est renvoyé
net<-niche_matrix(0.1,200)
sum(net$matrix)/(dim(net$matrix)[1]*(dim(net$matrix)[1]-1))## [1] 0.101005Food webs with prescribed degree sequence
Here, the idea is not to generate a degree distribution in general, but rather to use a degree sequence that is “typical” of the wanted degree distribution. To do that, from the number of nodes/species and the distribution, one can extract the typical degrees of the nodes using the quantiles of the distribution. Then, using the Viger-Latapy’s algorithm, it is possible to generate a random network, and then make it a food web.
Here, an example with the Poisson distribution.
net<-generate_DM_model(50,0.1,quant_fun="qpois", lambda=5)## Warning: `as.directed()` was deprecated in igraph 2.1.0.
## i Please use `as_directed()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## Warning: `subgraph.edges()` was deprecated in igraph 2.1.0.
## i Please use `subgraph_from_edges()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.plot(net,layout=layout_as_food_web3(net))
And now an example with the negative binomial distribution.
net<-generate_DM_model(50,0.1,quant_fun="qnbinom", mu = 4, size=50)
plot(net,layout=layout_as_food_web3(net))
It is also possible to check whether the degrees follow the desired distribution.
net<-generate_DM_model(450,0.1,quant_fun="qpois", lambda=7)
hist(degree(net),breaks=0:max(degree(net)),main="")
ks.test(degree(net),"ppois",lambda = mean(degree(net)))## Warning in ks.test.default(degree(net), "ppois", lambda = mean(degree(net))): ties should not be present for the one-sample Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: degree(net)
## D = 0.14916, p-value = 4.022e-09
## alternative hypothesis: two-sidedFood webs following the expected degree distribution (EDD) model
The EDD model assumes that each node has a latent parameter (its expected degree) which is used to assess whether edges exist or not (multiplicative model). These latent parameters are assigned randomly among nodes following a given distribution (actually, using an increasing function of a uniform variable).
The EDD model can generate networks with more or less heterogeneity in expected degrees (the “lambda” parameter of the function Simul_EDD).
net_EDD<-Simul_EDD(50,0.2,2)
graph_EDD<-graph_from_adjacency_matrix(net_EDD$net,mode="undirected")
graph_EDD<-largest_component(graph_EDD)
plotMyMatrix(net_EDD$net)
foodweb_EDD<-make_food_web_from_undirected(graph_EDD,5)
plot(foodweb_EDD,layout=layout_as_food_web3(foodweb_EDD))
Bipartite network generation using Erdos-Rényi
Here we simply use the sample_bipartite function of igraph to generate a random bipartite network in which all edges between level 1 and level 2 have the same probability of occurring.
net<-sample_bipartite(25,25,"gnp",0.1)
plot(net,layout=layout_as_bipartite)
plotweb(as_biadjacency_matrix(net))
Bipartite networks using the bipartite EDD (BEDD) model
The BEDD model is the analogue of the EDD model, but for bipartite networks. Each node receives a random parameter which determines the probability of edge occurrences. The “hyper-parameters” controlling the heterogeneity of node degrees can be adjusted separately for the two levels of the network.
net_bedd<-SimulB_EDD(25, 25, 0.1, 2, 2)
graph_bedd<-graph_from_biadjacency_matrix(net_bedd$net)
plot(graph_bedd,layout=layout_as_bipartite)
plotweb(net_bedd$net)
plotMyMatrix(net_bedd$net)
Randomization through the BEDD model: nestedness, BEDD vs. configuration
The BEDD model can be used to randomize a network and thus test the significance of some properties such as nestedness.
sample.bip.EDD<-lapply(1:1000, function(x) randomize.BEDD(net_mat))
nestedness_EDD<-sapply(1:1000, function(x) nested(sample.bip.EDD[[x]], method = "NODF2"))
p.val(net_nestedness,nestedness_EDD,method="two-sided",label="nestedness")
## [1] 0.1238761It is enlightening to compare the density of simulated values of nestedness according to two different null models, here the configuration model (in blue) and the BEDD model (in green).
plot(density(nestedness_configs),xlim=c(min(nestedness_configs,nestedness_EDD),max(nestedness_configs,nestedness_EDD)),xlab="",main="nestedness distributions",col="blue")
lines(density(nestedness_EDD),col="darkgreen")
It is apparent that the BEDD allows for a wider variability in nestedness in simulated networks than the configuration model. So p-values obtained from each model have a different interpretation…
Randomization using the EDD model: modularity, EDD vs. configuration
Another example of use of the EDD model, this time to assess the significance of modularity.
foodweb_EDD<-lapply(1:1000, function(x) randomize.EDD(one_comp_foodweb))
mods_EDD<-sapply(1:1000, function(x) cluster_leading_eigen(as.undirected(foodweb_EDD[[x]]))$mod)
p.val(one_comp_foodweb_LE.mod$mod,mods_EDD,method="larger",label="modularity")
## [1] 0.000999001Again, we can compare the distribution of modularity scores generated by the EDD model (green) and the configuration model (blue).
plot(density(mods_configs),xlim=c(min(mods_configs,mods_EDD),max(mods_configs,mods_EDD)),xlab="",main="modularity distributions",col="blue")
lines(density(mods_EDD),col="darkgreen")
This time, it is not only the extent of the distribution, but also its mean and mode that are completely shifted by the change in underlying null model.
References
Barrett, S. C. & Helenurm, K. (1987) The reproductive biology of boreal forest herbs. I. Breeding systems and pollination. Canadian Journal of Botany, 65, 2036-2046.
Blüthgen, N., Menzel, F. & Blüthgen, N. (2006) Measuring specialization in species interaction networks. BMC ecology, 6, 9.
Clauset, A., Shalizi, C. R. & Newman, M. E. J. (2009) Power-law distributions in empirical data. SIAM Review, 51, 661-703.
Cruz-Escalona, V. H., Arreguín-Sánchez, F. & Zetina-Rejón, M. (2007) Analysis of the ecosystem structure of Laguna Alvarado, western Gulf of Mexico, by means of a mass balance model. Estuarine, Coastal and Shelf Science, 72, 155-167.
Dunne, J. A., Williams, R. J. & Martinez, N. D. (2002) Food-web structure and network theory: The role of connectance and size. Proceedings of the National Academy of Sciences, 99, 12917-12922.
Dunne, J. A., Williams, R. J. & Martinez, N. D. (2002) Network structure and biodiversity loss in food webs: robustness increases with connectance. Ecology Letters, 5, 558-567.
James, A., Pitchford, J. W. & Plank, M. J. (2012) Disentangling nestedness from models of ecological complexity. Nature, 487, 227-230.
Leger, J.-B., Daudin, J.-J. & Vacher, C. (2015) Clustering methods differ in their ability to detect patterns in ecological networks. Methods in Ecology and Evolution, 6, 474-481.
Lewinsohn, T. M., Prado, P. I., Jordano, P., Bascompte, J. & Olesen, J. M. (2006) Structure in plant-animal interaction assemblages. Oikos, 113, 174-184.
MacKay, R. S., Johnson, S. & Sansom, B. (2020) How directed is a directed network? Royal Society Open Science, 7, 201138.
Massol, F., Dubart, M., Calcagno, V., Cazelles, K., Jacquet, C., Kéfi, S. & Gravel, D. (2017) Island biogeography of food webs. Advances in Ecological Research vol. 56 - Networks of Invasion: A Synthesis of Concepts (eds D. A. Bohan, A. J. Dumbrell & F. Massol), pp. 183-262. Academic Press.
Newman, M. E. J. (2006) Modularity and community structure in networks. Proceedings of the National Academy of Sciences, 103, 8577-8582.
Pocock, M. J. O., Evans, D. M. & Memmott, J. (2012) The robustness and restoration of a network of ecological networks. Science, 335, 973-977.
Payrató-Borràs, C., Hernández, L. & Moreno, Y. (2019) Breaking the Spell of Nestedness: The Entropic Origin of Nestedness in Mutualistic Systems. Physical Review X, 9, 031024.
Podani, J. & Schmera, D. (2012) A comparative evaluation of pairwise nestedness measures. Ecography, 35, 889-900.
Strona, G., Nappo, D., Boccacci, F., Fattorini, S. & San-Miguel-Ayanz, J. (2014) A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals. Nature Communications, 5.
Stumpf, M. P. H. & Porter, M. A. (2012) Critical truths about power laws. Science, 335, 665-666.
von Luxburg, U. (2007) A tutorial on spectral clustering. Statistics and Computing, 17, 395-416.
Williams, R. J. & Martinez, N. D. (2000) Simple rules yield complex food webs. Nature, 404, 180-183.
Yang, Z., Algesheimer, R. & Tessone, C. J. (2016) A Comparative Analysis of Community Detection Algorithms on Artificial Networks. Scientific Reports, 6, 30750.