Pakiet R w analizie układów złożonych

LABORATORIUM 8

Sieci złożone są pewnym przedstawieniem powiązań pomiędzy wyodrębnionymi jednostkami, przy czym pojęcie jednostki jest bardzo szerokie - może być to człowiek, zwierzę, komputer, przystanek autobusowy czy słowo. Podejście sieci opiera się na teorii grafów: w tym ujęciu graf \(\mathcal{G}\) to zbiór węzłów (wierzchołków, nodes, verices) \(\mathbb{V}\), połączonych krawędziami (edges, links) zawartymi w zbiorze \(\mathbb{E}\). W przypadku pakietu R tworzenie o obsługa sieci może być oparta o np. bibliotekę igraph.

Tworzenie sieci

Rozpoczynamy od stworzenia pustego grafu (make_empty_graph()), do którego następnie dodamy pięć wierzchołków oraz trzy krawędzie. Za każdym razem możliwe jest narysowanie grafu przeciążoną funkcją plot(), jak również wypisanie atrybutów obiektu poprzez wywołanie komendy print() lub po prostu wypisując zawartośc zmiennej. Za pomocą operatora + można dodawać zarówno węzły jak i krawędzie.

## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## IGRAPH b3592ca D--- 0 0 -- 
## + edges from b3592ca:
## IGRAPH df624ff DN-- 5 0 -- 
## + attr: name (v/n)
## + edges from df624ff (vertex names):

## IGRAPH f685acc DN-- 5 3 -- 
## + attr: name (v/n)
## + edges from f685acc (vertex names):
## [1] 1->2 3->4 3->5

Najbardziej użytecznymi i wykorzystywanymi funkcjami są V() oraz E() - umożliwiają one wypisanie, odpowiednio, wszystkich węzłów i krawędzi w sieci. Operator - można zastosować do usuwania węzłów i krawędzi.

## + 5/5 vertices, named, from f685acc:
## [1] 1 2 3 4 5
## + 3/3 edges from f685acc (vertex names):
## [1] 1->2 3->4 3->5
## IGRAPH 4256555 DN-- 5 2 -- 
## + attr: name (v/n)
## + edges from 4256555 (vertex names):
## [1] 1->2 3->4
## IGRAPH e7ff3e6 DN-- 3 1 -- 
## + attr: name (v/n)
## + edge from e7ff3e6 (vertex names):
## [1] 1->2

Sieć można także utworzyć podając listę krawędzi

## IGRAPH 2fee21d DN-- 4 3 -- 
## + attr: name (v/c)
## + edges from 2fee21d (vertex names):
## [1] EUR->PLN PLN->USD PLN->GBP
## + 4/4 vertices, named, from 2fee21d:
## [1] EUR PLN USD GBP
## + 3/3 edges from 2fee21d (vertex names):
## [1] EUR->PLN PLN->USD PLN->GBP

Warto tu od razu zaznaczyć, że funkcja plot() ma sporą listę opcji, związanych z formatowaniem węzłów i krawędzi - można ją zobaczyć wywołując helpa dla igraph.plotting. Efekty dla poprzedniego grafu mogą być np. takie:

Oprócz grafu pustego możliwe jest także wygenerowanie innych standardowych konstrukcji: grafu pełnego (make_full_graph()), pierścienia (make_ring()), siatki lub kostki (make_lattice(dimvector = c(n1, n2))), gwiazdy (make_star()) czy drzewa (make_tree()). Dla podanych konstruktorów możliwe jest określenie czy graf ma być skierowany czy nie (directed=TRUE/FALSE), w przypadku drzew i gwiazdy parametrem in/out określamy, w którą stronę ma odbywać się skierowanie.

Modele sieci

Model grafów przypadkowych Erdősa-Rényi

Jednym z podstawowych i dobrze znanych modeli sieci jest model grafow przypadkowych Erdősa-Rényi (skracany jako ER, lata 50-te). Jego idea jest bardzo prosta: każdą parę z \(N\) wierzchołków sieci łączymy z prawdopodobieństwem \(p\). Za tworzenie grafu ER odpowiedzialna jest funkcja erdos.renyi.game(), gdzie jako parametry podajemy liczbę wierzchołków oraz prawdopodobieństwo \(p\) (opcja "gnp") lub też po prostu liczbę krawędzi (opcja "gnm").

## IGRAPH b1f919b U--- 10 23 -- Erdos renyi (gnp) graph
## + attr: name (g/c), type (g/c), loops (g/l), p (g/n)
## + edges from b1f919b:
##  [1] 2-- 3 3-- 4 3-- 5 1-- 6 4-- 6 5-- 6 1-- 7 2-- 7 3-- 7 5-- 7 2-- 8 3-- 8
## [13] 4-- 8 5-- 8 6-- 8 1-- 9 2-- 9 4-- 9 1--10 2--10 5--10 7--10 8--10

