Olhando e comparando outros modelos para descrever propriedades de redes de interações

01

Continuando o que estávamos vendo no post anterior, vamos olhar novamente somente a rede do exemplo 4, acima.
Aqui temos bastante interações e desde polinizadores muito generalistas a poucos generalistas.
Apesar de essa rede ter sido gerada por acaso, é muito interessante o fato que as redes biológicas vão ser fruto das estratégias das espécies que a compõe, e ao conseguir descrever essa rede e falar ela é assim, podemos passar a pensar quais forças evolutivas estão agindo como um todo na comunidade. Claro que a pressão evolutiva vai ser exercida nos indivíduos resultando no que é o jeitão comum da espécie e não vemos esse processo, mas aqui estamos olhando para o resultado dele.

Então continuando, nos ajustamos o modelo exponencial para a “escadinha” que se forma ao contarmos quantas espécies tem um determinado número de interações.

Lembrando que esse é a função que a gente fala o degrau da escadinhas e a função da a altura do degrau:

 p(k) = b \cdot exp(-\gamma \cdot k)

E ajustando o modelo no R usando a função nls temos como resultado

Formula: P ~ b * exp(-gamma * k) Parameters: Estimate Std. Error t value Pr(>|t|) gamma 0.14839 0.01197 12.39 2.15e-07 *** b 1.26079 0.06602 19.10 3.37e-09 *** ---

E olhando esses números no gráfico temos

02

Temos nossos 2 parâmetros, que o mais interessante é o gamma (lembre-se que gamma é essa letra grega \gamma), que é a inclinação da escadinha. Mais interessante no sentido que ele esta informando muita coisa sobre a comunidade, ele que esta informando se temos mais especialistas, menos especialistas, mais generalistas, quantos de cada tipo, isso fica impresso naquele gráfico.

Existem ainda duas possibilidades como proposto pelo Jordano (2003).
A primeira é a power law, aqui o power é de potencia, um número elevado ao outro, que vemos na formula da função que é o modelo.

 p(k) = b \cdot k^{-\gamma}

Novamente ajustando no R temos:

Formula: P ~ b * k^(-gamma) Parameters: Estimate Std. Error t value Pr(>|t|) gamma 0.51585 0.09301 5.546 0.000245 *** b 1.18185 0.14298 8.266 8.83e-06 *** ---

E olhando no gráfico temos

03

Uma coisa legal de se pensar em relação a esse power law é que ele não enverga muito bem, ele tem uma capacidade de previsão muito pior que o modelo exponencial aqui nesse exemplo, mas talvez ele seja muito mais adequado para descrever vários generalistas, ou uma predominância de generalismo, já que não “envergar” bem o modelo quer dizer que os degraus (lembre-se do gráfico de barras do post anterior) vão continuar do mesmo tamanho, ou seja o mesmo número de espécies varias vezes com o mesmo número de interações, generalistas.

Agora vamos olhar a terceira opção, Truncated Power Law.
Bem se é power law, é algo a potência do outro, agora com esse “Truncated” que eu não sei o que faz direito. Mas vamos olhar ele melhor.

 p(k) = b \cdot (k^{-\gamma}) \cdot exp(-k/kx)

E temos um novo parâmetro, que deve ser o tal do “Truncated”, e vamos la e ajustamos ele no R e temos…

Formula: P ~ b * (k^(-gamma)) * exp(-k/kx) Parameters: Estimate Std. Error t value Pr(>|t|) gamma -0.29097 0.08258 -3.523 0.00648 ** b 1.22771 0.04615 26.605 7.24e-10 *** kx 4.48568 0.47416 9.460 5.67e-06 *** ---

Agora temos um gamma negativo, mas olhando o gráfico que é mais simples de ver as coisas que esses números, de inicio vemos que…

04

