HW04
JUNG KYU HYEONG
2020 6 9
1. Implement all commands for ERGMs.
library(sand)## Loading required package: igraph##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## union## Loading required package: igraphdata##
## Statistical Analysis of Network Data with R
## Type in C2 (+ENTER) to start with Chapter 2.data(lazega)
lazega=upgrade_graph(lazega)
A=get.adjacency(lazega)
V.attrs=get.data.frame(lazega,what="vertices")
library(ergm)## Loading required package: network## network: Classes for Relational Data
## Version 1.16.0 created on 2019-11-30.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Martina Morris, University of Washington
## Skye Bender-deMoll, University of Washington
## For citation information, type citation("network").
## Type help("network-package") to get started.##
## Attaching package: 'network'## The following objects are masked from 'package:igraph':
##
## %c%, %s%, add.edges, add.vertices, delete.edges, delete.vertices,
## get.edge.attribute, get.edges, get.vertex.attribute, is.bipartite,
## is.directed, list.edge.attributes, list.vertex.attributes,
## set.edge.attribute, set.vertex.attribute##
## ergm: version 3.10.4, created on 2019-06-10
## Copyright (c) 2019, Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Carter T. Butts, University of California -- Irvine
## Steven M. Goodreau, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Martina Morris, University of Washington
## with contributions from
## Li Wang
## Kirk Li, University of Washington
## Skye Bender-deMoll, University of Washington
## Chad Klumb
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").## NOTE: Versions before 3.6.1 had a bug in the implementation of the bd()
## constriant which distorted the sampled distribution somewhat. In
## addition, Sampson's Monks datasets had mislabeled vertices. See the
## NEWS and the documentation for more details.## NOTE: Some common term arguments pertaining to vertex attribute and
## level selection have changed in 3.10.0. See terms help for more
## details. Use 'options(ergm.term=list(version="3.9.4"))' to use old
## behavior.lazega.s=network::as.network(as.matrix(A),directed=F)
network::set.vertex.attribute(lazega.s,"Office",V.attrs$Office)
network::set.vertex.attribute(lazega.s,"Practice",V.attrs$Practice)
network::set.vertex.attribute(lazega.s,"Gender",V.attrs$Gender)
network::set.vertex.attribute(lazega.s,"Seniority",V.attrs$Seniority)my.ergm.bern=formula(lazega.s~edges)
my.ergm.bern## lazega.s ~ edgessummary(my.ergm.bern)## edges
## 115my.ergm=formula(lazega.s~edges+kstar(2)+kstar(3)+triangle)
summary(my.ergm)## edges kstar2 kstar3 triangle
## 115 926 2681 120my.ergm=formula(lazega.s~edges+gwdegree(1,fixed=T)+gwesp(1,fixed=T))
summary(my.ergm)## edges gwdeg.fixed.1 gwesp.fixed.1
## 115.0000 79.2000 213.1753lazega.ergm=formula(lazega.s~edges+gwesp(log(3),fixed=T)+nodemain("Seniority")+nodemain("Practice")+match("Practice")+match("Gender")+match("Office"))
summary(lazega.ergm)## Warning: `set_attrs()` is deprecated as of rlang 0.3.0
## This warning is displayed once per session.## edges gwesp.fixed.1.09861228866811
## 115.0000 222.9108
## nodecov.Seniority nodecov.Practice
## 4687.0000 359.0000
## nodematch.Practice nodematch.Gender
## 72.0000 99.0000
## nodematch.Office
## 85.0000lazega.ergm.fit=ergm(lazega.ergm)## Starting maximum pseudolikelihood estimation (MPLE):## Evaluating the predictor and response matrix.## Maximizing the pseudolikelihood.## Finished MPLE.## Starting Monte Carlo maximum likelihood estimation (MCMLE):## Iteration 1 of at most 20:## Optimizing with step length 1.## The log-likelihood improved by 0.3309.## Step length converged once. Increasing MCMC sample size.## Iteration 2 of at most 20:## Optimizing with step length 1.## The log-likelihood improved by 0.02566.## Step length converged twice. Stopping.## Finished MCMLE.## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check
## for degeneracy, use the mcmc.diagnostics() function.anova(lazega.ergm.fit)## Analysis of Variance Table
##
## Model 1: lazega.s ~ edges + gwesp(log(3), fixed = T) + nodemain("Seniority") +
## nodemain("Practice") + match("Practice") + match("Gender") +
## match("Office")
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 630 873.37
## Model 1: 7 414.07 623 459.30 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lazega.ergm.fit)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: lazega.s ~ edges + gwesp(log(3), fixed = T) + nodemain("Seniority") +
## nodemain("Practice") + match("Practice") + match("Gender") +
## match("Office")
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -6.980114 0.696985 0 -10.015 < 1e-04 ***
## gwesp.fixed.1.09861228866811 0.591482 0.086461 0 6.841 < 1e-04 ***
## nodecov.Seniority 0.024443 0.006228 0 3.925 < 1e-04 ***
## nodecov.Practice 0.393004 0.106972 0 3.674 0.000239 ***
## nodematch.Practice 0.773543 0.198434 0 3.898 < 1e-04 ***
## nodematch.Gender 0.717003 0.250599 0 2.861 0.004221 **
## nodematch.Office 1.159229 0.192486 0 6.022 < 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: 459.3 on 623 degrees of freedom
##
## AIC: 473.3 BIC: 504.4 (Smaller is better.)gof.lazega.ergm=gof(lazega.ergm.fit)
par(mfrow=c(2,2))
plot(gof.lazega.ergm)
2. Fit Karate network using ERGMs with following statistics, edge count + triangle count. You need to explain your results in detail.
data(karate)
karate=upgrade_graph(karate)
A=get.adjacency(karate)
v.attrs=get.data.frame(karate,what="vertices")먼저 karate 네트워크를 인접행렬로 분리하였다. v.attrs에서는 각 vertice의 Faction,name,label,color의 값이 들어가있다.
library(statnet.common)##
## Attaching package: 'statnet.common'## The following object is masked from 'package:base':
##
## orderkarate.s=network::as.network(as.matrix(A),directed=F)
karate.s=network::as.network(as.matrix(A),directed=F)
network::set.vertex.attribute(karate.s, "Faction", v.attrs$Faction)
#network::set.vertex.attribute(karate.s, "name", v.attrs$name)
network::set.vertex.attribute(karate.s, "label", v.attrs$label)
network::set.vertex.attribute(karate.s, "color", v.attrs$color)ergm 패키지는 1번에서 로드하였으므로 따로 불러오지 않았다. karate.s는 인접행렬을 네트워크 객체로 바꾼 형태이다.(ergm 데이터 분석을 위해서) 그리고 각 vertex의 특징을 네트워크 객체에 더하였다.
my.ergm.bern=formula(karate.s~edges)
my.ergm.bern## karate.s ~ edgessummary(my.ergm.bern)## edges
## 78summary()를 통해 위처럼 edge의 개수를 얻을 수 있다.
fit0=ergm(my.ergm.bern)## Starting maximum pseudolikelihood estimation (MPLE):## Evaluating the predictor and response matrix.## Maximizing the pseudolikelihood.## Finished MPLE.## Stopping at the initial estimate.## Evaluating log-likelihood at the estimate.summary(fit0)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: karate.s ~ edges
##
## Iterations: 5 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -1.823 0.122 0 -14.94 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 777.7 on 561 degrees of freedom
## Residual Deviance: 452.4 on 560 degrees of freedom
##
## AIC: 454.4 BIC: 458.7 (Smaller is better.)simple model을 통한 추정결과이다. edges의 변수만 존재하고 있으며 위 모델은 mcmc대신 로지스틱 회귀를 사용하여 풀었다.edges는 네트워크의 밀도를 의미하고 있어 연결 확률을 나타낸다고 볼 수 있다.
plogis(-1.823)## [1] 0.1390743connection probablity는 약 14%정도이다.
my.ergm=formula(karate.s~edges+triangle)
summary(my.ergm)## edges triangle
## 78 45문제에 나와있는 식을 통해 fit하고자 하였다. edge와 triangle의 개수를 확인할 수 있다.
karate.ergm.fit=ergm(my.ergm)## Starting maximum pseudolikelihood estimation (MPLE):## Evaluating the predictor and response matrix.## Maximizing the pseudolikelihood.## Finished MPLE.## Starting Monte Carlo maximum likelihood estimation (MCMLE):## Iteration 1 of at most 20:## Optimizing with step length 0.955074906238642.## The log-likelihood improved by 6.759.## Iteration 2 of at most 20:## Error in ergm.MCMLE(init, nw, model, initialfit = (initialfit <- NULL), : Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.위와 같은 모델을 사용하였을때 error가 발생하는 이유는 무엇인가?
수업시간에 교수님이 말씀하신 것처럼 # of triangles 은 dyadic dependent statistic이기 때문이다.3개의 edge가 있다면 하나의 edge는 나머진 2개의 edge와 depend하다.또한, 네트워크의 구조를 설명하는 많은 네트워크 statistics들은 dyadic dependent statisitcs이다. dyadic dependent statistics에서 발생하는 문제를 model degeneracy라 한다.
“dyadic dependent statistics make the other parameter estimation problem that is called as a model degeneracy”
만약 모델에서 하나의 edge를 지우게 된다면 아래와 같이 값이 변하게 된다.
Caption for the picture.
이러한 이유로 n이 커질수록 변화는 매우 커지게 될 것이다. 동시에 모수를 추정을 돕는 네트워크 sample generation은 매우 어려워진다. 이러한 model degeneracy를 막기위해 아래 3번에서 쓰이는 Geometrically Weighted Statistic를 사용한다.
3. Fit Karate network using ERGMs with following statistics, edge count + GWD(= 0.2) + GWESP(=0.2), and generate GOF plot. You need to describe the performance of GOF plot.
karate.ergm=formula(karate.s~edges+gwdegree(0.2,fixed=T)+gwesp(0.2,fixed=T))
summary(karate.ergm)## edges gwdeg.fixed.0.2 gwesp.fixed.0.2
## 78.00000 40.81246 73.438553번의 식을 통해 fitting 하였다. edge의 개수, gwdegree, gwesp의 값을 확인할 수 있다.
gwdegree: Geometrically weighted degree statistic with fixed 0.2
gwesp: Geometrically weighted edge-wise shared partnership statistic with fixed 0.2
karate.ergm.fit=ergm(karate.ergm)## Starting maximum pseudolikelihood estimation (MPLE):## Evaluating the predictor and response matrix.## Maximizing the pseudolikelihood.## Finished MPLE.## Starting Monte Carlo maximum likelihood estimation (MCMLE):## Iteration 1 of at most 20:## Optimizing with step length 0.719720659893507.## The log-likelihood improved by 3.53.## Iteration 2 of at most 20:## Optimizing with step length 1.## The log-likelihood improved by 0.7206.## Step length converged once. Increasing MCMC sample size.## Iteration 3 of at most 20:## Optimizing with step length 1.## The log-likelihood improved by 0.0394.## Step length converged twice. Stopping.## Finished MCMLE.## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check
## for degeneracy, use the mcmc.diagnostics() function.summary(karate.ergm.fit)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: karate.s ~ edges + gwdegree(0.2, fixed = T) + gwesp(0.2, fixed = T)
##
## Iterations: 3 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -3.7259 0.3959 0 -9.410 <1e-04 ***
## gwdeg.fixed.0.2 3.6672 2.5128 0 1.459 0.144
## gwesp.fixed.0.2 1.3001 0.2676 0 4.859 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 777.7 on 561 degrees of freedom
## Residual Deviance: 416.0 on 558 degrees of freedom
##
## AIC: 422 BIC: 435 (Smaller is better.)위 추정치를 통해 모델에서 개별변수들의 기여도를 확인할 수 있다. 1번의 단순모델과 다르게 위의 모델에서는 mcmc를 이용하였다. 또한, 2번과 다르게 error가 발생하지 않았다. 대략적인 MLE 계수를 찾기 위해 시뮬레이션을 수행한 것으로 보인다. 이는 위의 과정에서도 확인 할 수 있다.edges와 gwesp의 경우 p-value가 유의한 값을 가지고 있지만 gwdeg의 경우 그렇지 않다. 다른 변수를 적용하여 더 낮은 AIC를 가지는 모델은 없는지 찾아봐야할 것이다.
anova(karate.ergm.fit)## Analysis of Variance Table
##
## Model 1: karate.s ~ edges + gwdegree(0.2, fixed = T) + gwesp(0.2, fixed = T)
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 561 777.71
## Model 1: 3 361.71 558 416.00 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model1의 Deviance가 361.48로 변화한 것을 확인할 수 있다. 또한 p-value값이 매우 작게 나와 위 모델은 유의미한 모델이라고 해석할 수 있을 것이다. 이는 3개의 변수를 이용한 karate.ergm 모델이 karate 네트워크의 연결성을 잘 나타내는 강한 증거라고 볼 수 있다.다만 위에서 gwdeg의 p-value가 낮게 나왔기 때문에 이 변수를 제외한 모델은 어느정도의 성능을 보이는지 확인해보는 과정이 필요하다고 생각한다.
gof.karate.ergm=gof(karate.ergm.fit)
par(mfrow=c(2,2))
plot(gof.karate.ergm)
프린트 8p를 참고하여 해석한 결과는 아래와 같다. boxplot은 simulated 된 모델들이며 굵은 선은 original network 그래프를 나타낸다. degree,edge-wise shared partners,minimum geodesic distance plot을 확인하면 boxplot이 굵은 직선위에 적절하게 fitting되어 있는 것을 확인 할 수 있다. 몇몇 boxplot은 굵은 직선위에 있지 않지만 회색 선으로 나타난 신뢰구간 안에 분포하고 있다. 이러한 이유로 위 모델의 적합성은 꽤 적절한 것으로 보인다. model statisitcs plot의 경우 gwdeg가 original 값과 약간 다른 것으로 보인다. 이는 위에서 확인한 바와 같이 다른 2 변수에 비해 높은 p-value를 가지고 있기 때문일 것이다.