## IGRAPH c84f257 U--- 10 10 -- Erdos renyi (gnm) graph
## + attr: name (g/c), type (g/c), loops (g/l), m (g/n)
## + edges from c84f257:
##  [1] 1-- 3 3-- 5 2-- 6 4-- 6 2-- 7 4-- 8 2-- 9 7-- 9 4--10 9--10

Model Wattsa-Strogatza

Istnieje możliwość płynnego przechodzenia pomiędzy periodyczną strukturą sieci (łańcuch, siatka, kostka) a grafem przypadkowym. Zapewnia to model Wattsa-Strogatza (watts.strogatz.game()) - startujemy od struktury periodycznej, a następnie z prawdopodobieństwem \(p\) dokonujemy tzw. przełączenia krawędzi poprzez wylosowanie nowego sąsiada przypadkowo wybranego węzła.

Model Barabasi'ego-Albert

Trzecim podstawowym modelem jest sieć rosnąca Barabasi'ego-Albert. W tym przypadku mamy do czynienia z układem dynamicznym: w każdym kroku czasowym do istniejących węzłów dołącza się nowy za pomocą \(m\) krawędzi, przy czym prawdopodobieństwo dołączenia się do jakiegoś węzła jest wprost proporcjonalne do liczby krawędzi, które węzeł już posiada. W pakiecie R sieć BA tworzymy za pomocą funkcji barabasi.game().

Własciwości sieci

Podstawową własnościa dowolnej sieci jest jej macierz sąsiedztwa \(A\) (adjacency matrix) - \(A_{ij}\), gdzie \(i\) oraz \(j\) są numerami węzłów, przyjmuje wartość 1 w przypadku, gdy pomiędzy wspomnianymi węzłami istnieje krawędź oraz 0 w przeciwnym wypadku. Za pomocą funkcji get.adjacency() można otrzymać tę macierz w formie "rzadkiej", przy wybraniu dodatkowo opcji sparse=FALSE otrzymujemy macierz standardową.

## IGRAPH f4d5b68 D--- 10 17 -- Barabasi graph
## + attr: name (g/c), power (g/n), m (g/n), zero.appeal (g/n), algorithm
## | (g/c)
## + edges from f4d5b68:
##  [1]  2->1  3->1  3->2  4->1  4->2  5->2  5->4  6->2  6->1  7->1  7->3  8->2
## [13]  8->1  9->4  9->7 10->7 10->2
## 10 x 10 sparse Matrix of class "dgCMatrix"
##                          
##  [1,] . . . . . . . . . .
##  [2,] 1 . . . . . . . . .
##  [3,] 1 1 . . . . . . . .
##  [4,] 1 1 . . . . . . . .
##  [5,] . 1 . 1 . . . . . .
##  [6,] 1 1 . . . . . . . .
##  [7,] 1 . 1 . . . . . . .
##  [8,] 1 1 . . . . . . . .
##  [9,] . . . 1 . . 1 . . .
## [10,] . 1 . . . . 1 . . .
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    0    0    0    0    0    0    0    0     0
##  [2,]    1    0    0    0    0    0    0    0    0     0
##  [3,]    1    1    0    0    0    0    0    0    0     0
##  [4,]    1    1    0    0    0    0    0    0    0     0
##  [5,]    0    1    0    1    0    0    0    0    0     0
##  [6,]    1    1    0    0    0    0    0    0    0     0
##  [7,]    1    0    1    0    0    0    0    0    0     0
##  [8,]    1    1    0    0    0    0    0    0    0     0
##  [9,]    0    0    0    1    0    0    1    0    0     0
## [10,]    0    1    0    0    0    0    1    0    0     0

Zsumowanie poszczególnych kolumn i wierszy odpowiada wyznaczeniu dla danej sieci stopni wszystkich wierzchołków, czyli liczby krawedzi (wchodzących i wychodzących) związanych z wierzchołkami (w R możemy się do tego posłużyć funkcją degree()). Rozkład prawdopodobieństwa stopni wierzchołków w znaczący sposób różni się dla poszczególnych sieci i dlatego jest jedną z najczęściej wykorzystywanych własności sieci. W przykładzie poniżej wygenerujemy nieskierowaną sieć BA dla \(N=100\) węzłów i \(m=2\), a potem dla takiej samej ilości węzłów i krawędzi (funkcje vcount() i ecount() zliczają, odpowiednio, liczbę węzłów i krawędzi grafu) stworzymy sieć ER. Dla obu wypiszemy stopnie wierzchołków, a następnie policzymy ich rozkłady.