Ele parece que tem um ajuste melhor que todos os outros modelos, a linha passa bem pertinho dos pontos, mas por outro lado ele tem um parâmetro a mais. Parâmetros deixam a gente mexer na forma da linha logo, quanto mais parâmetros sempre a gente vai ter a possibilidade de um melhor ajuste, lembrando que possibilidade num é garantia.

Mas a gente pode usar o AIC para verificar a qualidade dos ajustes, e colocar em número isso que estamos olhando nesses gráficos.
Vamos fazer então um gráfico com os 3 modelos e olhar os respectivos AICs

05

#Exponencial AIC(m1.exe4) [1] -27.11152 #Power Law (m2.exe4) [1] -5.733472 #Truncated Power Law AIC(m3.exe4) [1] -36.24854

Bem o aic é uma medida de ajuste do modelo, o quanto a linha fica bem pertinho dos pontinhos do gráfico mais ou menos. Quanto menor é melhor, e ao olhar-mos aqui vemos o seguinte, olhando no gráfico da para ver facilmente que o Power Law é o pior nesse caso, talvez com mais generalistas ele fosse melhor, mas para esses dados que temos aqui ele é ruim, ou seja tem o Maior AIC, -5.

Agora olhando no gráfico, o modelo Exponencial e o Truncated Power Law tem um bom ajuste, e tem valores baixos de AIC, -27 e -36 respectivamente. No entanto podemos pensar um pouco ainda sobre simplesmente usar a regra de pegar o menor, o modelo truncado.
Bem o modelo truncado tem um parâmetro a mais. No exponencial as coisas são mais simples de imaginar, o b vai acabar empurrando a linha mais para cima ou mais para baixo no gráfico, e o gamma vai ser o quão íngreme é a subida da curva. Agora esse modelo truncando tem um outro parâmetro que deixa tudo mais difícil de pensar, ele parece que torce um pouco a curva, uma coisa a se fazer é ficar mudando valores de kx, b e gamma e tentar entender como tudo funciona se você como eu não conseguir imaginar como a curva vai ficar.

Outra coisa é que qualquer um ponto a mais no gráfico muda esse valor de AIC, qualquer interação que vimos a mais, ou que não conseguimos ver pode trocar quem é o modelo com ajuste melhor. Apesar da diferença dos dois modelos para o power law serem grandes, entre eles talvez esse 10 de diferença de AIC seja espurico. Mas podemos tentar fazer um grau de confiança se realmente esses modelos são diferentes, pois se ele não forem diferentes é mais fácil a gente ficar com o mais simples, que é mais fácil de pensar, sempre lembramos da navalha de Occam nesses momentos. Usando a estatística f vemos o seguinte ao comparar os modelos.

anova(m1.exe4,m3.exe4,test="f") Analysis of Variance Table Model 1: P ~ b * exp(-gamma * k) Model 2: P ~ b * (k^(-gamma)) * exp(-k/kx) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 10 0.044500 2 9 0.017591 1 0.026909 13.767 0.004842 ** ---

E parece que os modelos são diferentes, o modelo truncando é um melhor preditor e ficamos com ele. Mas esse teste f para comparar os modelos não é muito perfeito. Principalmente porque comparar modelo não linear é muito complicado. Nos fóruns as pessoas são bastante relutantes em usar essas coisas com modelos não lineares, como os que estamos ajustando aqui, mas a principio não estamos vendo nada que não esparavamos ao avaliar visualmente os gráficos e a analise dar um resultado que você já espera ao olhar o gráfico é uma coisa boa.

Então assim finalmente resolvemos usar o modelo truncado para descrever essa comunidade. Se tivessemos mais comunidades poderíamos comparar os parâmetros delas agora quanto a essa propriedade descrita por esse gamma do modelo e começar a pensar o que esta acontecendo, porque que esses valores fiquem assim, ou se tiver alguma diferença, o que poderia gerar essa diferença em comunidades.

Bem um ultimo detalhe é que se olhar o gráfico do bipartite ele fica diferente do que fizemos aqui.

06

