Introdução aos Grafos Aleatórios

Adrian Hinojosa
Departamento de Estatística, Universidade Federal de Minas Gerais
 Sumário

Capítulo 1. Introdução 5
Capítulo 2. Introdução aos Grafos: tipos de grafos e propriedades 7
1. Tipos de grafos 8
2. Análise descritiva dos grafos 16
3. Usando R 23

Capítulo 3. Grafos Aleatórios: Modelos e Inferência 27

1. Modelos de Grafos Aleatórios 27
2. Estimadores de Máxima Verosimilhança (ERGM) 36
3. Estimadores Bayesianos (ERGM) 39
4. Usando R 42

Capítulo 4. Modelos de Regressão para Grafos e Modelos Dinâmicos 45

1. Modelos de Regressão ou Previsão para Grafos 45
2. Outros Tópicos 54
3. Usando R 56
Referências 61

3
 CAPíTULO 1

Introdução

Estas notas tem por objetivo discutir os conceitos básicos de rede ou grafo aleatório
e das aplicações ao estudo das redes sociais. Modelos probabilísticos de redes aleató-
rias. Estatísticas descritivas das redes aleatórias: grau, subgrafos, conetividade, etc .
Estimadores de máxima verosimilhança e Bayesianos assim como a implementação no
software R.
Usaremos alguns algoritmos dos programas igraph, network, ergm e bergm. para
serem utilizados em alguns exemplos.
Os objetivos deste trabalho são:
• Introdução aos Grafos, tipos de grafos, propriedades.
• Análises Descritivas dos grafos.
• Modelos de grafos Aleatórios, Grafos Aleatórios Exponenciais.
• Estimadores de Máxima Verosimilhança para grafos aleatórios.
• Estimadores Bayesianos para grafos aleatórios.
• Tópicos de regressão linear e previsão. Outros Tópicos.
Estas notas estão baseadas no livro Statistical Analysis of Network Data: Methods
and Models, do Kolaczyk, Eric D., ver [7].

5
 CAPíTULO 2

Introdução aos Grafos: tipos de grafos e

propriedades

Um famoso problema histórico de matemática, que se tornou uma legenda popular

conhecida como as sete pontes de Königsberg.
Cortado pelo rio Prególia, a cidade de Königsberg (território de Prússia) possuía
duas grandes ilhas que, juntas, formaram um complexo que naquela época continha
sete pontes. Foi discutida a possibilidade de atravessar todas as pontes sem repetir
nenhuma ponte. Em 1736, o matemático Leonard Euler provou, de forma simples, que
não havia como satisfazer essas condições. Criando o primeiro gráfico da história, Euler
transformou as estradas em linhas retas e suas interseções em pontos. primeiro grafo
de la historia, Euler transformou os caminhos em retas e suas interseções em pontos.

Figura 1. As sete pontes de Königsberg

Algumas aplicações de grafos:

• Redes sociais: Facebook, Twitter, etc.
• Epidemiologia, saúde pública.
• Comunicações: celulares, redes de computadores
• Biologia molecular, genética.
• Inteligência Artificial e novos paradigmas da teoria da informação.
• Etc.
Para resolver problemas práticos, como os mencionados acima, um grafo é a melhor
solução para representar a relação entre os objetos de um determinado conjunto, onde
alguns pares de objetos, vértices, também chamados de nós ou pontos, são conectados
por bordas, também conhecidas como linhas ou arcos.
Tal como as equipes de futebol, ou estilos de música e gêneros cinematográficos,
os grafos também possuem características que os separam em diferentes categorias,
como serão apresentadas a seguir.
7
8 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

0.1. Definição de grafos. Considere o conjunto de vértices V , que é finito com

n elementos, V = {1, . . . , n}. Para este conjunto vamos definir o conjunto de arestas
E ⊆V ×V.
Definição 0.1. Um grafo G é um par ordenado de conjuntos: G = (V, E). O
conjunto de todos os possíveis grafos com n vértices será denotado por Gn , isto é:
Gn = {G : G = (V, E), |V | = n, E ⊆V ×V}
Exemplo 0.2. Considere o conjunto de vértices V = {1, . . . , 7} e as arestas
E = {(1, 4), (1, 5), (2, 3), (2, 7), (3, 6)(4, 1), (4, 7), (5, 6), (5, 3, )(7, 7)}, então, podemos
representar o grafo G = (V, E) por médio de pontos e setas:

Figura 2. Grafo G com n = 7 vértices

O tamanho de um grafo é o número de arestas do grafo e a ordem é o numero

de vértices, assim, no exemplo anterior o tamanho do grafo é |E| = 10 e a ordem
é |V | = 7. Usualmente os grafos classificam-se pelo tamanho em dispersos (sparce
graph) se |E| ≈ |V | ou densos se |E| ≈ |V |2 .

1. Tipos de grafos
Um grafo não dirigido é aquele em que os vértices não são ordenados, o seja, as
arestas não possuem orientação. Assim, por exemplo, se suprimimos no grafo anterior
as direciones, então obtemos:
 1. TIPOS DE GRAFOS 9

Figura 3. Grafo G não dirigido com n = 7 vértices

Este grafo é um exemplo de um multigrafo com loops, isto é, vértices com múltiplas
arestas e arestas que conectam o mesmo vértice. A seguir vamos considerar apenas
grafos que não tem estas características. Para isso vamos definir um grafo simples
como aquele grafo que não tem loops ou arestas múltiplas, em particular vamos também
considerar nestas notas somente grafos não dirigidos.

Figura 4. Grafo simples G (não dirigido) com n = 7 vértices

10 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

No caso dos grafos simples com n vértices, vale que o número destes grafos é:
n
|Gn | = 2( 2 ) , pois o número de possíveis pares de vértices é n2 e o número de sub-


n
conjuntos deste conjunto de pares é 2( 2 ) ; lembre que um grafo é um subconjunto de
pares. Vamos a definir o grau d(i) de um vértice i ∈ V como o número de vértices
conectados com este vértice por médio de uma aresta que pertence ao grafo.

X
d(i) = 1.
j:(i,j)∈E

Assim no grafo da Figura 4 temos que, d(1) = 2, d(2) = 2, etc. Podemos representar
esta função por médio do seguinte gráfico:

Figura 5. Grau do grafo G

A distribuição do grau de um grafo é a frequência relativa dos graus do grafo. No

exemplo anterior somente ha vértices com graus 2 (com frequência relativa de 0, 7) e 3
(com frequência relativa de 0, 3); através de um gráfico a representamos na Figura 6.
 1. TIPOS DE GRAFOS 11

Figura 6. A distribuição do grau do grafo G

Definimos a media do grau como:

1 X
d(G) = d(v).
|V |
v∈V

Então, vale que:

1
|E| = d(G)|V | (1.1)
2
Dizemos que dois vértices são adjacentes se uma aresta do grafo conecta a ambos.
Outra maneira de representar um grafo é usar esta noção de adjacência para definir
a matriz de adjacência A = (ai,j ), na cual as linhas e as colunas representam os
vértices do grafo, e os coeficientes ai,j = 1 ou 0 de acordo com as correspondentes
linha i e coluna j sejam adjacentes, isto é, formem uma aresta do grafo: (i, j) ∈ E.
No grafo do exemplo anterior a correspondente matriz de adjacência A seria:

> A
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0 1 1 0 0 0 0
[2,] 1 0 0 0 0 1 0
[3,] 1 0 0 0 1 0 1
[4,] 0 0 0 0 1 1 0
[5,] 0 0 1 1 0 0 1
[6,] 0 1 0 1 0 0 0
[7,] 0 0 1 0 1 0 0

Figura 7. Matriz de adjacência A

Observe que esta matriz A é simétrica pois o grafo G é simples, logo não é dirigido.
12 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

Outra matriz importante é a matriz Laplaciana, do grafo: L = D − A, com D a

matriz diagonal com o grau de cada vértice na diagonal: D =diag(d(i)).

> L
1 4 5 2 3 7 6
1 2 -1 -1 0 0 0 0
4 -1 2 0 0 0 -1 0
5 -1 0 3 0 -1 0 -1
2 0 0 0 2 -1 -1 0
3 0 0 -1 -1 3 0 -1
7 0 -1 0 -1 0 2 0
6 0 0 -1 0 -1 0 2

Figura 8. Laplaciano L

Dizemos que L está normalizado se