##   [1] 12 10 13  8  9 14  7 14  6  6  5 12 13  3  4  3  5  9  4  3  8  4  7  4  4
##  [26]  3  4  5  7  2  2  6  3  2  7  3  3  2  2  4  3  2  5  3  3  2  4  2  6  3
##  [51]  2  3  2  4  3  2  2  2  3  2  5  2  3  2  2  2  4  4  2  2  2  2  3  2  2
##  [76]  2  4  2  2  4  2  3  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
##   [1]  4  3  4  0  3  3  6  4  5  6  7  3  3  4  3  5  4  2  4  3  8  5 10  7  4
##  [26]  1  4  3  5  3  5  1  2  5  4  5  3  3  2  6  3  5  3  8  3  0  3  4  3  3
##  [51]  2  3  5  4  7  5  3  4  7  1  4  2  4  4  5  5  1  1  5  6  3  3  4  1  2
##  [76]  7  8  6  3  3  3  3  2  4  2  4  8  4  3  2  3  5  5  8  3  2  5  5  2  6

Ponadto własciwości sieci dobrze oddaje średnia najkrótsza droga (average path length), będąca liczbą najmniejszych kroków, za pomocą których możemy się dostać od jednego węzła do drugiego, uśredniona po wszystkich parach węzłów w sieci. Na poniższym przykładzie widać jak drastycznie zmiejsza się jej wartość podczas zwiększania wartości parametru p w modelu Wattsa-Strogatza.

##        p          l
## 1  0e+00 125.375375
## 2  1e-05 125.375375
## 3  5e-05 125.375375
## 4  1e-04 125.375375
## 5  5e-04 101.633934
## 6  1e-03  86.094545
## 7  5e-03  29.468687
## 8  1e-02  19.852713
## 9  5e-02   8.694096
## 10 1e-01   6.874995
## 11 5e-01   5.253807
## 12 1e+00   5.067366
## Warning: Transformation introduced infinite values in continuous x-axis

Współczynnik gronowania

Współczynnik gronowania (klasteryzacji) jest miarą określającą stopnień tego, jak bardzo węzły w sieci się "skupiają". Formalnie dla węzła i można go zdefiniować jako prawdopodobieństwo, że dowolna para jego sąsiadów ma pomiędzy sobą połączenie. Prawdopodobieństwo jest równe ilości faktycznych połączeń pomiędzy sąsiadami \(E_i\) podzielonej przez maksymalną ilość takich połączeń \(E_{max}\)

\(C_i = \frac{E_i}{E_{max}} = \frac{2 E_i}{k_i(k_i-1)}\)

W R liczenie współczynnika gronowania realizowane jest przez funkcję transitivity(g, "local"). W przypadku, gdy chcemy wyznaczyć średni współczynnik gronowania w całej sieci, korzystamy z opcji "localaverage".

## [1] 0.00000000 0.03571429 0.16666667 0.16666667 0.00000000
## [1] 0.7222222

Współczynnik gronowanie jest drugą (obok średniej najkrótszej drogi) cechą wyróżniającą przejście od silnie periodycznych struktur (np. udekorowany pierścień) do grafów przypadkowych.

##        p          l           c
## 1  0e+00 125.375375 0.500000000
## 2  1e-05 125.375375 0.500000000
## 3  5e-05 125.375375 0.499633333
## 4  1e-04  94.925642 0.499633333
## 5  5e-04  97.314891 0.499633333
## 6  1e-03  94.368555 0.497200000
## 7  5e-03  27.324831 0.487833333
## 8  1e-02  20.106951 0.463100000
## 9  5e-02   9.310428 0.395133333
## 10 1e-01   7.045063 0.277158110
## 11 5e-01   5.230180 0.009000551
## 12 1e+00   5.109328 0.004920808
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis

Współczynnik gronowania pozwala także na wykazanie, czy w dana sieć ma charakter hierarchiczy - w takim przypadku będzie można zaobserwować skalowanie wartości współczynnika gronowania w funkcji stopnia wierzchołka.

## Warning: Transformation introduced infinite values in continuous y-axis

Wspólnoty

Odrębnym działem "sieciologii" jest poszukiwanie w strukturze sieci wspólnoty (zbiorowości, community). W najprostszy sposób można ją określić jako własność, która grupuje węzły ze względu na jakieś kryterium. Oznacza to, że algorytmy to wykrywania wspólnot starają się znaleźć w sieci pewne ich silnie połączone podgrafy za pomocą optymalizacji określonego kryterium. W efekcie otrzymujemy przypisanie węzła do numeru określającego community. Jedną z funkcji, wykorzystywaną w R do wykrywania wspólnot jest fastgreedy.community()

## IGRAPH clustering fast greedy, groups: 6, mod: 0.7
## + groups:
##   $`1`
##   [1]  3  6 12 13 24 25 26 27
##   
##   $`2`
##   [1]  4  8 16 17 32 33 34 35
##   
##   $`3`
##   [1]  7 14 15 28 29 30 31
##   
##   $`4`
##   + ... omitted several groups/vertices