Talvez para uma comunidade muito maior que a desse exemplo, o gráfico ficaria muito grande do jeito que eu fiz, então ele usa a escala logarítmica para plotar o gráfico. Quer dizer que ele calcula o log do valor de x e do valor de y para plotar. Eu achei mais ilustrativo usar os dados brutos para explorar aqui o que esta acontecendo, mas se a gente fizer questão de ver a mesma coisa é so colocar log nos valores de x e y

07

#Semente para os resultados serem reproduzíveis
set.seed(13)
###################################################
#Função para gerar rede de interação baseado no modelo de Erdős–Rényi
#Creditos a Edwin Grappin, vi o post no R-bloggers
#http://probaperception.blogspot.com.br/2012/12/have-you-tried-to-understand-your.html
###################################################
generateER = function(n = 100, p = 0.5){
map = diag(rep(1, n))
link = rbinom(n*(n-1)/2, 1,p)
t = 1
for(j in 2:n){
for(i in 1:(j-1)){
map[i,j] = link[t]
t = t + 1
}
}
dimnames(map)<-list(paste("Pol",1:n),paste("Sp",1:n))
return(map)
}
###########################
#Minha função para transformar matriz em lista para o igraph
##############################
mat2lis<-function(matriz) {
lista<-dimnames(matriz)
dados<-data.frame(Pol=rep(lista[[1]],length(lista[[1]])),
SP=rep(lista[[2]],each=length(lista[[2]])),
int=c(matriz))
return(dados[dados$int==1,1:2])
}
############################
library(igraph)
 
#Gerando exemplos
exe4<-generateER(n=15,p=0.8)
 
#Transformando na forma de dados que o igraph usa
exe4.l<-mat2lis(exe4)
exe4.g <-graph.data.frame(exe4.l, directed=F)
 
#Plots da rede
#01
plot(exe4.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.5,
vertex.label.color="black"
)
title("Exemplo 4")
 
#Para ver os modelos consulte…
#Jordano, P., Bascompte, J. and Olesen, J. M. (2003)
#Invariant properties in coevolutionary networks of
#plant-animal interactions. Ecology Letters 6, 69–81
 
#Exemplo 4
k<-sum(exe4)
pol.n<-rowSums(exe4)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol4 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe4<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol4)
m2.exe4<-nls(P~b*k^(-gamma),start=list(gamma = 0, b = 1),data=dados.pol4)
m3.exe4<-nls(P~b*(k^(-gamma))*exp(-k/kx),start=list(gamma = 0, b = 1,kx=5),data=dados.pol4)
 
#Grafico modelo exponencial
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red")
legend("topright",lwd=2,lty=2,col=c("red"),
legend=c("Exponential"))
#grafico do modelo power law
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue")
legend("topright",lwd=2,lty=2,col=c("blue"),
legend=c("Power Law"))
#grafico do modelo truncado
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green")
legend("topright",lwd=2,lty=2,col=c("green"),
legend=c("Truncated Power Law"))
 
#os 3 modelos juntos
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red")
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue")
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green")
legend("topright",lwd=2,lty=2,col=c("red","blue","green"),
legend=c("Exponential","Power Law","Truncated Power Law"))
 
#resultados dos modelos e seus AICs
summary(m1.exe4)
summary(m2.exe4)
summary(m3.exe4)
AIC(m1.exe4)
AIC(m2.exe4)
AIC(m3.exe4)
 
#comparação dos modelos, leia nos foruns primeiros
#sobre comparar modelos não lineares antes de confiar
#nesse resultado
anova(m1.exe4,m3.exe4,test="f")
 
#Grafico e ajuste usando o bipartite para o exemplo 4
library(bipartite)
degreedistr(exe4, plot.it=T)
 
#grafico como apresentado no bipartite
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,log="xy")
abline(h=1,lty=2)
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red",log="xy")
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue",log="xy")
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green",log="xy")
legend(2,0.3,lwd=2,lty=2,col=c("red","blue","green"),
legend=c("Exponential","Power Law","Truncated Power Law"))
 
