Granovetter’s impossible triad
Neste post abordaremos a última palavra em modelos estatísticos para dados relacionais os exponential random graph models, ou modelos p* (pê-estrela). Iniciaremos com um post mais técnico focado nos comandos e possíveis parâmetros do pacote statnet. No próximo post, abordaremos questões mais conceituais relacionadas aos modelos bem como sua interpretação.
Além disso, este post está organizado neste belo e eficiente layout chamado Tufte Handout inspirado nos livros e comunicações do prof. Edward Tufte. Para mais informações sobre o uso, vantagens e questões técnicas, consulte este site.
Os modelos p* foram desenvolvidos para dar conta de dados que assumem interdependência das observações. Ora, os modelos econométricos clássicos (regressão linear, regressão logística, etc.) assumem independência das observações e portanto não devem ser aplicados a dados relacionais. Em outras palavras, se A conhece C, a existência de um laço entre B e C altera a probabilidade de existência de um laço entre B e A1 Granovetter, M. S. (1973). The strength of weak ties. American journal of sociology, 1360-1380..
Granovetter’s impossible triad
O modelo p* pode ser definido por
\[Pr(Y=y) \quad = \quad \left(\frac{1}{k}\right) exp \left\{ \sum_{A} \eta_A g_A (\textbf{y}) \right\}\]
onde Y é o grafo teórico simulado, y é o grafo observado, \(\sum_{A}\) é a soma de todas as configurações A, \(\eta_A\) é o parâmetro estimado correpondente à configuração A, \(g_A (\textbf{y})\) é a estatística da configuração A observada no grafo y e k é uma constante que assegura a distribuição de probabilidades2 Robins, Garry et al. An introduction to exponential random graph (p*) models for social networks. Social networks, Elsevier, v. 29, n. 2, p. 173–191, 2007..
Para estimar um modelo p* no R, usaremos os comandos da suite statnet e o pacote igraph para auxílio com a importação dos dados quando for o caso. Vamos ao primeiro exemplo.
Vamos tentar modelar a rede casamentos de famílias de Florença do famoso trabalho de Padgett (1994)3 Padgett, John F. (1994) Marriage and Elite Structure in Renaissance Florence, 1282-1500. Paper delivered to the Social Science History Association.. Primeiro, investiguemos a rede.
library(statnet)
library(sand)
library(intergraph)data(florentine)
flomarriage## Network attributes:
## vertices = 16
## directed = FALSE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## total edges= 20
## missing edges= 0
## non-missing edges= 20
##
## Vertex attribute names:
## priorates totalties vertex.names wealth
##
## No edge attributes
plot(flomarriage, displaylabels=T,
label.cex=0.8, main="Padgett Florentine Families")Investiguemos agora os atributos. Como acho mais fácil fazer isso usando o pacote igraph, vamos a ele.
marriage <- asIgraph(flomarriage)
V(marriage)$priorates## [1] 53 65 0 12 22 0 21 0 53 0 42 0 38 35 74 0V(marriage)$totalties## [1] 2 3 14 9 18 9 14 14 54 7 32 1 4 5 29 7V(marriage)$wealth## [1] 10 36 55 44 20 32 8 42 103 48 49 3 27 10 146 48Podemos plotar a rede representando algum atributo. Vamos representar a riqueza com o tamanho dos vértices. Como o pacote igraph lê os valores absolutos e o grafo fica muito poluído, vamos dividir o vetor de valores de riqueza por 4 para uma visualização mais limpa.
plot(marriage, vertex.label.cex=0.9,
vertex.label=V(marriage)$vertex.names,
vertex.size=V(marriage)$wealth/4,
vertex.color="lightblue",
vertex.label.color="black",
main = "Padgett Florentine Families",
xlab="Size = Wealth")
Agora vamos à modelagem. Primeiro, é comum estimarmos o modelo nulo4 Este modelo é também conhecido como \(P_1\) e foi desenvolvido por Holland e Leinhardt., ou seja, um modelo apenas com o parâmetro edges que corresponde ao intercepto dos modelos econométricos comuns. Este modelo toma a díade por unidade de análise e estima a “propensão dos atores a escolherem outros atores, serem escolhidos por outros, a fazer escolhas recíprocas e a tendência média (parâmetro de densidade) a interagir com os outros (Lazega & Higgins, 2014, p. 81)5 Lazega, E., & Higgins, S. S. (2014). Redes sociais e estruturas relacionais. Belo Horizonte, MG: Fino Traço..
fit1 <- ergm(flomarriage~edges)summary(fit1)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: flomarriage ~ edges
##
## Iterations: 5 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -1.6094 0.2449 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 166.4 on 120 degrees of freedom
## Residual Deviance: 108.1 on 119 degrees of freedom
##
## AIC: 110.1 BIC: 112.9 (Smaller is better.)O parâmetro negativo e significativo nos mostra que esta rede se desenvolve com menos laços entre seus atores do que esperaríamos numa rede aleatória6 No próximo post pretendo explicar melhor o conceito teórico de rede aleatória e como ele é mobilizado na interpretação dos ERGMs.. Em seguida, vamos estimar um modelo controlando por triângulos e por distribuição de laços compartilhados geometricamente pesada7 Geometrically weighted edgewise shared partner distribution. (gwesp). Vamos usar os comandos formula e summary.statistics para contar as estatísticas existentes na rede original.
model2 <- formula(flomarriage~edges+triangle+
gwesp(1,fixed=T))
summary.statistics(model2)## edges triangle gwesp.fixed.1
## 20.000000 3.000000 8.632121Em seguida, “rodamos” o modelo.
fit2 <- ergm(model2)summary(fit2)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: flomarriage ~ edges + triangle + gwesp(1, fixed = T)
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -1.7307 0.3754 0 <1e-04 ***
## triangle -4.3617 8.7671 0 0.620
## gwesp.fixed.1 1.7068 3.1747 0 0.592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 166.4 on 120 degrees of freedom
## Residual Deviance: 107.6 on 117 degrees of freedom
##
## AIC: 113.6 BIC: 122 (Smaller is better.)Às vezes é útil avaliar o processo de amostragem por Monte Carlo Markov Chain (MCMC). Para isso, usamos o comando plot().
plot(fit2)## Sample statistics summary:
##
## Iterations = 16384:4209664
## Thinning interval = 1024
## Number of chains = 1
## Sample size per chain = 4096
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## edges -0.03442 3.995 0.06242 0.06242
## triangle -0.02246 1.853 0.02896 0.02896
## gwesp.fixed.1 -0.05941 5.161 0.08064 0.08064
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## edges -8.000 -3 0.000e+00 3 8.00
## triangle -3.000 -1 0.000e+00 1 4.00
## gwesp.fixed.1 -8.632 -3 2.119e-12 3 11.26
##
##
## Sample statistics cross-correlations:
## edges triangle gwesp.fixed.1
## edges 1.0000000 0.7375073 0.742671
## triangle 0.7375073 1.0000000 0.998088
## gwesp.fixed.1 0.7426710 0.9980880 1.000000
##
## Sample statistics auto-correlation:
## Chain 1
## edges triangle gwesp.fixed.1
## Lag 0 1.000000000 1.000000000 1.000000000
## Lag 1024 0.001902441 0.006472082 0.005209743
## Lag 2048 0.006939257 0.011665822 0.011808763
## Lag 3072 0.016012409 0.010663944 0.012128877
## Lag 4096 0.030240606 0.005898352 0.005916047
## Lag 5120 0.031856924 0.005182487 0.004993409
##
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## edges triangle gwesp.fixed.1
## -0.08226 0.91539 0.91474
##
## Individual P-values (lower = worse):
## edges triangle gwesp.fixed.1
## 0.9344393 0.3599844 0.3603262
## Joint P-value (lower = worse): 0.6102331 .## Loading required namespace: latticeExtra##
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).Para verificar o ajuste e a convergência do modelo, usamos o comando gof().
gof2 <- gof(fit2)
par(mfrow=c(1,3))
plot(gof2)O modelo parece se ajustar bem aos dados. Este modelo convergiu e podemos confiar nos parâmetros estimados para explicar a emergência de nossa rede observada. Neste caso, observando os resultados, apesar de triângulos e distribuição de laços compartilhados terem estimadores respectivamente negativo e postivo, eles não são estatisticamente significativos indicando que essas configurações não são importantes para entendermos o surgimento da rede.
Agora, como outro exemplo, vamos trabalhar com os dados dos advogados de Lazega.
Estes dados estão no pacote sand. Precisamos transformá-los num objeto de classe network para realizar a estimação do p*.
lazega <- upgrade_graph(lazega)
lawyers <- asNetwork(lazega)fit1 <- ergm(lawyers~edges)summary(fit1)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: lawyers ~ edges
##
## Iterations: 5 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -1.4992 0.1031 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 873.4 on 630 degrees of freedom
## Residual Deviance: 598.8 on 629 degrees of freedom
##
## AIC: 600.8 BIC: 605.2 (Smaller is better.)Depois de estimar o modelo \(P_1\), vamos estimar um modelo com triângulos e distribuição de laços compartilhados.
model2 <- formula(lawyers~edges+triangle+
gwesp(1,fixed=T))
summary.statistics(model2)## edges triangle gwesp.fixed.1
## 115.0000 120.0000 213.1753fit2 <- ergm(model2)summary(fit2)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: lawyers ~ edges + triangle + gwesp(1, fixed = T)
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -3.6720 0.3211 0 < 1e-04 ***
## triangle -0.2216 0.2190 0 0.31215
## gwesp.fixed.1 0.9822 0.2538 0 0.00012 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 873.4 on 630 degrees of freedom
## Residual Deviance: 519.8 on 627 degrees of freedom
##
## AIC: 525.8 BIC: 539.2 (Smaller is better.)Vamos ver o ajuste desse modelo aos dados e sua convergência.
gof2 <- gof(fit2)
par(mfrow=c(1,3))
plot(gof2)Vamos tentar elaborar um modelo um pouco mais complexo usando algumas estatísticas a mais além dos atributos dos atores. Dessa vez vamos modelar ainda a distribuição de grau e os atributos gênero, escritório que pertence, escola de origem, idade e senioridade. Atributos qualitativos são inseridos no modelo por semelhança. Atributos quantitativos são inseridos, neste caso, somando-se os valores de díades para todos os atores na rede. A lista de estatísticas possíveis é longa e podemos investigá-la com o comando help("ergm-terms"). Vamos ao modelo:
model3 <- formula(lawyers~edges+triangle+
gwesp(1,fixed=T)+gwdegree(1,fixed=T)+
nodematch("Gender")+nodematch("Office")+
nodematch("School")+nodecov("Age")+
nodecov("Seniority"))
summary.statistics(model3)## edges triangle gwesp.fixed.1 gwdegree
## 115.0000 120.0000 213.1753 79.2000
## nodematch.Gender nodematch.Office nodematch.School nodecov.Age
## 99.0000 85.0000 36.0000 10526.0000
## nodecov.Seniority
## 4687.0000fit3 <- ergm(model3)summary(fit3)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: lawyers ~ edges + triangle + gwesp(1, fixed = T) + gwdegree(1,
## fixed = T) + nodematch("Gender") + nodematch("Office") +
## nodematch("School") + nodecov("Age") + nodecov("Seniority")
##
## Iterations: 3 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges 1.238827 2.407460 0 0.607031
## triangle -0.148988 0.224826 0 0.507781
## gwesp.fixed.1 0.943846 0.280924 0 0.000828 ***
## gwdegree 0.468504 0.633291 0 0.459706
## nodematch.Gender 0.342240 0.218447 0 0.117696
## nodematch.Office 0.985246 0.182096 0 < 1e-04 ***
## nodematch.School 0.004201 0.256795 0 0.986953
## nodecov.Age -0.052819 0.018697 0 0.004879 **
## nodecov.Seniority -0.025280 0.015257 0 0.098026 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 873.4 on 630 degrees of freedom
## Residual Deviance: 475.5 on 621 degrees of freedom
##
## AIC: 493.5 BIC: 533.6 (Smaller is better.)Para exemplos de interpretação, veja Lusher, Koskinen & Robins (2012) e Salej Higgins et al (2014)8 Salej Higgins, S., Ribeiro, A., Botrel de Vasconcellos, M., & Estevão Barbosa, J. (2014). L’émergence d’une structure cœur-périphérie dans un réseau brésilien de copublications en sciences comptables. Communiquer. Revue de communication sociale et publique, (12), 43-60.
Vamos agora investigar o ajuste do modelo aos dados e sua convergência.
gof3 <- gof(fit3)
par(mfrow=c(1,3))
plot(gof3)
Quando trabalhamos com modelos p*, as possibilidades são enormes. Devemos testar várias configurações de estatísticas modeladas de acordo com o que pretendemos investigar na rede em questão e com o ajuste do modelo aos dados.
O algoritmo default pelo qual o pacote statnet realiza a estimação dos modelos p* é o Geyer-Thompson (1992) (MCMLE)9 Lusher, D., Koskinen, J., & Robins, G. (2012). Exponential random graph models for social networks: Theory, methods, and applications. Cambridge University Press.. Este algoritmo costuma retornar os melhores resultados mas, por ser extremamente exigente (os modelos devem convergir 2 vezes, do contrário, o comando dá erro e para) pode não rodar para algumas redes mais extensas. O pacote statnet permite o uso de mais três algoritmos além do MCMLE: o algoritmo Robbins-Monro (usado pela família de softwares PNET), Metropolis-Hastings (aproximação estocástica) e Stepping. Vamos realizar a estimação do modelo 3 para os dados de Lazega com os outros 3 algoritmos e vamos compará-los.
fit3.rb <- ergm(model3,
control=control.ergm(main.method = "Robbins-Monro"))
fit3.sa <- ergm(model3,
control=control.ergm(main.method = "Stochastic-Approximation"))
fit3.step <- ergm(model3,
control=control.ergm(main.method = "Stepping"))library(texreg)htmlreg(list(fit3, fit3.rb, fit3.sa, fit3.step),
single.row=T, caption.above = T,
custom.model.names=c("Geyer-Thompson", "Robbins-Monro","Metropolis-Hastings","Stepping"))| Geyer-Thompson | Robbins-Monro | Metropolis-Hastings | Stepping | ||
|---|---|---|---|---|---|
| edges | 1.24 (2.41) | -0.51 | 0.17 (2.59) | 1.32 (2.51) | |
| triangle | -0.15 (0.22) | -0.07 | -0.05 (0.22) | -0.15 (0.25) | |
| gwesp.fixed.1 | 0.94 (0.28)*** | 0.97 | 0.82 (0.32)* | 0.94 (0.29)** | |
| gwdegree | 0.47 (0.63) | 1.28 | 1.01 (0.69) | 0.46 (0.65) | |
| nodematch.Gender | 0.34 (0.22) | 0.83 | 0.25 (0.24) | 0.33 (0.23) | |
| nodematch.Office | 0.99 (0.18)*** | 0.49 | 0.73 (0.22)*** | 0.99 (0.18)*** | |
| nodematch.School | 0.00 (0.26) | -0.19 | 0.10 (0.24) | -0.00 (0.24) | |
| nodecov.Age | -0.05 (0.02)** | -0.04 | -0.04 (0.02)* | -0.05 (0.02)** | |
| nodecov.Seniority | -0.03 (0.02) | -0.02 | -0.02 (0.02) | -0.03 (0.02) | |
| AIC | 493.54 | 511.64 | 507.21 | 494.06 | |
| BIC | 533.55 | 551.65 | 547.22 | 534.08 | |
| Log Likelihood | -237.77 | -246.82 | -244.61 | -238.03 | |
| ***p < 0.001, **p < 0.01, *p < 0.05 | |||||
Vamos verificar o ajuste dos modelos estimados com os outros algoritmos.
gof3.rb <- gof(fit3.rb)
gof3.sa <- gof(fit3.sa)
gof3.step <- gof(fit3.step)
par(mfrow=c(1,3))
plot(gof3.rb,
main = "Goodness-of-fit 'Robbins-Monro'")
plot(gof3.sa,
main = "Goodness-of-fit 'Metropolis-Hastings'")
plot(gof3.step,
main="Goodness-of-fit 'Stepping'")
A estimação de um modelo p* pode ser um processo longo e cansativo. Não há testes ou métodos para escolha prévia de um modelo a não ser o bom e velho “tentativa e erro”. A partir do seu problema de pesquisa, escolha os parâmetros que deseja estimar e vá melhorando o modelo escolhendo outros parâmetros. Muita atenção: o modelo deve convergir e os parâmetros estimados devem estar dentro dos limites dos plots de Goodness-of-fit. O algoritmo Geyer-Thompson tem a vantagem (ou não) de não rodar quando o modelo não converge. Os outros algoritmos rodam de qualquer modo e checamos se o modelo está degenerado pelo GOF. É importantíssimo salientar que o modelo deve convergir completamente. Um modelo que não converge em algum de seus parâmetros não quer dizer absolutamente nada (Lusher, Koskinem & Robins, 2012)10 Lusher, D., Koskinen, J., & Robins, G. (2012). Exponential random graph models for social networks: Theory, methods, and applications. Cambridge University Press..