( a
− √ i,j si i 6= j
Li,j = d(i)d(j)
1 si i = j
p
neste caso, L = D−1/2 (I − A)D−1/2 com D−1/2 =diag(1/ d(i)) (onde colocamos
zero se d(i) = 0), do exemplo anterior temos:

> L
1 4 5 2 3 7 6
1 1.0000000 -0.5 -0.4082483 0.0000000 0.0000000 0.0 0.0000000
4 -0.5000000 1.0 0.0000000 0.0000000 0.0000000 -0.5 0.0000000
5 -0.4082483 0.0 1.0000000 0.0000000 -0.3333333 0.0 -0.4082483
2 0.0000000 0.0 0.0000000 1.0000000 -0.4082483 -0.5 0.0000000
3 0.0000000 0.0 -0.3333333 -0.4082483 1.0000000 0.0 -0.4082483
7 0.0000000 -0.5 0.0000000 -0.5000000 0.0000000 1.0 0.0000000
6 0.0000000 0.0 -0.4082483 0.0000000 -0.4082483 0.0 1.0000000

Figura 9. Laplaciano normalizado L

O espectro do Laplaciano é o conjunto de auto-valores do Laplaciano normalizado,

neste exemplo:

> Espectro
[1] 0.00 0.34 0.53 1.20 1.39 1.62 1.91
Figura 10. Espectro de L

É possível provar que os autovalores estão limitados entre 0 e 2, e que o primeiro

autovalor é sempre 0. O número de subgrafos conexos do grafo corresponde ao número
de autovalores zero. Os maiores autovalores estão associados a formação de ‘anéis‘ no
grafo.
 1. TIPOS DE GRAFOS 13

1.1. Classes importantes de grafos: Um grafo regular é um grafo no qual todos

os vértices possuem o mesmo número de ligações com os outros vértices, o seja, todos os
vértices tem o mesmo grau. Um grafo regular com vértices de grau k é chamado grafo
k-regular ou gráfico regular de grau k. Por exemplo, G com 5 vértices é 2-regular:

Figura 11. Grafo 2-regular G

Este grafo também é chamado de 2-ciclo, C2 , e é um caminho fechado que não

tem interseções. Agora considere G com 5 vértices é 4-regular:

Figura 12. Grafo 4-regular G

Este grafo é um grafo completo, K5 , é dizer que cada vértice tem uma aresta
ligando com todos os outros vértices. Os grafos completos tem a característica de que
14 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

cada par de vértices tem uma aresta que os une. Finalmente, observe nestes exemplos
que um grafo com n vértices e k-regular não podemos ter n e k impares, um deles tem
que ser par (use a relação (1.1)).
Um caminho é um conjunto de vértices v1 , . . . , vk ∈ V tal que o conjunto de
arestas (v1 , v2 ), (v2 , v3 ), . . . , (v`−1 , v` ), (v` , v`+1 ), . . . , (vk−1 , vk ) está em E; uma ge-
odésica entre dos vértices é o caminho de menor tamanho entre eles. Um grafo é
chamado de grafo conexo o conectado se todo par de vértices do grafo pode unir-se
por um caminho. O grafo da Figura (4) é conexo:
> is.connected(G)
[1] TRUE
Mais, se removemos as arestas (1, 5) e (2, 3), obtemos um grafo desconexo:

Figura 13. Grafo desconexo G

verificamos que é desconexo:

> is.connected(G)
[1] FALSE
Também podemos verificar isto usando o espectro deste grafo, é possível mostrar
que neste caso terá dois autovalores zero:

> Espectro
[1] 0.0 0.0 0.5 1.5 1.5 1.5 2.0
Figura 14. Espectro de L

Um grafo bipartido é um grafo no cual os vértices estão separados em dois con-

juntos, V1 e V2 e as arestas somente conectam vértices de conjuntos separados. Na
Figura (15) podemos observar um grafo bipartido:
 1. TIPOS DE GRAFOS 15

Figura 15. Grafo bipartido g

Uma árvore é um grafo que não contem ciclos será k-árvore se for k-regular. Por
exemplo:

Figura 16. 3-árvore com 13 vértices T3

Os vértices de grau 1 são chamados folhas.

Árvores desconexos são conhecidos como florestas
Uma n-estrela, Sn , é um grafo que é uma árvore com n vértices e n − 1 folhas.
Por exemplo:
16 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

Figura 17. 6-estrela S6

Finalmente, um triângulo é um grafo completo com 3 vértices, K3 :

Figura 18. Triângulo

2. Análise descritiva dos grafos

As propriedades estruturais dos grafos são características de interesse em dois sen-
tidos. Em primeiro lugar, desde o ponto de vista das aplicações as redes sociais se
interessam com a conectividade do grafo e o papel de certos vértices centrais. Desde
o ponto de vista matemático, é possível mostrar que a quantidade de certos subgrafos
caracteriza a um grafo. Vamos a descrever três destas características: centralidade,
transitividade e conectividade.
 2. ANÁLISE DESCRITIVA DOS GRAFOS 17

2.1. Centralidade. A centralidade de um vértice é caracterizada pela importância

desse vértice num grafo. Podemos considerar, em principio, que basta observar seu grau,
pois esse é um indicativo de quantos vértices estão conectados com ele. Observe que
0 ≤ d(v) ≤ |V | − 1. Assim, teríamos que os vértices com maior grau são mais centrais.
Vejamos um exemplo, considere o grafo G da Figura (19 ), com |V | = 20 vértices:

Figura 19. O grafo G

O grau de cada vértice é:

> degree(G)
6 7 1 8 9 4 10 5 12 11 13 14 15 17 18 16 2 19 3 20
3 2 2 3 4 2 1 3 2 1 1 3 2 4 2 1 1 1 0 0

e a distribuição do grau é:
18 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

Figura 20. Grau do grafo G

como podemos observar existem 2 vértices com grau 4, que é o máximo. Estos
dois vértices são os mais conectados.
Outra medida central para vértices é a centralidade por proximidade:

1
cC (v) = P , v∈V (2.2)
u∈V d(v, u)

Donde d(v, u) é a distancia (número de arestas) da geodésica entre u e v, aqui

iremos supor que o grafo é conexo.
Os valores da proximidade para o maior subgrafo conexo do exemplo anterior (G1 ⊂
G, V (G1) = V (G) − {2, 3, 19, 20}) são:

> closeness(G)
6 7 1 8 9 4 10
0.01923077 0.01612903 0.01851852 0.02325581 0.02380952 0.02000000 0.01515152
5 12 11 13 14 15 17
0.01886792 0.01538462 0.01265823 0.01492537 0.02380952 0.02173913 0.02222222
18 16
0.01351351 0.01136364

ou graficamente também pode ser observado:

 2. ANÁLISE DESCRITIVA DOS GRAFOS 19

Figura 21. Proximidade dos vértices do grafo G1 ⊂ G

A Figura (21) mostra que os vértices 9 e 14 são os que estão mais próximos aos
demais vértices.
Agora, podemos considerar que um vértice é central usando a media de todos os
caminhos que usam esse vértice, a centralidade por mediação está definida como:

X σ(s, t|v)
cM (v) = , v ∈ V. (2.3)
σ(s, t)
s6=u6=v∈V

Onde σ(s, t|v) é o número de geodésicas entre sPe t que passam por v e σ(s, t) é
uma constante normalizadora definida por σ(s, t) = v σ(s, t|v).
Os valores da mediação para G1 são:

> betweenness(G)
6 7 1 8 9 4 10
47.0000000 26.0000000 0.3333333 23.5000000 31.5000000 0.0000000 0.0000000
5 12 11 13 14 15 17
38.0000000 14.0000000 0.0000000 0.0000000 54.3333333 50.0000000 45.3333333
18 16
14.0000000 0.0000000

graficamente:
20 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

Figura 22. Mediação dos vértices do grafo G1 ⊂ G

Aqui, observamos que os vértices 9 e 14 são os que tem um papel central maior
que os demais vértices.
Estas duas definições de centralidade podem ser estendidas para arestas e ainda
existe uma noção de centralidade usando o espectro.

2.2. Transitividade. A transitividade está relacionada com uma característica das

redes sociais que poderíamos resumir como: “o amigo de meu amigo é meu amigo”;
esta propriedade costuma aparecer em redes coesivas, é dizer em redes com grau alto
de agregação.
Num grafo simples a transitividade é também chamada o coeficiente de agre-
gação e está definida como:

3τK3 (G)
cl(G) = . (2.4)
τS2 (G)

Onde τK3 (G) é o número de triângulos, K3 , contidos no grafo e τS2 (G) é o

número de 2-estrelas, S2 , contidas no grafo. Observe que cada 2-estrelas é um possível
triângulo, basta que esteja a aresta que une as dois pontas da estrela para formar
o triângulo. logo o coeficiente de agregação, cl(G), mede a proporção de possíveis
triângulos que são triângulos. Como cada 2-estrela é contabilizada 3 vezes, em cada
possível triângulo, teríamos que dividir para 3 o número delas, τS2 (G)/3.
O coeficiente de agregação é um dos mais conhecidos índices para grafos. Vejamos
qual é o coeficiente de agregação nos grafos, G e G1 dos exemplos:
 2. ANÁLISE DESCRITIVA DOS GRAFOS 21

Figura 23. Transitividade dos grafos G e G1, cl(G) = 0.27,

cl(G1) = 0.1.

Claramente o grafo G1 é menos agregado, pois tem más vértices e porque entre
estes vértices extras tem más 2-estrelas que não são triângulos.

2.3. Conectividade. É de interesse nas aplicações em grafos a localização das

comunidades dentro de um grafo, isto é, subgrafos que tenham seus vértices conexos,
em algum sentido. Vamos a ver como se faz isto usando uma técnica de hierarquização.
Suponha que para um grafo G = (V, E) existe uma partição de V em subconjuntos
disjuntos V = {V1 , . . . , VK }. A partir de V definimos fi,j = fi,j (V) como a fração de
arestas de E que tem um vértice em Vi e outro vértice no Vj . Usando fi,j definimos a
modularidade de V como:

K
X  ∗
2
mod (V) = fk,k − fk,k (2.5)
k=1

∗
PK PK
Onde fi,j = fk,+ f+,k e fk,+ = j=1 fk,j , f+,k = i=1 fi,k . Observe que num
grafo simples, isto é, não dirigido, vale que fk,+ = f+,k . Esta seleção de f ∗ corresponde
a um grafo com o mesmo grau que o grafo original, mas com as arestas colocadas
aleatoriamente, sem levar em consideração a partição V. O problema consiste agora
em encontrar entre todas as partições V, com um número fixo K de elementos, aquela
que minimiza mod (V).
Segue um exemplo, considere o grafo G1 do exemplo da seguinte Figura (24):
22 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

Figura 24. Grafo G1.

Usando o algoritmo ‘fastgreedy’ do igraph, obtemos:

Figura 25. Partição do grafo G1 em quatro ‘comunidades’.

 3. USANDO R 23

3. Usando R
Nesta seção descrevemos os comandos do R que foram utilizados para representar
as figuras desta seção Observe que é necessário carregar a libraria do pacote igraph do
R para poder usar os comandos.
(1) Figura 2:
G <- graph.formula(1-+4, 1-+5, 2-+3, 2-+7, 3-+6, 4-+1,
4-+7, 5-+6, 5-+3, 7-+7,simplify = FALSE)
plot(G,edge.color="black",layout=layout.fruchterman.reingold)
(2) Figura 3:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3, 7-7,simplify = FALSE)
plot(G,edge.color="black",layout=layout.fruchterman.reingold)
(3) Figura 4:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
plot(G,edge.color="black",layout=layout.fruchterman.reingold)
(4) Figura 5:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
degree(G)
plot(V(G),degree(G))
(5) Figura 6:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
degree_distribution(G)
max(degree(G))
degree_distribution(G)
plot(0:max(degree(G)),degree_distribution(G),
xlab = "grau(G)",ylab = "frecuencia")
(6) Figura 7:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
A<-as_adj(G, type = "both", attr = NULL,edges = FALSE,
names = FALSE, sparse =FALSE)
A
(7) Figura 8:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
L<-laplacian_matrix(G, normalized = FALSE, weights = NULL,
sparse = FALSE)
L
(8) Figura 9:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1,
4-7, 5-6, 5-3)
L<-laplacian_matrix(G, normalized = TRUE, weights = NULL,
sparse = FALSE)
L
(9) Figura 10:
24 2. INTRODUÇÃO AOS GRAFOS: TIPOS DE GRAFOS E PROPRIEDADES

G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1, 4-7, 5-6, 5-3)
L<-laplacian_matrix(G, normalized = TRUE, weights = NULL,
sparse = FALSE)
Espectro<-sort(round(eigen(L)$values,2))
Espectro
(10) Figura 11:
G<-sample_k_regular(5, 2, directed = FALSE, multiple = FALSE)
plot(G)
(11) Figura 12:
G<-sample_k_regular(5, 4, directed = FALSE, multiple = FALSE)
plot(G)
(12) Figura 13:
G <- graph.formula(1-4, 2-7, 3-6, 4-1, 4-7, 5-6, 5-3)
plot(G,edge.color="black",layout=layout.fruchterman.reingold,
vertex.size=20,edge.arrow.size=20)
(13) Figura 15:
g <- sample_bipartite(3,4, type = "gnp", p=0.5, directed = FALSE)