#########################################

Propriedades em rede de interações, como descrever as comunidades de plantas e seus polinizadores?

Eu ja falei uma vez sobre como fazer gráficos de rede aqui, mas podemos explorar também várias propriedades de redes interessantes para a biologia.

Primeiro vamos simular aqui quatro redes para ajudar a pensar, como eu não sei muito bem como começar, vamos usar o modelo de Erdos para a simulação da rede de interações. Alias Paul Erdos era um cara muito interessante, diria até fascinante quanto ao seus estilo de vida, vale a pena uma olhadinha na história dele.

Mas por enquanto vamos olhar as redes que criamos onde bolinhas verdades são espécies de plantas e quadradinhos azuis são espécies de polinizadores.

01

Bem, eu fiz 4 redes, todas tem 15 polinizadores (linhas) e 15 espécies de plantas (colunas). As plantas dependem dos polinizadores para reprodução, e os polinizadores também dependem das flores para viver (de forma geral, existem excessões). Isso nos leva a conclusão óbvia que todo mundo que a gente ta vendo ai tem pelo menos 1 interação, cada polinizador visita pelo menos um espécie de flor.

Mas além disso podemos pensar em alguma possibilidades pensando nos polinizadores. Uma possibilidade é o polinizador ser generalista, visitar várias flores. Isso é bom já que ele deve encontrar flores mais facilmente, porém ele tem que estar preparado para ser um “pau para toda obra” já que cada flor tem uma forma (morfologia), as vezes disponibiliza seus recursos de uma forma diferente umas das outras, o que pode tornar esse polinizador meia sola em muitas espécies de plantas.

Outra estratégia que pode ser interessante se pensarmos bem é ser especializado. Ja que melhor que pegar só um pouquinho de comida de porta em porta que se bate, é ganhar um montão de comida em uma porta. Mas diferente do generalista, o especialista só consegue o que quer de uma única espécie de planta.

Em termo de energia ou tanto de comida que você consegue coletar, as duas estratégias podem ser equivalente, já que passar bastante tempo voando para pegar bastante pólen em uma única flor ou voar menos e ir pegando um pouquinho de pólen em cada flor pode dar na mesma no final das contas, mas existe outra coisa a se levar em conta, um polinizador generalista pode chegar em novos ambientes facilmente, já que ele consegue comida em varias espécies, enquanto o especialista só vai chegar no ambiente depois que sua “fonte de comida” chegou.

Até agora só consideramos comida mas podem ter flores outras utilidades para os polinizadores e existem outros tipos de interações. Existem ainda os interesses das plantas em todo esse jogo relações.

Mas talvez antes de ficar pensando em possibilidades e mais possibilidades, casos e mais casos, vamos olhar como é essa rede de interações como um todo.

Vamos olhar primeiro o exemplo 1 da perspectiva dos polinizadores.
Nele a gente vê bastante espécies de polinizadores especialista, que só visitam uma espécie de planta, como 2 casos de exclusividade polinizador e planta no meio do gráfico, 3 casos de grupinho envolvendo 2 polinizadores e 2 plantas e um grupão maior no canto inferior direito mais ou menos.

Mas a gente pode contar quantas plantas cada polinizador visita e fazer o seguinte gráfico;

02

Então vamos la, no eixo x temos o número de interações e no eixo y temos a contagem de polinizadores. Ou seja exitem 15 polinizadores, e como todo mundo interage com pelo menos uma planta (o obvio, senão nem na matriz de interação eles estariam), nos temos que com 1 planta os 15 polinizadores interagem, desses 15 polinizadores, apenas 7 interagem com 2 espécies de plantas e um apenas interagiu com 3 plantas.

Certo mas esse gráfico vai ser dependente do número de polinizadores no sistema, mas podemos resolver esse problema transformando tudo em porcentagem, ou seja dividir pelo número de polinizadores total. Então ficamos com…