## IGRAPH clustering fast greedy, groups: 5, mod: 0.37
## + groups:
##   $`1`
##   [1]  6  7 10 15 20 25 36
##   
##   $`2`
##    [1]  1  4  5 11 12 16 17 19 22 35 37 38
##   
##   $`3`
##   [1]  2  3 14 21 27 33 39
##   
##   $`4`
##   + ... omitted several groups/vertices

W wielu przypadkach oprócz samego istnienia linku, może być istotna jego waga - czyli siła oddziaływania pomiędzy węzłami \(i\) oraz \(j\). W poniższym przypadku wagi te zotsały nałożone losowo, ale w danych rzeczywistych mogą wynikać np. z ilości wystąpień połączenia pomiędzy węzłami. Aby można było zwizualizować wagi, które są przypisane do atrybutu weight funkcji E(g) należy je przekopiować do atrybutu E(g)$width (grubość krawędzi).

Wagi linków

##  [1] 2.4517986 0.9702205 3.6095231 2.6328141 2.8292734 0.4004281 1.1744146
##  [8] 3.2194118 2.9377264 1.4820081 0.2315733 0.5183020 5.1284077 6.8785059
## [15] 6.8263576 3.3739493 7.1504995 0.2846422 0.1423622 2.2515223 3.6311662
## [22] 1.1569477 0.3338656 0.9331507 2.0514359 1.0240020 1.5900451 1.5477286
## [29] 0.8429402 0.3654166 1.9161312 0.2452568 1.5265519 1.8337940 0.4195793
## [36] 0.7978760 4.0307303 6.2033995 0.1344158

Dla sieci ważonej znalezione wspólnoty będą wyglądać inaczej, niż dla przypadku nieważonego.

## IGRAPH clustering fast greedy, groups: 7, mod: 0.76
## + groups:
##   $`1`
##   [1]  5 10 11 20 21 22 23 40
##   
##   $`2`
##   [1]  8 16 17 32 33 34 35
##   
##   $`3`
##   [1]  3  6 13 26 27
##   
##   $`4`
##   + ... omitted several groups/vertices

## IGRAPH clustering fast greedy, groups: 6, mod: 0.7
## + groups:
##   $`1`
##   [1]  3  6 12 13 24 25 26 27
##   
##   $`2`
##   [1]  4  8 16 17 32 33 34 35
##   
##   $`3`
##   [1]  7 14 15 28 29 30 31
##   
##   $`4`
##   + ... omitted several groups/vertices

Analogiczna operacją do wyznaczania stopnia węzła jest dla sieci ważonej otrzymanie jego siły (funkcja graph.strength()), czyli sumy wag na wszystkich linkach węzła

\(s_i = \sum\limits_jw_{ij}\)

Rozkład siły \(p(s)\) jest podobny do rozkładu stopni wierzchołków \(p(k)\) lecz w określony sposób "rozmyty".

##   [1] 62.4180561 19.4268854 59.3780648  8.5399243 16.5462727 28.0770631
##   [7]  8.5382641 41.3889050 21.0984970  9.3258696 28.3535951 26.0288996
##  [13] 13.4052866  2.6597193  4.8269864  4.7921899  7.5163703  0.7929933
##  [19]  2.1897834  7.0793272  5.9525101 10.4623810 11.8117316  4.6814616
##  [25]  3.1570319  2.8813833  2.7606192  4.4973656  5.0106780  7.3815302
##  [31]  4.2846032  4.0520840 10.4642160  6.2862452 11.5033267  3.2407062
##  [37]  2.4226225  3.4084923  1.7259216  8.1168065 12.9327736  6.8316595
##  [43]  2.5609128  2.6737257  3.3265314  3.2075637  8.7128582  4.4329243
##  [49] 14.3962643  3.7254788  2.3294697  1.8166567  7.2718625 12.4080155
##  [55]  3.0235519  0.5687345  0.8976952  6.9527726  1.1283782  2.6591385
##  [61]  9.3853760  4.4932628  6.3848087  5.9135851  6.5610157  1.5112289
##  [67]  4.2488150  1.6271446 10.4144662  6.7056781  8.2341059  3.3395781
##  [73]  1.4351111  1.6993626  2.5194869  8.0375690  8.5472361  1.0672907
##  [79]  6.2043162  0.3742316  3.9587842  2.8327015 10.2694782  5.8883986
##  [85] 12.1019077  4.8603331  5.5028110  3.0120892  2.4128934  2.1191618
##  [91]  5.0017046  4.3149835  6.6474233  4.9493975  6.5972020  1.8506928
##  [97]  3.4806744  2.8817294  4.3632639  0.9548689
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.