plot(g, layout = layout_as_bipartite,

vertex.color=c("yellow","cyan")[V(g)$type+1],
vertex.size=20)
(14) Figura 14:
G <- graph.formula(1-4, 2-7, 3-6, 4-1, 4-7, 5-6, 5-3)
L<-laplacian_matrix(G, normalized = TRUE, weights = NULL,
sparse = FALSE)
Espectro<-sort(round(eigen(L)$values,2))
Espectro
(15) Figura 16:
T<-make_tree(13, 3, mode = "undirected")
plot(T,vertex.size=20,layout=layout_as_tree(T, root=c(1)))
(16) Figura 17:
S<-make_star(6, mode = "undirected")
plot(S,layout=layout.fruchterman.reingold,vertex.size=20 )
(17) Figura 18:
G<-graph.formula(1-- 2,2-- 3,3-- 1)
plot(G,edge.color="black",layout=layout.fruchterman.reingold,
vertex.size=20,edge.arrow.size=20)
(18) Figura 19:
G<-graph.formula(6-- 7,1-- 8,1-- 9,4-- 9,6--10,5--12,11--12,
5--13,8--14,9--14,6--15,14--15,4--17,5--17,8--17,9--17,7--18,
16--18, 2--19)
G<-add_vertices(G,2)
V(G)$name[19:20]<-c(3,20)
plot(G,layout=layout.fruchterman.reingold,vertex.size=20 )
(19) Figura 20:
grau<-degree(G)
hist(grau,col="cyan",
xlab="grau", ylab="Frecuencia",main = "",
breaks = c(0:(max(grau)+1)-0.5), xlim = c(-1,(max(grau)+1)))
 3. USANDO R 25

(20) Figura 21:

G1<-graph.formula(6-- 7,1-- 8,1-- 9,4-- 9,6--10,5--12,11--12,5--13,
8--14,9--14, 6--15,14--15,4--17,5--17,8--17,9--17,7--18,16--18)
closeness(G1)
plot(V(G1)$name,closeness(G1),col="blue",
ylab="cercanía(G1)", xlab="vértices",main = "")
(21) Figura 22:
G1<-graph.formula(6-- 7,1-- 8,1-- 9,4-- 9,6--10,5--12,11--12,5--13,
8--14,9--14,6--15,14--15,4--17,5--17,8--17,9--17,7--18,16--18)
betweenness(G1)
plot(V(G1)$name,closeness(G),col="blue",
ylab="cercanía(G1)", xlab="vértices",main = "")
(22) Figura 23:
G <- graph.formula(1-4, 1-5, 2-3, 2-7, 3-6, 4-1, 4-7, 5-6, 5-3)
transitivity(G)
G1<-graph.formula(6-- 7,1-- 8,1-- 9,4-- 9,6--10,5--12,11--12,5--13,
8--14,9--14, 6--15,14--15,4--17,5--17,8--17,9--17,7--18,16--18)
transitivity(G1)
(23) Figura 25:
G1<-graph.formula(6-- 7,1-- 8,1-- 9,4-- 9,6--10,5--12,11--12,5--13,
8--14,9--14,6--15,14--15,4--17,5--17,8--17,9--17,7--18,16--18)
plot(G1)
kc<-fastgreedy.community(G1)
length(kc)
sizes(kc)
membership(kc)
plot(kc, G1)
 CAPíTULO 3

Grafos Aleatórios: Modelos e Inferência

Neste capítulo vamos a discutir alguns modelos clássicos de grafos aleatórios. Os

primeiros modelos foram introduzidos por Erdös e Rényi em 1959, eram grafos simples
em que cada aresta era incluída no grafo com igual probabilidade.
Existem basicamente dois tipos de modelos; um é mais teórico-matemático e é
usado para representar, qualitativamente, fenômenos de grandes redes como por exem-
plo, os grafos de ‘escala livre’ ou ‘o mundo pequeno’, como no modelo de ‘conexão
preferencial’. Outro tipo de modelos mais usados nas aplicações são modelos paramé-
tricos, como o modelo exponencial, que é o modelo que vamos discutir aqui.

1. Modelos de Grafos Aleatórios

Vamos a considerar os grafos simples, é possível estender estas definições para
outros grafos dirigidos, multígrafos, bipartidos.
Definição 1.1. Um modelo aleatório para grafos esta definido pela coleção:
{Pθ (G), G ∈ Gn , θ ∈ Θ}.
Onde Gn = {G = (E, V ), |V | = n, G simples} é uma coleção de grafos simples de
ordem n que é chamada espaço amostral, e para cada θ ∈ Θ:
Pθ : Gn → [0, 1],
é uma distribuição de probabilidades em Gn . Usualmente o espaço paramétrico Θ é
um conjunto ‘regular’ contido em algum Rk
A principal dificuldade de trabalhar com modelos aleatórios para grafos consiste
em que devemos definir o valor da probabilidade P para cada grafo G = (E, V ), (com
|V | = n), do espaço amostral, Gn , e o número de grafos deste espaço é finito mas
n
muito grande, |Gn | = 2( 2 ) ≈ 2n e o número de vértices n é grande.
2

Outra característica determinante nestos modelos é o tamanho dos grafos que

estamos considerando, grafos dispersos (sparce graph) se |E| ≈ |V | ou densos se
|E| ≈ |V |2 .
Esta dificuldade e outras considerações nos levam a restringir a definição de P à
um conjunto menor, para isso consideramos certas funções chamadas estatísticas (ou
observações) η do grafo:
η : Gn → R.
Por exemplo, considerando η(G) = |E| = ‘o número de arestas do grafo’, podemos
restringir Gn aos grafos G que tenham um tamanho fixo m, isto é, o novo espaço
amostral G ∗ sera:
G ∗ = {G = (E, V ), |V | = n, η(G) = m, G simples}.
27
28 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

Exemplo 1.2. Considere n = 3 vértices, neste caso G3 = {G1, . . . , G8} está

3
formado por 2(2) = 8 grafos:

Figura 1. Os 8 grafos contidos em G3

Podemos definir uma distribuição de probabilidades em G3 fazendo:

θi θi (1 − θ)
Pθ (Gi) = P8 = , i = i = 1, . . . , 8, θ ∈ R+ .
k=1 θ
k θ(1 − θ7 )

5
θ (1−θ)
Assim P(G5) = θ(1−θ 7 ) . Observe que se θ < 1 então os grafos mais densos (G8

por exemplo) tem menos probabilidade. No entanto se θ > 1 ocorre o oposto.

Agora se fazemos η = |E| e fixamos o conjunto de grafos que tem η = 2 vértices,
isto é, estamos interessados em uma característica ‘estrutural’: ter dois vértices. Assim
obtemos G ∗ = {G5, G6, G7}
 1. MODELOS DE GRAFOS ALEATÓRIOS 29

Figura 2. G5, G6, G7

Neste caso poderíamos calcular as probabilidades destes grafos, usando a probabi-

lidade condicional, por exemplo:
θ 7 (1−θ)
∗ Pθ (G7 ∩ η = 2) θ(1−θ 7 ) θ7
P (G7) = Pθ (G7|η = 2) = = P7 θk (1−θ) = 5 .
Pθ (η = 2) 7
θ + θ6 + θ7
k=5 θ(1−θ )

Observe que o evento {η = 2} está formado pelos grafos com dois arestas e que nosso
espaço amostral tem agora somente 3 grafos.
Poderíamos também escolher outras características η como o número de triângulos,
o grau, etc. e fixar este número.
1.1. Modelos Clássicos. Nesta seção descrevemos quatro modelos: Uniforme,
Erdös-Rényi, Conexão Preferencial e Exponencial.
1.1.1. Modelo Uniforme. Neste modelo todos los grafos tem igual probabilidade:
1
PU (Gi) = .
|Gn |
Para este modelo η é uma Variável Aleatória com função de distribuição de
probabilidade:
#{G ∈ Gn , η(G) ≤ t}
Fη,Gn (t) = ,
|Gn |
onde #{·} é a cardinalidade de um conjunto e 0 ≤ t ∈ R. Esta Variável Aleatória
consiste em que vamos à observar esta característica η para cada grafo de Gn e medimos
que tão provável é ela.
Exemplo 1.3. Considere o espaço amostral G3 dos grafos com três vértices. Neste
caso é fácil ver que:
1
PU (Gi) = , i = 1, . . . , 8
8
Por outro lado, observe que η(G) = |E(G)| = 0, 1, 2, 3 é uma variável aleatória
discreta com distribuição:
1 3
P(η = 0) = , P(η = 1) = ,
8 8
30 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

3 1
P(η = 2) = P(η = 3) = .
8 8
Em geral, suponha que o número de vértices é fixo, |V | = n e que também
fixamos o número de arestas, |E| = m. Definimos o modelo Uniforme com m arestas,
G(n, m), como a distribuição de probabilidades que atribui igual probabilidade aos
grafos pertencentes à Gn,m = {G = (E, V ) : |V | = n, |E| = m}, isto é:
1
P(G) = N

, G ∈ Gn,m ,
m

com N = n2 o número de arestas possíveis num grafo com n vértices.

Para simular grafos do modelo G(n, m) podemos usar o programa R. Um exemplo,

simulado, de um grafo aleatório Uniforme com |V | = 30 e |E| = 100 fixo:

Figura 3. Grafo G Uniforme com |V | = 30 e |E| = 100.

1.1.2. Modelo Erdös-Rényi. O modelo de Erdös-Rényi é uma generalização do mo-

delo uniforme e é o modelo mais simples de grafo aleatório. Foi introduzido em 1959
por Paul Erdös em um dos trabalhos mais citados sobre combinatória. Erdös formulou
neste trabalho a ideia de usar probabilidade positiva de um evento para provar um
resultado de existência na matemática combinatória, para isso introduz pela primeira
vez na matemática a ideia de um grafo aleatório. Este método de prova será conhecido
depois como o ‘Método Probabilístico’.
A construção simplesmente consiste em dar uma probabilidade p ao evento de que
cada aresta pertença ao grafo, e isto independentemente das outras arestas: P(e ∈
G) = p.
 1. MODELOS DE GRAFOS ALEATÓRIOS 31

A probabilidade para G = (E, V ) ∈ Gn no modelo de Erdös, G(n, p), é então:

n
Pp (G) = p|E| (1 − p)( 2 )−|E| .
Observe que a probabilidade de cada grafo G com m arestas seria P(G) = pm (1 −
n
p)( 2 )−m .
Podemos demostrar uma relação entre os modelos Uniforme G(n, m) e Erdös-Rényi
G(n, p). Se a ‘densidade de vértices’ se mantêm igual a p = m então ambos modelos
(n2 )
coincidem quando n é suficientemente grande, isto ocorre pela lei dos grandes números.
Para simular grafos do modelo G(n, p) podemos usar o programa R. Um exemplo
de um grafo simulado para o modelo de Erdös-Rényi, com |V | = 10:

Figura 4. Grafo G Erdös-Rényi com |V | = 10, p = 0.3, p = 0.5 e

p = 0.7

n

Observe que o número de arestas neste modelo é aleatório, em media é E|E| =
2 p. Assim na figura (4) o número médio de arestas é:
Ep=0.3 |E| = 13.5, Ep=0.5 |E| = 22.5 e Ep=0.7 |E| = 31.5
1.1.3. Modelo Conexão Preferencial. Esta família de modelos é uma das mais im-
portantes na literatura de modelos de redes, pois representa de maneira aproximada o
comportamento de redes como a internet. Este é um modelo dinâmico, no sentido de
que o número de vértices não é fixo, neste modelo vamos dinamicamente adicionando
arestas e vértices com um mecanismo que descreveremos a continuação.
Uma das propriedades importantes deste modelo é que a distribuição do grau ‘decai’
como uma lei de potencias, (lembre que nossos grafos são aleatórios logo o grau dos
vértices de cada grafo não é fixo), isto é, se d(vk ) é o grau do k-essimo vértice, quando
os graus estão ordenados em forma crescente, então:
1
d(vk ) ≈ α .
k
Este fenômeno é chamado de escala libre (scale-free), e significa que com probabi-
lidade pequena, não nula, existem vértices com um grau muito grande, isto é, vértices
32 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

que estão conectados com uma grande quantidade de outros vértices. Este fenômeno
é conhecido como ‘mundo pequeno’ (small-world), pois significa que usando caminhos
pequenos é possível conectar a todos os vértices. Em estúdios empíricos, tem-os en-
contrado que o tamanho desses caminhos é 6, como consequência nas redes sociais
necessitamos apenas 6 pessoas para conectarnos com qualquer outra pessoa.
A seguir veremos um modelo de grafos que possuiem estas características, é o
modelo de ‘Barabási–Albert’, ou modelo BA. A ideia é criar vértices que estejam co-
nectados aos vértices com grau maior:

• A rede começa com n0 vértices (|V (Gn )| = n0 .

• Um novo vértice in é criado em cada passo n. Com probabilidade proporcional
ao grau de cada uno dos vértices já existentes, assim, criamos uma aresta que
une este vértice com os anteriores:

d(j)
P((in , j) ∈ V (Gn+1 )) = P , j ∈ V (Gn ).
k∈V (Gn ) d(k)

Vejamos um exemplo simulado de um grafo BA com n = 500 vértices e

com poderes α = 0.1 e α = 0.001, α é a potencia del decaimento.:

Figura 5. Grafo G do modelo BA com |V | = 500, α = 0.1 e α = 0.001

Histogramas dos graus para os grafos anteriores, são mostrados a seguir,

 1. MODELOS DE GRAFOS ALEATÓRIOS 33

Figura 6. Histogramas dos Grafos G do modelo BA com |V | =

500, α = 0.1 e α = 0.001

1.1.4. Modelo Exponencial Aleatório ERGM. O Modelo Exponencial Aleatório

para Grafos (ERGM) é o modelo paramétrico mais usado nos estudos de redes soci-
ais, sua importância radica em que permite introduzir observações das características
estruturais dos grafos na distribuição de probabilidade, isto é, a probabilidade de um
grafo depende de características tais como o número de arestas, triângulos, estrelas,
etc. que estão presentes no grafo.
Este modelo foi introduzido por Holland e Leinhardt em 1981, ver [5], e em sua
forma general por Frank e Strauss em 1986, ver [4]. O estudo matematico e rigoroso
deste modelo foi desenvolvido por Diaconis et. al. em 2013, ver [2], nesse trabalho
mostra-se a lei dos grades números para distribuições de grafos , assim como também
a transição de fase para uma família importante de ERGMs.
Para descrever este modelo, considere um conjunto de características: η1 , . . . , ηk de
um grafo, estas podem ser o número de arestas, triângulos, estrelas, etc. e suponha que
associamos parâmetros: θ1 , . . . , θk a cada uma delas. Então definimos a distribuição
de probabilidades deste modelo como:

Pk
θi ·ηi (G)
e i=1
PΘ (G) = , G ∈ Gn ,
Z(Θ)
onde Θ = (θ1 , . . . , θk ) são os parâmetros e Z(Θ) é uma constante normalizadora
chamada função de partição:

X Pk
θi ·ηi (G)
Z(Θ) = e i=1 .
G∈Gn

O grafo de Erdös-Rényi é um caso particular de este modelo, ele corresponde a

usar somente o número de vértices do grafo no modelo:

eθ·|E(G)|
PΘ (G) = , G ∈ Gn .
Z(θ)
34 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

Observe que η(g) = |E(G)|. Neste caso, é fácil achar o valor da constante norma-
lizadora: Z(θ) = (1 + eθ )n . Assim, temos que a probabilidade de ter uma aresta é
eθ
p = 1+e θ.

Em geral é um problema difícil encontrar o valor de Z(Θ), no lugar disso, são

utilizados algoritmos aproximados, ao menos quando n é pequeno.

Exemplo 1.4. Considere os grafos com n = 10 vértices, usaremos duas caracte-

rísticas: o número de arestas (η1 (G)) e o número de triângulos (η2 (G)) presentes no
grafo.
Vejamos que ocorre se fixamos o parâmetro que controla o número de triângulos
θ2 = 0.3 e variamos θ1 (θ1 = −1, θ1 = 1 ) que controla o número de arestas. Lembre
eθ1
que a probabilidade de ter uma aresta é p = 1+e θ1 (obtemos p1 = 0.27 e p2 = 0.73)

Usando o programa ‘ergm’ de R obtemos os grafos G1 (θ1 = −1, θ2 = 0.3) e G2

(θ1 = −1, θ2 = 0.3):

Figura 7. Grafos G1 e G2 do modelo ERGM com θ1 = −1, 1, θ2 = 0.3

Observe o número de arestas e triângulos que se formaram:

> summary(G1 ~ edges + triangle)

edges triangle
13 2

> summary(G2 ~ edges + triangle)

edges triangle
41 90

Isto significa que ao aumentar as arestas (edges) aumentamos o número de triân-

gulos, isto é ambas características estão correlacionadas no grafo
Os histogramas dos graus para os grafos anteriores:
 1. MODELOS DE GRAFOS ALEATÓRIOS 35

Figura 8. Distribuição do grau de G1 e G2

Como vemos o grafo G2 é mas conectado, seu índice de agregação (transitividade

o densidade de triângulos) é quatro vezes maior que o de G1:

> transitivity(G1)
[1] 0.2
> transitivity(G2)
[1] 0.9060403

1.1.5. Simulação de Grafos Aleatórios. Os algoritmos para simular modelos de gra-

fos aleatórios tem o problema relacionado com o número extremadamente grande de
grafos que pertencem ao espaço amostral. Por exemplo, no caso do modelo ERGM,
pode-se percever na função de partição Z(Θ):

X Pk
θi ·ηi (G)
Z(Θ) = e i=1 .
G∈Gn

2
A soma tem 2n termos. Por isto existem algoritmos aproximados para encontrar
Z(Θ), vamos a explicar brevemente em que consiste um deles.
Neste algoritmo usamos a Lei de los grandes números para encontrar Z(Θ), usando
o fato de que podemos simular com facilidade grafos de Erdös-Rényi, com p = 0, 5,
estes grafos correspondem à Θ0 = (0, . . . , 0).
36 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

Agora observamos a seguinte relação:

X e ki=1 θi ·ηi (G)
P
Z(Θ)
= ,
Z(Θ0 ) Z(Θ0 )
G∈Gn
Pk Pk
θi ·ηi (G)
X e i=1 e i=1 0·ηi (G)
= ,
Z(Θ0 )
G∈Gn
X Pk
θi ·ηi (G)
= PΘ0 (G)e i=1 ,
G∈Gn
X t
·η(G)
= PΘ0 (G)eΘ ,
G∈Gn
t
·η(G)
= EΘ0 eΘ .
Observe que para p = 0, 5, vale que Z(Θ0 ) = (1 + e0 )n = 2n , logo: Z(Θ) = 2n ·
t
EΘ0 eΘ ·η(G) .
Suponha que simulamos m de esses grafos: G1 , . . . Gm , pois é fácil simular de
PΘ0 (·) que são grafos de Erdös-Rényi, então, usando a lei de los grandes números, vale
que:
t t
eΘ ·η(G1 ) + · · · + eΘ ·η(Gm ) t
lim = EΘ0 eΘ ·η(G) .
m→∞ m
Em principio funciona, mas ocorrem outros problemas quando n é grande e deixa
de funcionar esta ideia.
Exemplo 1.5. Considere o exemplo anterior com Θ = (−1, 0.3) , com m = 100
grafos de Erdös-Rényi com p = 0.5, usando o seguinte código R para implementar as
simulações:
m=100
g.sim <- simulate(network(10,,density=0.1,directed=FALSE)
~ edges + triangle, coef=c(0, 0),nsim=m)
S<-print(g.sim,stats.print=TRUE)
eta<-attr(S,"stats")
theta=c(-1,0.3)
Theta<-matrix(rep(theta, m),nrow = m)
s<-c()
for(i in 1:m){
s[i]<-exp(Theta[i,]%*%eta[i,])
}
Z_Theta<-2^10*(sum(s)/m)
Z_Theta

Obtemos Z(Θ) = 2462486387

2. Estimadores de Máxima Verosimilhança (ERGM)

Após a construção de modelos de grafos aleatórios vem o problema, empírico,
de encontrar o valor dos parâmetros destes modelos a partir das observações, que
na Estatística se conhece como Inferência Paramétrica. Inicialmente vamos descrever o
paradigma clássico de estimação paramétrica que se conhece como ‘Método da Máxima
 2. ESTIMADORES DE MÁXIMA VEROSIMILHANÇA (ERGM) 37

Verosimilhança’, depois, discutiremos como aplicamos este método para o parâmetro θ

del modelo ERGM.
2.0.1. Estimadores de Máxima Verosimilhança. Suponha que temos m observa-
ções independentes: x1 , . . . xm de um modelo aleatório com distribuição Pθ , para esta
‘amostra’ definimos a verosimilhança como a função:
m
Y
L := L(θ, x1 , . . . xm ) = Pθ (xi )
i=1

O método da máxima verosimilhança consiste em encontrar um θ que maximize a

função L. θ também seria o máximo do logaritmo de L, então, definimos o estimador
de máxima verosimilhança (EMV) do parâmetro θ como o valor:
θ̂ = arg max log L(θ, x1 , . . . xm ) = arg max `(θ, X ).
θ∈Θ θ∈Θ ˜
Para distribuições (da família exponencial) que tem a forma, Pθ (x) = ceθ·T (x) ,
vale a seguinte propriedade:
Eθ̂ T (X) = T (X ) (2.6)
˜
2.0.2. Estimadores de Máxima Verosimilhança para ERGM. Para o modelo ERGM
temos uma única observação, o grafo Go , com verosimilhança dada por:
eθ·η(Go )
Pθ (Go ) = ,
Z(θ)
e o EMV de θ é
θ̂ = arg max `(θ, Go ).
θ∈Θ
O problema de avaliar `(θ, Go ) é que necessitamos o valor da função de partição
Z(θ), vamos a usar a mesma ideia de simulações para achar una expressão aproximada
para `(θ, Go ), e suponha que podemos simular da distribuição Pθ0 . Definimos r(θ, θ0 ) =
`(θ, Go ) − `(θ0 , Go ).
Observe que:
 θ·η(Go )   θ0 ·η(Go ) 
e e
`(θ, Go ) − `(θ0 , Go ) = log − log
Z(θ) Z(θ0 )
 
Z(θ)
= (θ − θ0 ) · η(Go ) − log .
Z(θ0 )
Como visto antes simulamos m grafos G1 , . . . Gm usando Pθ0 e utilizando a lei dos
grandes números podemos aproximar o valor de
  m
!
Z(θ) 1 X (θ−θ0 )·η(Gi )
log ≈ log e .
Z(θ0 ) m i=1
Finalmente„ usamos o valor aproximado para definir:
m
!
1 X (θ−θ0 )·η(Gi )
rm (θ, θ0 ) = (θ − θ0 ) · η(Go ) − log e .
m i=1
Assim, o EMV aproximado para o modelo é:
θ̂m = arg max rm (θ, θ0 ).
θ∈Θ
38 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

Como o modelo pertence à família exponencial, então, o EMV θ̂ satisfaz a relação (2.6),
Eθ̂ η(X) = η(Go ).
O programa ergm do R implementa esta solução. Veja o seguinte exemplo.

Exemplo 2.1. Usaremos como exemplo o grafo G2 simulado a partir de um ERGM

com estatísticas: o número de vértices η1 (G) = |E| e η1 (G) = |triângulos(G)| com
θ = (θ1 , θ2 ) = (1, 0.3):

Figura 9. Grafo G2 do modelo ERGM com θ1 = 1, θ2 = 0.3

Para este grafo a saída do programa é:

theta_est<-ergm(G2 ~ edges+ triangle,control=control.ergm(MCMC.burnin=100))
summary(theta_est)

==========================
Summary of model fit
==========================

Formula: G2 ~ edges + triangle

Iterations: 2 out of 20

Monte Carlo MLE Results:

Estimate Std. Error MCMC % p-value
edges 0.0197 2.4155 0 0.994
triangle 0.3822 0.3916 0 0.334
 3. ESTIMADORES BAYESIANOS (ERGM) 39

O valor do EMV é θ̂ = (0.0197, 0.3822). A estimação não foi boa (p-valor para θ1
é 0.994, muito alto, e p-valor para θ2 é 0.334, baixo, mais não aceitável, pois aceitamos
valores menores que 0.025), isto era esperado pois, estamos utilizando uma aproximação
bastante regular para a verosimilhança. Podemos ver que utilizando técnicas Bayesianas
podemos melhorar o resultado.
Avaliamos a estimação usando a distribuição de outras características: grau, arestas
comuns e geodésicas de um conjunto de grafos simulados com os parâmetros estimados.
gofG2 <- gof(theta_est)

par(mfrow=c(1,3))
par(oma=c(0.5,2,1,0.5))
plot(gofG2)
o resultado é dado por:

Figura 10. Bondade do ajuste (goodness of fit: gof) para o valor

estimado θ̂ = (0.0197, 0.3822)

Podemos observar na Figura (10) que estas outras características do grafo G2 (a

linha mais obscura dos gráficos correspondentes ao grau, arestas comuns e geodésicas)
se encontram dentro dos valores esperados das simulações (as linhas cinça dos gráficos).

3. Estimadores Bayesianos (ERGM)

A estadística Bayesiana se baseia no famoso teorema de Bayes:
P(B|A)P(A)
P(A|B) = .
P(B)
Supondo que o parâmetro que será estimado θ é aleatório e que sua distribuição
é P(θ), podemos expressar a verosimilhança como L(θ, X ) = P(X |θ) e reformular o
˜ ˜
problema original de achar θ ao problema de achar a distribuição de θ condicional à
informação observada X ,
˜
P(X |θ)P(θ)
˜
P(θ|X ) = .
˜ P(X )
˜
40 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

Nesta expressão π(θ) := P(θ) é chamada de distribuição a priori e, o objetivo

consiste em encontrar a distribuição a posteriori π(θ|X ) := P(θ|X ). Observe que
˜ ˜
P(X ) não depende de θ, logo é constante, assim, a relação de Bayes pode ser expressa
˜
como:

P(θ|X ) = cL(θ, X )π(θ),

˜ ˜

onde c é uma constante. Existem dois métodos para encontrar uma expressão para a
a posteriori. O primeiro consiste em encontrar uma distribuição que tenha a mesma
expressão que a posteriori, neste caso dizemos que ambas, a priori e a posteriori são
conjugadas, o que em general não ocorre.
O segundo método, mais general, consiste em encontrar um mecanismo que nos
permite simular m observações da posteriori: θ1 , . . . , θm , θi ∼ π(θ|X ). Dependendo
˜
do mecanismo de simulação (usualmente simulamos um processo estocástico) vale que:
m
1 X
π(θ|X ) ≈ Fm (θ|X ) := δθ (θ),
˜ ˜ m i=1 i

onde δa (x) = 1 se a = x e δa (x) = 0, caso contrario. A função Fm (·) é chamada

distribuição empírica. Usualmente a relação anterior é consequência da ergodicidade
do processo (isto é da unicidade da distribuição estacionaria).
Um método que nos permite simular os θ’s é o método de Metropolis-Hastings,
que consiste em usar uma cadeia de Markov que seja fácil de simular e, digamos, que
tenha probabilidade de transição q(θ0 |θ), então o objetivo consiste em perturbar esta
cadeia de modo tal que a nova cadeia tenha a distribuição estacionaria desejada, isto é,
a distribuição a posteriori. Suponha que a nova cadeia tem probabilidade de transição
k(θ0 |θ) = q(θ0 |θ)α(θ0 |θ). Se vale que:

k(θ0 |θ)π(θ|X ) = k(θ|θ0 )π(θ0 |X ),

˜ ˜

isto é, que a cadeia modificada, com probabilidade de transição k(θ0 |θ), é ‘reversível’
para π(θ|X ). Então, temos necessariamente que π(θ|X ) é a única distribuição estaci-
˜ ˜
onaria para k(θ0 |θ). Esta é uma propriedade das cadeias de Markov muito utilizada na
simulação de processos.
Usando as ideias anteriores o algoritmo de Metropolis-Hastings, que nos permite
simular os θ’s necessários para encontrar a distribuição empírica Fm (·), resume-se aos
seguintes passos:
Passo 1 Simulamos θ0 da distribuição a priori π(θ),
Passo 2 Supondo que tenhamos simulado θn , simulamos um candidato θ∗ usando
q(θ, θn ).
Passo 3 Avaliamos:

π(θ∗ |X )q(θ∗ , θn ) cL(θ∗ , X )π(θ∗ )q(θ∗ , θn )

   
˜ ˜
α(θ0 |θ) = min  , 1 = min  , 1 .
π(θn |X )q(θn , θ∗ ) cL(θn , X )π(θn )q(θn , θ∗ )
˜ ˜
 3. ESTIMADORES BAYESIANOS (ERGM) 41

Observe que a constante c cancela-se. Geralmente usamos uma cadeia tal

que q(u, u0 ) = q(u0 , u) pelo que a expressão anterior se simplifica:
L(θ∗ , X )π(θ∗ )
 
˜
α(θ0 |θ) = min  , 1
L(θn , X )π(θn )
˜
Passo 4 Simulamos uma observação u da distribuição Uniforme e fazemos:
(
θ∗ , se u ≤ α(θ0 |θ),
θn+1 =
θn , se u > α(θ0 |θ).
Passo 5 Volvemos ao passo 2.
Os passos são iterados até que os θ’s estejam ‘estacionários’. O algoritmo foi
chamado da ‘revolução computacional’ por Persi Diaconis, ver [8].
3.0.1. Estimadores Bayesianos (ERGM). Consideramos um grafo observado Go
que segue o modelo ERGM com verosimilhança:
eθ·η(Go )
L(θ, Go ) = Pθ (Go ) =
Z(θ)
e distribuição a priori π(θ) queremos encontrar a distribuição a posteriori π(θ|Go )
usando o algoritmo de Metropolis-Hastings. Suponha que os θ’s sejam gerados usando
uma cadeia de Markov simétrica: q(θ, θ0 ) = q(θ0 , θ). Então, para avaliar o α(θ0 |θ) do
algoritmo temos que avaliar:
∗
L(θ∗ , Go )π(θ∗ ) π(θ∗ )e(θ −θn )·η(Go ) Z(θn )
= ,
L(θn , Go )π(θn ) π(θn ) Z(θ∗ )
Como podemos ver a dificuldade de implementar isto, consiste em avaliar Z(θn )/Z(θ∗ )
e como foi visto, isso da origem a um error muito grande. Caimo et. al. ( 2012) mo-
dificaram este algoritmo, usando uma simulação auxiliar para θ, veremos no exemplo
este resultado.
Exemplo 3.1. Usaremos como exemplo o grafo G2 simulado a partir de um ERGM
com estadísticas: o número de vértices η1 (G) = |E| e η1 (G) = |triângulos(G)| com
θ = (θ1 , θ2 ) = (1, 0.3); em lugar de procurar θ̂, vamos a supor que, ‘a priori’, este valor
é representado por uma distribuição, a priori, π(θ1 , θ2 ) ∼ Uniforme em [0, 2] × [0, 1] e
queremos a distribuição empírica da a posteriori π((θ1 , θ2 )|Go ). Para isso usaremos o
programa bergm do R.
library("ergm")
library("bergm")
g.sim <- simulate(network(10,,density=0.1,directed=FALSE) ~ edges + triangle,
coef=c(0, 0))

G2 <- simulate( ~ edges + triangle, coef=c(1, 0.3),directed=FALSE,nsim=1,

basis=g.sim,seed=786,control=control.simulate(
MCMC.burnin=10000,
MCMC.interval=100))
theta_post <- bergm(G2 ~ edges + triangle,
burn.in=10,
aux.iters=500,
main.iters=500,
gamma=1)
42 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

bergm.output(theta_post)

e os resultados são:
Overall posterior density estimate:
theta1 (edges) theta2 (triangle)
Post. mean 3.377747 -0.02407964
Post. sd 5.224419 0.73015466

Overall acceptance rate: 0.31

Podemos observar que as medias da distribuição empírica Fm (θ1 , θ2 ) são,

(θ¯1 , θ¯2 ) = (3, 377747; −0, 02407964),

não muito diferentes dos valores originais, lembre que agora os θ’s são aleatórios. Ob-
serve que a probabilidade de aceitação de cada θ foi de 0,31.
A seguir pode ser observada a convergência das cadeias para a estacionaridade:

Figura 11. Distribuição a posteriori e convergência do processo (θ1 (n), θ2 (n))

4. Usando R
Neste capitulo utilizamos os programas do R: igraph, intergraph, network, ergm e
bergm.
(1) Figura 1:
G1 <- graph.formula(1,2,3)
G2 <- graph.formula(1--2,3)
G3 <- graph.formula(1,2--3)
G4 <- graph.formula(1--3,2)
G5 <- graph.formula(1--2--3)
G6 <- graph.formula(1--3--2)
 4. USANDO R 43

G7 <- graph.formula(2--1--3)
G8 <- graph.formula(1--2--3,3--1)
graphics.off()
par(mfrow=c(3,3))
plot(G1, xlab=’G1’,vertex.size=55,layout=layout_in_circle)
plot(G2, xlab=’G2’,vertex.size=55,layout=layout_in_circle)
plot(G3, xlab=’G3’,vertex.size=55,layout=layout_in_circle)
plot(G4, xlab=’G4’,vertex.size=55,layout=layout_in_circle)
plot(G5, xlab=’G5’,vertex.size=55,layout=layout_in_circle)
plot(G6, xlab=’G6’,vertex.size=55,layout=layout_in_circle)
plot(G7, xlab=’G7’,vertex.size=55,layout=layout_in_circle)
plot(G8, xlab=’G8’,vertex.size=55,layout=layout_in_circle)
plot(5,5,type="n",axes=FALSE,ann=FALSE,xlim=c(0,10),
ylim =c(0,10))
(2) Figura 2:
graphics.off()
par(mfrow=c(1,3))
plot(G5, xlab=’G5’,vertex.size=35,layout=layout_in_circle)
plot(G6, xlab=’G6’,vertex.size=35,layout=layout_in_circle)
plot(G7, xlab=’G7’,vertex.size=35,layout=layout_in_circle)
(3) Figura 3:
G<-erdos.renyi.game(30, 100, type ="gnm", directed = FALSE,
loops = FALSE)
plot(G)
(4) Figura 4:
G1<-erdos.renyi.game(10, 0.3, type ="gnp", directed = FALSE,
loops = FALSE,)
G2<-erdos.renyi.game(10, 0.5, type ="gnp", directed = FALSE,
loops = FALSE,)
G3<-erdos.renyi.game(10, 0.7, type ="gnp", directed = FALSE,
loops = FALSE)

graphics.off()
par(mfrow=c(1,3))
plot(G1,xlab=’G1’)
plot(G2,xlab=’G2’)
plot(G3,xlab=’G3’)
(5) Figuras 5 e 6:
graphics.off()
G1 <- sample_pa(500, power=0.1,directed = FALSE)
G2 <- sample_pa(500, power=0.001,directed = FALSE)
graphics.off()
par(mfrow=c(1,2))
plot(G1,vertex.label=NA,vertex.size=3)
plot(G2,vertex.label=NA,vertex.size=3)
graphics.off()
par(mfrow=c(1,2))
par(mfrow=c(1,2))
hist(degree(G1),main=’’)
hist(degree(G2),main=’’)
(6) Figuras 7 e 8:
44 3. GRAFOS ALEATÓRIOS: MODELOS E INFERÊNCIA

library("ergm")
g.sim <- simulate(network(10,,density=0.1,directed=FALSE)
~ edges + triangle, coef=c(0, 0))
summary(g.sim ~ edges + triangle)
G1 <- simulate( ~ edges + triangle, coef=c(-1, 0.3),
directed=FALSE,nsim=1,basis=g.sim,seed=234,
control=control.simulate(
MCMC.burnin=10000,
MCMC.interval=100))
summary(G1 ~ edges + triangle)
G2 <- simulate( ~ edges + triangle, coef=c(1, 0.3),
directed=FALSE,nsim=1,basis=g.sim,seed=786,
control=control.simulate(
MCMC.burnin=10000,
MCMC.interval=100))
summary(G2 ~ edges + triangle)
g1<-asIgraph(G1)
g2<-asIgraph(G2)

par(mfrow=c(1,2))
plot(g1,vertex.label=NA,vertex.size=25,xlab=’G1’)
plot(g2,vertex.label=NA,vertex.size=25,xlab=’G2’)
par(mfrow=c(1,2))
hist(degree(g1),main=’’, xlab=’grau(G1)’,
xlim = c(0,10),ylim = c(0,7))
hist(degree(g2),main=’’, xlab=’grau(G2)’,
xlim = c(0,10),ylim = c(0,7))
 CAPíTULO 4

Modelos de Regressão para Grafos e Modelos

Dinâmicos

Neste capítulo discutiremos alguns modelos de regressão para grafos: vizinhos

próximos, campos de Markov, espectrais e por Kernels. Um exemplo sobre genética
sera apresentado e ser utilizados parte dos códigos de R do libro [7].
Além disso discutiremos um modelo dinâmico para grafos, o modelo de epidemio-
logia SIR (susceptível, infectado e recuperado ), neste modelo veremos como influi na
velocidade de propagação da doença quando se considera a infeção propagando-se num
grafo.

1. Modelos de Regressão ou Previsão para Grafos

Estes modelos são preditivos no sentido de que uma variável aleatória Y é explicada,
ou prevista, por outra variável aleatória X, isto é, como função de X. Esta função
usualmente é linear, de aí o nome de ‘regressão linear’. Existem muitas variedades de
modelos de regressão e para o caso multilinear a variável X é trocada por um vetor,
usando um grafo podemos generalizar este modelo.
No caso de aplicações de regressão usando grafos, as variáveis estão indexadas com
os vértices do grafo, como num campo de vectores, a variável tem um valor em cada
vértice. A previsão consiste em que saber o valor da variável num vértice a partir dos
valores dos vértices conectados a ele.

1.1. Regressão Linear. No caso mais simples, não aleatório, temos um conjunto
de pontos (x1 , y1 ), . . . , (xn , yn ) no plano e queremos encontrar uma reta Y = α̂ + β̂X
cujos pontos minimizam a distância ao conjunto de pontos observados. α̂ é chamado
o interceptor e β̂ a inclinação.
Encontrar (α̂, β̂) corresponde a resolver o problema de minimização:

min {(Y − α − βX)t (Y − α − βX)},

(α,β)

onde X = (x1 , . . . , xn ) e Y = (y1 , . . . , yn ).

Como (Y − α + βX)t (Y − α − βX) = ||Y − α − βX||22 , chama-se a este problema
de minimização quadrática, este valorP é a distancia euclídea entre Y e α + βX. Lembre
que distância(a, b) = ||a − b||2 = i (ai − bi )2 , a, b ∈ Rn .
No caso aleatório, os pontos (x1 , y1 ), . . . , (xn , yn ) são observações do vector alea-
tório (X, Y ), cujas marginais X e Y tem uma relação de dependência: Y = α+βX +ε,
aqui, ε é uma variável aleatória com esperança E(ε) = 0 e variância V(ε) = σ 2 > 0,
logo E(Y |X = x) = α + βx por isso dizemos que este modelo é o melhor preditor lineal
para Y .
45
46 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

A solução deste problema consiste em minimizar:

min E{(Y − α − βX)t (Y − α − βX)},
(α,β)

Observe que ECM:= E{(Y − α − βX)t (Y − α − βX)} = E(Y − α − βX)2 , o erro

quadrático médio é a quantia a ser minimizada.
1.2. Regressão Vizinhos Próximos. A regressão por vizinhos próximos, NN (Ne-
arest Neirghbor), é um tipo de regressão local. Suponha que temos como dado um
grafo G = (V, E) e que em cada vértice i ∈ V temos uma variável Xi com a mesma
distribuição . Mais desconhecemos a relação de dependência entre estas variáveis,
sabemos que estão ligadas pelo grafo G.
Se temos observações {xi , i ∈ V } nosso objetivo é predizer o valor de cada xi a
partir dos valores de x no seus vizinhos Ni = {j ∈ V : (i, j) ∈ E}, o preditor NN para
xi neste caso está definido por:
P
j∈Ni xj
x̂i = .
|Ni |
Exemplo 1.1. Vamos a usar um exemplo da genética, os dados consistem de
uma rede de 241 interações (arestas) entre 134 proteínas (vértices) pertencentes à um
organismo, Saccharomyces cerevisiae, cada proteína tem um atributo ((ICSC)i = 1 o
0) que indica se a ‘cascada de sinalização intracelular’ (intracellular signaling cascade)
está presente o não, esta é uma característica de comunicação intracelular. Estos dados
foram colectados por [6]. Este grafo é chamado ppi.CC,

Figura 1. Grafo ppi.CC de interações entre proteínas de Saccha-

romyces cerevisiae: ICSC = 1 (amarelo) e ICSC = 0 (azul)

Se usamos o grafo ppi.CC e a relação de vizinhos próximos descrita antes obtemos

um grafo ppi.CC.nn igual ao anterior pero com valores de ICSC previstos por esta
relação NN, podemos comparar os dois gráficos:
 1. MODELOS DE REGRESSÃO OU PREVISÃO PARA GRAFOS 47

Figura 2. Grafos ppi.CC e ppi.CC.nn

Para verificar que tanto coincidem os valores dos vértices vamos a contar as coin-
cidências em cada vértice:

Figura 3. Coincidências entre os valores dos vértices ppi.CC e ppi.CC.nn

48 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

Observe que valores próximos a 1 (cor amarelo) tem uma frequência alta de coinci-
dência, isto é, poucos dos vértices que tem valor perto de zero são amarelos em ppi.CC.
(o mesmo acontece para (ICSC)i = 0 em ppi.CC). Vejamos qual a porcentagem de
mal classificados:
#### %de error: não coincidencias
nn.pred <- as.numeric(nn.ave > 0.5)
mean(as.numeric(nn.pred != V(ppi.CC.gc)$ICSC))
[1] 0.2598425

O erro de predição é 25,9 %, que é baixo.

1.3. Regressão Linear para Campos de Markov. Nesta subseção vamos definir
um modelo aleatório para as observações xi dos vértices, este modelo é chamado
Campo Aleatório Marcoviano .
Seja G = (E, V ) um grafo com |V | = n e X = (X1 , . . . , Xn ) um conjunto de
variáveis aleatórias definidas em V , dizemos que X é um Campo Aleatório Markoviano
(MRF: Markov Random Field em inglês) se:
(1) P(X = x) > 0 para todo x ∈ Rn .
(2) É valida a propriedade de Markov:
P(Xi = xi |X(−i) = x(−i) ) = P(Xi = xi |XNi = xNi )
com X(−i) = (X1 , . . . , Xi−1 , Xi+1 , . . . , Xn ) e XNi é X restrito aos vértices
j ∈ Ni = {j ∈ V : (i, j) ∈ E}
A propriedade (2) diz que a distribuição de Xi , depende somente dos vizinhos do
vértice i, isto é, ele é Markov, esta distribuição que depende somente dos vizinhos é
chamada de ‘especificação local’ do modelo.
O modelo MRF é comum na estatística espacial ou nos modelos para o proces-
samento digital de imagens, origina-se no trabalho de Besag (1976), ver [1]. Uma
propriedade dos MRF, o teorema de Hammersley–Clifford, que permite representar X
como:
eU (x)
P(X = x) = (1.7)
Z P U (x)
onde U (x) é uma função real, chamada ‘Potencial’, e Z = x e é a função de
partição, podemos ver que este modelo inclui os modelos exponenciais. Através desta
propriedade é possível que posamos fazer modelos de regressão lineal. Distribuições do
tipo (1.7) são chamadas medidas de Gibbs em honor ao físico J. W. Gibbs que as
descreveu no seculo XIX. P
É possível mostrar que U (x) = c∈C Uc (x) onde C é o conjunto de cliques ou
grafos completos Kn contidos em G. Por exemplo, um vértice é K1 ( o grafo completo
com 0 arestas e um vértice), uma aresta é K2 ( o grafo completo com 1 aresta e 2
vértices), um triângulo é K3 ( o grafo completo com 3 arestas e 3 vértices), etc. Em
consequência a equação (1.7) diz que a distribuição de X depende dos valores de U (x)
nestes ‘cliques’.
 1. MODELOS DE REGRESSÃO OU PREVISÃO PARA GRAFOS 49

O modelo MRF que vamos utilizar é chamado de Modelo Auto-Logístico e cor-

responde a um potencial que somente depende dos valores nos vértices (K1 ) e nas
arestas (K2 ), além do mais, os valores de xi podem ser só 1 o 0, isto é,

X X
U (x) = αxi + βxi xj .
i∈V (i,j)∈E

È possível provar que nesse caso a especificação local do Modelo Auto-logístico é

dada por:

P
eα+β j∈Ni xj
Pα,β (Xi = 1|XNi = xNi ) = P .
1 + eα+β j∈Ni xj

Então, vale que:

Pα,β (Xi = 1|XNi = xNi ) X

log =α+β xj .
Pα,β (Xi = 0|XNi = xNi )
j∈Ni

Assim é possível fazer regressão linear para as especificações locais.

Exemplo 1.2. Continuando com o exemplo das interações com as proteínas,

vamos usar o modelo auto-logístico para predizer os valores de ICSC no grafo ppi.CC
e logo comparar os resultados com o grafo realmente observado; primeiro, temos que
encontrar α e β, para isto usamos o programa ngspatial de estatística espacial do R ,
usando o comando autologistic obtemos:

m1.mrf <- autologistic(formula1, A=A,control=list(confint="none"))

m1.mrf$coefficients
(Intercept) eta
0.2004949 1.1351942

Logo, os coeficientes são α = 0, 2 e β = 1, 13 uma comparação gráfica dos grafos:

50 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

Figura 4. Grafos ppi.CC e mrf 1

como antes, vamos verificar quantos vértices coincidem:

Figura 5. Coincidências entre os valores dos vértices ppi.CC e mrf 1

A porcentagem de não coincidências neste caso é de 20%, que é mais baixa que
do modelo de vizinhos próximos.
 1. MODELOS DE REGRESSÃO OU PREVISÃO PARA GRAFOS 51

> mean(as.numeric(mrf1.pred != V(ppi.CC.gc)$ICSC))

[1] 0.2047244

Finalmente, como em todos os modelos de regressão linear podemos pesquisar o

que tanto é explicado pelo modelo:
####explicação do modelo
> assortativity(ppi.CC.gc, X+1, directed=FALSE)
[1] 0.3739348

Os parâmetros do modelo explicam 37,39% da informação contida nos dados do

grafo ppi.CC.gc.
1.4. Regressão por Kernels para Grafos. Os modelos de predição que usam
Kernels são conhecidos nas pesquisas sobre modelos de auto-aprendizado ou Machine
Leraning, são regressões locais, usualmente se usa uma rede para estimar os parâmetros
locais e logo se faz predição do comportamento da nova informação. Um Kernel
K = (ki,j ){i,j∈V } é uma matriz simétrica semi-definida positiva, isto é, ki,j = kj,i e
xt Kx ≥ 0 para todo x ∈ R|V | ; além disso satisfaz para todo subconjunto de m vértices
V (m) = {vi1 , . . . , vim } a restrição de K a esse conjunto K (m) = (ki,j ){i,j∈V (m) } , ela
também é simétrica e semi-definida positiva.
Uma propriedade das matrizes definidas positivas (respetivamente semi-definidas
positivas), um Kernel é semi-definida positiva, é que elas definem uma distância (res-
petivamente quase-distância) entre os vectores:
d(x, y) = (x − y)t K(x − y),
esta distância é uma norma, pois K é simétrica.
Isto nos permite considerar o problema de regressão local:

min{(X − K (m) α)t (X − K (m) α) + λαt (K (m) )−1 α}. (1.8)

α
Para usar esta regressão em um grafo G, vamos usar o Laplaciano L = D − A do
grafo que é uma matriz simétrica e semi-definida positiva , onde A é a matriz de
Adjacência e D = diag(d(i)). Fazemos K = L−1 (esta é uma pseudo inversa pois
existem autovalores zero) que também é simétrica y semi-definida positiva. Vale a
‘decomposição espectral’:
K = Φt ∆Φt ,
donde Φ esta formado pelos auto-vetores de K e ∆ = diag(δi ) são os auto-valores de
K = L−1 , isto é, γi = δi−1 são os auto-valores de L o Laplaciano de G ( se existem
auto-valores zero, fazemosδi := 0):
X
L−1 = f (γi )φt φ,
i∈V

com f (x) = 1/x, x 6= 0 e f (x) = 0, x = 0.

Usando esta decomposição podemos mostrar que,
X
αt K −1 α = ht Lh = (h1 − hj )2 ,
(i,j)∈E

onde h = Φα é uma combinação linear dos auto-vetores. Logo na minimização serão

favorecidos vértices que sejam adjacentes.
52 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

Exemplo 1.3. Vamos usar a regressão local para predizer os valores de ICSC no
grafo ppi.CC . Primeiro, vejamos a distribuição dos f (γi ):

Figura 6. Distribuição dos f (γi )

Podemos ver o efeito de selecionar os auto-vetores maiores em cada vértice do

grafo:

Figura 7. Auto-vetores maiores no radio de cada vértice do grafo

Finalmente o preditor usando kernels no grafo ppi.CC :

 1. MODELOS DE REGRESSÃO OU PREVISÃO PARA GRAFOS 53

Figura 8. Grafos ppi.CC.gc e knl

a comparação entre os dos:

Figura 9. Coincidências entre ppi.CC.gc e knl

54 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

a porcentagem de errores foi 11%

> mean(as.numeric(m1.svm.fitted != V(ppi.CC.gc)$ICSC))
[1] 0.1102362

2. Outros Tópicos
Entre inúmeros tópicos sobre aplicações em grafos aleatórios podemos mencionar
a análise de fenômenos como os grafos de escala-libre, como é o comportamento dos
diferentes modelos de previsão, nos exemplos só temos considerado grafos finitos e
relativamente pequenos, embora que é conhecido que a medida que os grafos aumentam
de tamanho, fenômenos tais como: degeneração, transição de fase o criticalidade são
relevantes, isto é, algumas das características mudam completamente.
Um dos tópicos que discutiremos é o efeito do uso de um grafo aleatório na disper-
são de alguma doença. A continuação discutimos o modelo de infeção SIR. A seguir
vamos usar o exemplo e os códigos de R do libro de [7] para este modelo.

2.1. Modelos Dinâmicos. Um dos modelos padrão da Epidemiologia é o processo

SIR, Susceptíveis-Infectados-Recuperados. Neste processo observamos ao longo do
tempo a evolução de uma doença, primeiro, o número de infectados na população
cresce, para logo ir decrescendo a medida que vão ficando menos pessoas susceptíveis
de serem infectadas pois se tem imunizado depois da recuperação, como está descrito
esquematicamente no seguinte gráfico:

Figura 10. Evolução de uma infecção (vermelho), gráfico tomado

de [7]

Para construir o processo, que descreve a evolução no gráfico, considere uma po-
pulação de N + 1 indivíduos, e digamos que em qualquer instante essa população
decompõe-se em três tipos: Susceptíveis (S), Infectados (I) e Recuperados (R), assim
teríamos que: N + 1 = NS (t) + NI (t) + NR (t), onde NS (t) (resp.NI (t) NR (t)) é o
número de Susceptíveis (resp. Infectados, Recuperados) no instante t.
No instante inicial NI (0) = 1 e NS (0) = N , a evolução subsequente é governada
por uma cadeia de Markov à tempo contínuo, isto é, por um processo estocástico que
tem mudanças somente em certos tempos aleatórios, as mudanças estão governados
 2. OUTROS TÓPICOS 55

pelas seguintes probabilidades de transição:

P(NS (t + δt) = s − 1, NI (t + δt) = i + 1|NS (t) = s, NI (t) = i) ≈ βsiδt,
P(NS (t + δt) = s, NI (t + δt) = i − 1|NS (t) = s, NI (t) = i) ≈ γiδt e
P(NS (t + δt) = s, NI (t + δt) = i|NS (t) = s, NI (t) = i) ≈ 1 − (βs + γ)iδt.
(2.9)
onde δt significa que esta transição acontece em um tempo infinitesimal depois de t, β
é a taxa de infeção e γ a taxa de recuperação. O comportamento deste processo esta
descrito na Figura (10).
O problema do modelo assim descrito é que as interações espaciais não existem,
o que significa que um infectado pode contagiar a qualquer susceptível, o que não é
realista. Uma maneira de melhorar o modelo é introduzir efeitos espaciais por meio
de um grafo G = (E, V ) que descreva as relações entre os indivíduos. Para isso,
consideramos o processo Xi (t) = 0, 1, 2, com i ∈ V , onde 0 (respectivamente 1, 2) se
o vértice i esta susceptível (respectivamente infectado, recuperado). O processo Xi (t)
obedece:

βMi (x)δt
 se xi = 0 e x0i = 1
0
P(X(t + δt) = x |X(t) = x) = γδt se xi = 1 e x0i = 2 (2.10)
1 − [βMi (x) + γ]δt se xi = 2 e x0i = 2



onde Mi (x) é o número de vizinhos do vértice i que estão infectados xj = 1, j ∈ Ni .

Finalmente NS (t), NI (t) e NR (t) são a soma dos vértices com valores 0, 1 ou 2.
56 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

Simulamos este processo usando o programa SIR do R, com β = 0.5 e γ = 1 e com

grafos simulados dos modelos de Erdös-Rényi, Barabasi-Albert (BA) e Watts-Strogatz
(este é outro modelo de escala-livre), para obter:

Figura 11. Evolução de uma infeção num grafo aleatório, exemplo

tomado de [7]

3. Usando R
Utilizamos neste capitulo os programas do R: sand, ngspatial, kernlab e sir.
(1) Figura 1:
###datos
set.seed(42)
library(sand)
data(ppi.CC)
 3. USANDO R 57

ppi.CC<-upgrade_graph(ppi.CC)
summary(ppi.CC)
### 10 primeros vertices
V(ppi.CC)$ICSC[1:10]
###visualizacion
V(ppi.CC)[ICSC == 1]$color <- "yellow"
V(ppi.CC)[ICSC == 0]$color <- "blue"
graphics.off()
plot(ppi.CC, vertex.size=5, vertex.label=NA)
(2) Figura 2:
####vizinhos proximos
clu <- clusters(ppi.CC)
ppi.CC.gc <- induced.subgraph(ppi.CC,
clu$membership==which.max(clu$csize))
nn.ave <- sapply(V(ppi.CC.gc),
function(x) mean(V(ppi.CC.gc)[nei(x)]$ICSC))
class(nn.ave)
###grafo predecido
nn<-ppi.CC.gc
nn$ICSC<-nn.ave
V(nn)[ICSC == 1]$color <- "yellow"
V(nn)[ICSC == 0]$color <- "blue"
graphics.off()
par(mfrow=c(1,2))
plot(ppi.CC, vertex.size=5, vertex.label=NA,xlab="Grafo ppi.CC ")
plot(nn, vertex.size=5, vertex.label=NA, xlab="Grafo ppi.CC.nn ")

(3) Figura 3:
#### porcentage de coincidencias
par(mfrow=c(2,1))
hist(nn.ave[V(ppi.CC.gc)$ICSC == 1], col="yellow",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices com ICSC")
hist(nn.ave[V(ppi.CC.gc)$ICSC == 0], col="blue",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices sin ICSC")

(4) Figura 4:
##### campos markovianos
set.seed(42)
library(sand)
data(ppi.CC)
#### grafo G aqui corresponde a ppi.CC
ppi.CC<-upgrade_graph(ppi.CC)
V(ppi.CC)[ICSC == 1]$color <- "yellow"
V(ppi.CC)[ICSC == 0]$color <- "blue"
library(ngspatial)

clu <- clusters(ppi.CC)

ppi.CC.gc <- induced.subgraph(ppi.CC,
clu$membership==which.max(clu$csize))
X <- V(ppi.CC.gc)$ICSC
58 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

A <- get.adjacency(ppi.CC.gc, sparse=FALSE)

formula1 <- X ~ 1
#### Modelo autologistico
m1.mrf <- autologistic(formula1, A=A,control=list(confint="none"))
m1.mrf$coefficients

#####prediccion
mrf1.pred <- as.numeric((m1.mrf$fitted.values > 0.5))
####comparacion
mrf1<-ppi.CC.gc
mrf1$ICSC<-mrf1.pred
V(mrf1)[ICSC == 1]$color <- "yellow"
V(mrf1)[ICSC == 0]$color <- "blue"
graphics.off()
par(mfrow=c(1,2))
plot(ppi.CC, vertex.size=5, vertex.label=NA,xlab="Grafo ppi.CC ")
plot(mrf1, vertex.size=5, vertex.label=NA, xlab="Grafo mrf1 ")

(5) Figura 5:
par(mfrow=c(2,1))
hist(mrf1.pred[V(ppi.CC.gc)$ICSC == 1], col="yellow",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices com ICSC")
hist(mrf1.pred[V(ppi.CC.gc)$ICSC == 0], col="blue",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices sin ICSC")

(6) Figura 6:
###### autovectores mayores
e.vec1 <- e.L$vectors[, (nv-1)]
v1.colors <- character(nv)
v1.colors[e.vec1 >= 0] <- "red"
v1.colors[e.vec1 < 0] <- "blue"
v1.size <- 15 * sqrt(abs(e.vec1))
l1 <- layout.fruchterman.reingold(ppi.CC.gc)
#####
e.vec2 <- e.L$vectors[, (nv-2)]
v2.colors <- character(nv)
v2.colors[e.vec2 >= 0] <- "red"
v2.colors[e.vec2 < 0] <- "blue"
v2.size <- 15 * sqrt(abs(e.vec2))
l2 <- layout.fruchterman.reingold(ppi.CC.gc)
plot(ppi.CC.gc, layout=l2, vertex.color=v2.colors,
vertex.size=v2.size, vertex.label=NA, xlab=c("(n-1) mayor autovector"))
#####
e.vec3 <- e.L$vectors[, (nv-3)]
v3.colors <- character(nv)
v3.colors[e.vec3 >= 0] <- "red"
v3.colors[e.vec3 < 0] <- "blue"
v3.size <- 15 * sqrt(abs(e.vec3))
l3 <- layout.fruchterman.reingold(ppi.CC.gc)
#####
 3. USANDO R 59

graphics.off()
par(mfrow=c(1,3))
plot(ppi.CC.gc, layout=l1, vertex.color=v1.colors,
vertex.size=v1.size, vertex.label=NA, xlab=c("(n) mayor autovector"))
plot(ppi.CC.gc, layout=l2, vertex.color=v2.colors,
vertex.size=v2.size, vertex.label=NA, xlab=c("(n-1) mayor autovector"))
plot(ppi.CC.gc, layout=l3, vertex.color=v3.colors,
vertex.size=v3.size, vertex.label=NA, xlab=c("(n-2) mayor autovector"))

(7) Figura 7:
library(kernlab)
K1.tmp <- e.L$vectors %*% diag(f.e.vals) %*%
t(e.L$vectors)
K1 <- as.kernelMatrix(K1.tmp)
m1.svm <- ksvm(K1, X, type="C-svc")
m1.svm.fitted <- fitted(m1.svm)
####comparación
knl<-ppi.CC.gc
knl$ICSC<-m1.svm.fitted
V(knl)[ICSC == 1]$color <- "yellow"
V(knl)[ICSC == 0]$color <- "blue"
graphics.off()
par(mfrow=c(1,2))
plot(ppi.CC, vertex.size=5, vertex.label=NA,xlab="Grafo ppi.CC ")
plot(knl, vertex.size=5, vertex.label=NA, xlab="Grafo knl ")

(8) Figura 8:
par(mfrow=c(2,1))
hist(m1.svm.fitted[V(ppi.CC.gc)$ICSC == 1], col="yellow",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices com ICSC")
hist(m1.svm.fitted[V(ppi.CC.gc)$ICSC == 0], col="blue",
ylim=c(0, 30), xlab="Proporcion de Vecinos com ICSC",
main="Vertices sin ICSC")

(9) Figura 11:

#########modelo SIR com grafos
####simulacion de grafos
library(sir)
gl <-list()
gl$ba<- barabasi.game(250, m=5, directed=FALSE)
gl$er<- erdos.renyi.game(250, 1250, type=c("gnm"))
gl$ws<- watts.strogatz.game(1, 100, 12, 0.01)
######tasas de infeccion e recuperacion
beta <- 0.5
gamma <- 1
#### numero de simulaciones
ntrials <- 100
######simulacion de SIR
sim <- lapply(gl, sir, beta=beta, gamma=gamma,
no.sim=ntrials)
#####porcentajes de infectados
60 4. MODELOS DE REGRESSÃO PARA GRAFOS E MODELOS DINÂMICOS

graphics.off()
par(mfrow=c(3,1))
plot(sim$er, xlab="tiempo", main="Erdos")
plot(sim$ba, color="palegoldenrod",
median_color="gold", quantile_color="gold",
xlab="tiempo", main="Barabasi")
plot(sim$ws, color="pink", median_color="red",
quantile_color="red",xlab="tiempo",
main="Watts-Strogatz")
 Referências

[1] J. Besag. Spatial Interaction and the Statistical Analysis of Lattice Systems. Journal
of the Royal Statistical Society Series B (Methodological), 36, 1974.
[2] S. Chatterjee and P. Diaconis. Estimating and understanding exponential random
graph models. Ann. Statist., 41(5):2428–2461, 2013.
[3] R. Diestel. Graph theory. Graduate Texts in Mathematics. Springer, 3rd edition,
2006.
[4] O. Frank and D. Strauss. Markov Graphs. Journal of the American Statistical
Association, 81, 09 1986.
[5] P. W. Holland and S. Leinhardt. A Method for Detecting Structure in Sociometric
Data. American Journal of Sociology, 76, 11 1970.
[6] X. Jiang, E. D. Kolaczyk, N. Nariai, M. Steffen, and S. Kasif. Integration of
relational and hierarchical network information for protein function prediction. BMC
Bioinformatics, 9, 12 2008.
[7] E. D. Kolaczyk. Statistical Analysis of Network Data: Methods and Models. Sprin-
ger Series in Statistics. Springer, 2009.
[8] P. D. L. Saloff-Coste. What Do We Know about the Metropolis Algorithm? Journal
of Computer and System Sciences, 57, 1998.

61