03

Exatamente o mesmo gráfico, mas notem que agora o eixo X vai de 0 a 1, ou seja é independente do número de polinizadores, isso facilita se quisermos comparar características que vemos nesse gráfico com outras comunidades, independente do número de espécies de polinizadores. Mas ao invéz de barrinhas de contagem, podemos simplesmente usar bolinhas e ficamos com isso.

04

Agora sim, bolinha é bom para ajustar modelos, mas na verdade ainda estamos olhando a mesma coisa do primeiro gráfico de barras, apenas alguns ajustes aqui e ali.

Nesse caso um dos modelos que ajuda a descrever essa comunidade segundo o Jordano (2003) é o exponencial, é um modelo com dois parâmetros que ele chama de b e gamma. Gamma é a letra grega  gamma

E o modelo é assim:

 p(k) = b cdot exp(-gamma cdot k)

Mas notem que no artigo o Jordano (2003) não cita o b, que não fala muita coisa, mas ele não usa o símbolo de igual, ele usa o símbolo de aproximadamente, acho que para não citar o b, não entendo muito o porque disso, mas vamos em frente, é vivendo e aprendendo.

Bem podemos ajustar esse modelo não linear no R usando a função nls(). E ajustando temos…

05

Agora que porra é essa? Bem esse modelo (a linha vermelha) ta nos dizendo como muda a distribuição acumulada de interações, ou seja, numa comunidade onde todo mundo é especialista como essa do exemplo 1, as barrinhas dos gráficos la em cima descem muito rápido, formam uma escada muito íngreme, de degraus altos, enquanto tivessem mais generalistas, a escadinha seria bem menos íngreme e subir essa escadinha seria fácil. E é isso que o gamma do modelo ta falando.
Podemos fazer o mesmo procedimento para os 4 exemplos de comunidade e ver o gamma neles.

coef(m1.exe1) gamma b 0.9372338 2.5877177 coef(m1.exe2) gamma b 0.5020859 1.6452600 coef(m1.exe3) gamma b 0.2602011 1.3744055 coef(m1.exe4) gamma b 0.1417438 1.2410299

E sem olhar para o b por enquanto, mas somente o gamma que nos interessa agora, sabemos que a comunidade 1 é a que tem mais especialistas e menos generalistas, enquanto a comunidade do exemplo 4 tinha muitos mais generalistas e poucos especialistas, olhem o emaranhado que são as interações no exemplo 4.

Isso se reflete nos valores de gamma, olhem para o exemplo 4, ele é bem próximo de zero, bem baixinho 0.14, mas o exemplo 1, cheio de especialistas e quase nenhum generalista forte temos um valor de gamma de 0.93, quase 1.

Ou seja valor alto de gamma é uma escada íngreme naquele gráfico de barras, sinal de especialistas na comunidade, muitos especialistas.
Valor baixo de gamma é sinal de mais generalistas.
Agora temos um pouco de intuição sobre o que o modelo esta ajustando não?

Mas o Jordano citou três modelos no artigo dele, ajustar o modelo desse jeito vai falhar em um caso que imagino. No caso de todo mundo ser generalista. Isso eu imagino pois a gente pode tentar olhar o que acontece com as curvas de acordo com a mudança de gamma.
Como olhando a formula não me traz muitos insights pois sou péssimo de matemática, eu resolvo ela para tudo que é valor de gamma e fico olhando o que acontece.

06

E olha ai, vemos como quanto maior o valor de gamma, mais íngreme é a subida, e biologicamente isso é mais gente com 1 interação apenas, ou seja, especialistas dominam na comunidade. Mas o gamma vai ser no máximo 1, quando ele for 1 o modelo vai ser uma reta inclinada, e no caso de somente haver generalistas, precisaríamos de uma reta não inclinada, ou seja esse modelo ia falhar. Mas o Jordano foi legal deixou outras opções. Mas investigaremos elas depois 🙂

Tudo muito legal mas foram um milhão de comandos, eu não sei fazer tudo isso denovo. Tire a preocupação do seu coração que o senhor Carsten foi muito gentil em resumir todo esse processo numa função do pacote bipartite. La é só digitar degreedistr(sua comunidade) que toda essa magia se realiza para os polinizadores e para as plantas de uma vez so, apesar que não abordamos plantas aqui mas a idéia continua a mesma.

O Resultado vai sair assim e com um gráfico que podemos comparar com o que fizemos até aqui, mas lembrando que aqui estamos olhando para o exemplo 4, a função do bipartite não vai te deixar ajustar o modelo se vc tiver somente 3 barrinhas como nos temos no exemplo 1, o que é fato ja que são poucos graus de liberdade, mas que dificultaria nossa didática aqui.

07

degreedistr(exe4, plot.it=T) $`lower trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1417439 0.01486731 1.210787e-05 0.9717753 -18.632825 power law 0.4722450 0.09654721 1.206915e-03 0.8745991 -4.366322 truncated power law -0.3072618 0.12245915 4.045036e-02 0.9856348 -23.467338 $`higher trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1377984 0.01249272 1.573018e-06 0.9746622 -23.366621 power law 0.4814095 0.08664691 3.536921e-04 0.8791545 -7.092832 truncated power law -0.2837947 0.11314399 3.646958e-02 0.9863076 -28.219390

E lembrando que polinizadores que eram nossas linhas são o “lower trophic level” nesse caso, e olhem que o modelo exponential tem um estimate de 0.14…

coef(m1.exe4) gamma b 0.1417438 1.2410299

E nosso estimate é o mesmo valor, 0.14, ohhhhh, será que porque eu copiei o procedimento dele? 🙂

Bem é isso, fica no script todo o procedimento que fiz para gerar as figuras e a estimativa, depois eu preparo mais posts com os outros modelos, principalmente que essa propriedade tem outro grande interesse em biologia, se um gamma maior significa mais especialistas, um gamma maior também traz algo a respeito da susceptibilidade da comunidade a extinção em cascata.

#Semente para os resultados serem reproduziveis
set.seed(123)
###################################################
#Função para gerar rede de interação baseado no modelo de Erdős–Rényi
#Creditos a Edwin Grappin, vi o post no R-bloggers
#http://probaperception.blogspot.com.br/2012/12/have-you-tried-to-understand-your.html
###################################################
generateER = function(n = 100, p = 0.5){
map = diag(rep(1, n))
link = rbinom(n*(n-1)/2, 1,p)
t = 1
for(j in 2:n){
for(i in 1:(j-1)){
map[i,j] = link[t]
t = t + 1
}
}
dimnames(map)<-list(paste("Pol",1:n),paste("Sp",1:n))
return(map)
}
###########################
#Minha função para transformar matriz em lista para o igraph
##############################
mat2lis<-function(matriz) {
lista<-dimnames(matriz)
dados<-data.frame(Pol=rep(lista[[1]],length(lista[[1]])),
SP=rep(lista[[2]],each=length(lista[[2]])),
int=c(matriz))
return(dados[dados$int==1,1:2])
}
############################
library(igraph)
 
#Gerando exemplos
exe1<-generateER(n=15,p=0.1)
exe2<-generateER(n=15,p=0.2)
exe3<-generateER(n=15,p=0.4)
exe4<-generateER(n=15,p=0.8)
 
#Transformando na forma de dados que o igraph usa
exe1.l<-mat2lis(exe1)
exe1.g <-graph.data.frame(exe1.l, directed=F)
 
exe2.l<-mat2lis(exe2)
exe2.g <-graph.data.frame(exe2.l, directed=F)
 
exe3.l<-mat2lis(exe3)
exe3.g <-graph.data.frame(exe2.l, directed=F)
 
exe4.l<-mat2lis(exe4)
exe4.g <-graph.data.frame(exe4.l, directed=F)
 
#Plots das redes
par(mfrow=c(2,2))
 
plot(exe1.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 1")
 
plot(exe2.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 2")
 
plot(exe3.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 3")
 
plot(exe4.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.5,
vertex.label.color="black"
)
title("Exemplo 4")
 
#Ajustando o modelo exponencial para polinizadores
#Jordano, P., Bascompte, J. and Olesen, J. M. (2003)
#Invariant properties in coevolutionary networks of
#plant-animal interactions. Ecology Letters 6, 69–81
 
#Rede do Exemplo 1
#Número de interações
k<-sum(exe1)
#Número de interações por espécie de polinizador
pol.n<-rowSums(exe1)
#Número de polinizadores por quantidade de interação
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
#Juntando os dados
dados.pol1 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
#Ajustando modelo exponencial aos dados
m1.exe1<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol1)
 
#Grafico 2
barplot(P.pol,names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 3
barplot(dados.pol1[,2],names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 4
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
 
#Grafico 5
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
curve(coef(m1.exe1)[2]*exp(-coef(m1.exe1)[1]*x),1,5,add=T,lwd=2,lty=2,
cex=1.5,col="red")
 
#Exemplo 2
k<-sum(exe2)
S<-sum(dim(exe2))
pol.n<-rowSums(exe2)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol2 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe2<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol2)
 
#Exemplo 3
k<-sum(exe3)
pol.n<-rowSums(exe3)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol3 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe3<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol3)
 
#Exemplo 4
k<-sum(exe4)
pol.n<-rowSums(exe4)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol4 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe4<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol4)
 
#Extraindo coeficientes
coef(m1.exe1)
coef(m1.exe2)
coef(m1.exe3)
coef(m1.exe4)
 
#Niveis de Gamma
graus<-seq(0.1,0.9,0.1)
 
#GRafico 6
jpeg("06.jpg")
plot(0,0,type="n",xlim=c(0,10),ylim=c(0,1),frame.plot=F,
xlab="Número de interações",
ylab="Distribuição Cumulativa do número de interações")
for(i in 1:length(graus)) {
curve(1*exp(-graus[i]*x),0,10,add=T,col=i,lty=2,lwd=2)
}
legend("topright",lty=2,lwd=2,col=c(1:length(graus)),legend=paste("Gamma é ",graus))
dev.off()
 
#Grafico e ajuste usando o bipartite para o exemplo 4
library(bipartite)
degreedistr(exe4, plot.it=T)

Gráficos de Rede no R com o Pacote Igraph

Quando fui olhar pela primeira vez como fazer gráficos de rede no R, a primeira parada foi o pacote bipartite.
O pacote bipartite além de calcular as métricas e estatísticas mais comuns para entender como redes funcionam, tem duas opções para gráficos, as funções plotweb e visweb.

Então quando temos nossa matriz de interações de forma pouco visual:

> table(dados) flor pol Flor 1 Flor 2 Flor 3 Flor 4 Flor 5 Polinizador 1 0 0 1 1 0 Polinizador 2 0 1 0 1 0 Polinizador 3 1 0 0 0 0 Polinizador 4 3 0 0 0 2 Polinizador 5 0 0 2 1 2

Onde o valor é o número de vezes que observamos determinada interação. Assim fica difícil responder facilmente qual polinizador é mais especialista, qual é mais generalista e assim vai. Mas a função plotweb() resolveu bem esse problema.

Mas, apesar de essa figura ficar muito bonita, eu senti dificuldade de conseguir, por exemplo, grudar uma filogenia na rede, já que a posição dos polinizadores não é simétrica, e não é possível controlar de forma simples essa característica. Mas o pacote bipartite tem como dependência o pacote igraph, e fui olhar do que se tratava.

O pacote igraph é exatamente especifico para plotar redes, além de cálculo de métricas e “coisas” de rede. Ele é na verdade uma biblioteca de funções e esta implementado para mais de uma linguagem de computador, como por exemplo para python também.
Mas ele trabalha com uma classe própria de dados, que não é difícil de entender. A matriz acima nele seria representada mais ou menos assim:

> unique(dados) pol flor 1 Polinizador 4 Flor 5 2 Polinizador 1 Flor 3 3 Polinizador 5 Flor 5 4 Polinizador 4 Flor 1 5 Polinizador 2 Flor 4 7 Polinizador 5 Flor 4 8 Polinizador 2 Flor 2 10 Polinizador 5 Flor 3 11 Polinizador 1 Flor 4 13 Polinizador 3 Flor 1

Ou seja, quem interage com quem. Então convertemos esses dados para a classe graph usando a função graph.data.frame() nos dados organizados dessa forma e estamos prontos para fazer os gráficos.

Podemos dar plot nos dados e vemos o que acontece…

Ele faz um gráfico de rede, espalhando todas as plantas e polinizadores e ligando conforme as interações. Na verdade ele tem varias formas de espalhar os nodes, que no nosso caso são polinizadores e flores, e isso é muito interessante pois proporciona varias possibilidades de exaltar varias características da rede, como quem são os especialistas, como são as divisões de flores e polinizadores, além disso é possível modificar formas, cores e muitas outras opções.

Bem mas são funções que fazem esse posicionamento, então podemos ao invés de fornecer uma função, podemos simplesmente dar o posicionamento que desejamos, por exemplo filas emparelhadas como na função plotweb(), claro que a função plotweb ele tem algumas opções de reorganizar a ordem para diminuir a sobreposição entre os riscos que representam as interações. Mas de forma simples podemos fazer no igraph algo similar ao plotweb().

A vantagem que eu achei em usar o igraph é em controlar o tamanho dos nodes, dos símbolos que representam polinizadores e flores, e assim fica mais simples a tarefa de por exemplo colar a filogenia no gráfico de rede, tarefa que não consegui fazer usando o pacote bipartite.

library(bipartite)
library(igraph)
 
#gerando dados de exemplo
set.seed(15)
pol.s<-paste("Polinizador",1:5)
flor.s<-paste("Flor",1:5)
dados<-data.frame(pol=sample(pol.s,15,replace=T),flor=sample(flor.s,15,replace=T))
 
#matriz de interações
table(dados)
 
#Rede com bipartite
plotweb(table(dados))
 
#Convertendo dados
dados.g <-graph.data.frame(unique(dados), directed=F)
 
#Figura 2
plot(dados.g,vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)),vertex.size=25)
 
#Figura 3
par(mfrow=c(2,2))
plot(dados.g, layout=layout.circle,main="Círculo")
plot(dados.g, layout=layout.sphere,vertex.shape="square",main="Esfera")
plot(dados.g, layout=layout.fruchterman.reingold,vertex.color="red",main="Fruchterman Reingold")
plot(dados.g, layout=layout.fruchterman.reingold.grid,
main="Fruchterman reingold Grid",edge.color="red")
 
#Criando posicionamento dos nodes
posicao<-matrix(c(rep(1:5,2),rep(1:2,each=5)),ncol=2,nrow=10)
 
#Rede similar ao bipartite
plot(dados.g, layout=posicao,vertex.color=rep(c("green","yellow"),each=5),
vertex.shape=rep(c("square","circle"),each=5),vertex.size=23,
vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)))
 
#Ligando rede com filogenia
library(ape)
 
par(mfrow=c(2,1),mar=c(1,9,1,9))
plot(rtree(5),show.tip.label=F,direction="downwards",use.edge.length = F,edge.width = 2)
title("Filogenia")
par(mar=c(0,0,0,0))
plot(dados.g, layout=posicao+1,vertex.color=rep(c("green","yellow"),each=5),
vertex.shape=rep(c("square","circle"),each=5),vertex.size=35,
vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)),edge.color="black",
edge.width=3)