## ----setup, include=FALSE------------------------------------------------ library(ergm) library(knitr) knitr::opts_chunk$set(cache=F, comment=NA) ## ----eval=FALSE---------------------------------------------------------- ## install.packages('statnet') ## library(statnet) ## ----eval=FALSE---------------------------------------------------------- ## install.packages('ergm') # will install the network package ## install.packages('sna') ## ----eval=FALSE---------------------------------------------------------- ## update.packages('name.of.package') ## ----eval=FALSE---------------------------------------------------------- ## install.packages('coda') ## ----results='hide', message=FALSE--------------------------------------- library(statnet) ## ----results='hide', message=FALSE--------------------------------------- library(ergm) library(sna) library(coda) ## ----eval=FALSE---------------------------------------------------------- ## # latest versions: ergm 3.6.0 and network 1.13.0 (as of 4/03/2016) ## sessionInfo() ## ------------------------------------------------------------------------ set.seed(0) ## ------------------------------------------------------------------------ data(package='ergm') # tells us the datasets in our packages ## ------------------------------------------------------------------------ data(florentine) # loads flomarriage and flobusiness data flomarriage # Let's look at the flomarriage network properties par(mfrow=c(1,2)) # Setup a 2 panel plot (for later) plot(flomarriage, main="Florentine Marriage", cex.main=0.8) # Plot the flomarriage network summary(flomarriage~edges) # Look at the $g(y)$ statistic for this model flomodel.01 <- ergm(flomarriage~edges) # Estimate the model summary(flomodel.01) # The fitted model object ## ------------------------------------------------------------------------ summary(flomarriage~edges+triangle) # Look at the g(y) stats for this model flomodel.02 <- ergm(flomarriage~edges+triangle) summary(flomodel.02) ## ------------------------------------------------------------------------ class(flomodel.02) # this has the class ergm names(flomodel.02) # the ERGM object contains lots of components. ## ------------------------------------------------------------------------ flomodel.02$coef # you can extract/inspect individual components ## ------------------------------------------------------------------------ wealth <- flomarriage %v% 'wealth' # %v% references vertex attributes wealth summary(wealth) # summarize the distribution of wealth plot(flomarriage, vertex.cex=wealth/25, main="Florentine marriage by wealth", cex.main=0.8) # network plot with vertex size proportional to wealth summary(flomarriage~edges+nodecov('wealth')) # observed statistics for the model flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth')) summary(flomodel.03) ## ------------------------------------------------------------------------ data(faux.mesa.high) mesa <- faux.mesa.high ## ------------------------------------------------------------------------ mesa par(mfrow=c(1,1)) # Back to 1-panel plots plot(mesa, vertex.col='Grade') legend('bottomleft',fill=7:12,legend=paste('Grade',7:12),cex=0.75) ## ------------------------------------------------------------------------ fauxmodel.01 <- ergm(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T)) summary(fauxmodel.01) ## ------------------------------------------------------------------------ table(mesa %v% 'Race') # Frequencies of race mixingmatrix(mesa, "Race") ## ------------------------------------------------------------------------ summary(mesa ~edges + nodematch('Grade',diff=T) + nodematch('Race',diff=T)) ## ------------------------------------------------------------------------ data(samplk) ls() # directed data: Sampson's Monks samplk3 plot(samplk3) summary(samplk3~edges+mutual) ## ------------------------------------------------------------------------ sampmodel.01 <- ergm(samplk3~edges+mutual) summary(sampmodel.01) ## ------------------------------------------------------------------------ missnet <- network.initialize(10,directed=F) missnet[1,2] <- missnet[2,7] <- missnet[3,6] <- 1 missnet[4,6] <- missnet[4,9] <- missnet[5,6] <- NA summary(missnet) # plot missnet with missing edge colored red. tempnet <- missnet tempnet[4,6] <- tempnet[4,9] <- tempnet[5,6] <- 1 missnetmat <- as.matrix(missnet) missnetmat[is.na(missnetmat)] <- 2 plot(tempnet,label = network.vertex.names(tempnet),edge.col = missnetmat) summary(missnet~edges) summary(ergm(missnet~edges)) ## ------------------------------------------------------------------------ missnet_bad <- missnet missnet_bad[4,6] <- missnet_bad[4,9] <- missnet_bad[5,6] <- 0 summary(missnet_bad) summary(ergm(missnet_bad~edges)) ## ----eval=FALSE---------------------------------------------------------- ## help('ergm-terms') ## ------------------------------------------------------------------------ flomodel.03.sim <- simulate(flomodel.03,nsim=10) class(flomodel.03.sim) summary(flomodel.03.sim) length(flomodel.03.sim) flomodel.03.sim[[1]] plot(flomodel.03.sim[[1]], label= flomodel.03.sim[[1]] %v% "vertex.names") ## ------------------------------------------------------------------------ flo.03.gof.model <- gof(flomodel.03 ~ model) flo.03.gof.model plot(flo.03.gof.model) ## ------------------------------------------------------------------------ flo.03.gof.global <- gof(flomodel.03 ~ degree + esp + distance) flo.03.gof.global plot(flo.03.gof.global) ## ------------------------------------------------------------------------ mesamodel.b <- ergm(mesa~edges) plot(gof(mesamodel.b ~ model, nsim=10)) plot(gof(mesamodel.b ~ degree + esp + distance, nsim=10)) ## ------------------------------------------------------------------------ summary(flobusiness ~ edges+degree(1)) fit <- ergm(flobusiness ~ edges+degree(1)) mcmc.diagnostics(fit) ## ---- eval=FALSE--------------------------------------------------------- ## ergm(flobusiness ~ edges+degree(1), ## control=control.ergm(MCMC.interval=1)) ## ------------------------------------------------------------------------ data('faux.magnolia.high') magnolia <- faux.magnolia.high plot(magnolia, vertex.cex=.5) summary(magnolia ~ edges+triangle) ## ---- eval=F------------------------------------------------------------- ## ergm(magnolia ~ edges+triangle) ## ----eval=T-------------------------------------------------------------- fit.mag.01 <- ergm(magnolia ~ edges+triangle, control=control.ergm(MCMLE.maxit=2)) ## ---- eval=T, results='hide',fig.show='asis'----------------------------- mcmc.diagnostics(fit.mag.01) ## ------------------------------------------------------------------------ fit.mag.02 <- ergm(magnolia ~ edges+gwesp(0.25,fixed=T)) mcmc.diagnostics(fit.mag.02) ## ------------------------------------------------------------------------ gof.mag.02.model <- gof(fit.mag.02, GOF = ~model) gof.mag.02.model plot(gof.mag.02.model) ## ------------------------------------------------------------------------ fit.mag.03 <- ergm(magnolia ~ edges+gwesp(0.25,fixed=T) +nodematch('Grade')+nodematch('Race')+nodematch('Sex'), control = control.ergm(MCMC.interval=20000), eval.loglik = F) summary(fit.mag.03) ## ------------------------------------------------------------------------ mcmc.diagnostics(fit.mag.03) ## ------------------------------------------------------------------------ plot(gof(fit.mag.03, GOF=~model, control=control.gof.ergm(nsim=200))) ## ------------------------------------------------------------------------ plot(gof(fit.mag.03, GOF = ~ degree + esp + distance)) ## ------------------------------------------------------------------------ ego.net <- network.initialize(500, directed=F) ego.net %v% 'sex' <- c(rep(0,250),rep(1,250)) ## ------------------------------------------------------------------------ ego.deg <- c(180, 245, 60, 15) # node distn ego.mixmat <- matrix(c(164,44,26,176)/2, nrow=2, byrow=T) # adjusted tie distn ## ------------------------------------------------------------------------ ego.edges <- sum(ego.mixmat) ego.sexmatch <- ego.mixmat[1,1]+ego.mixmat[2,2] ## ------------------------------------------------------------------------ ego.target.stats <- c(ego.edges, ego.sexmatch) ego.target.stats ## ------------------------------------------------------------------------ ego.fit <- ergm(ego.net ~ edges + nodematch('sex'), target.stats = ego.target.stats) ## ------------------------------------------------------------------------ summary(ego.fit) ## ------------------------------------------------------------------------ ego.sim1 <- simulate(ego.fit) plot(ego.sim1, vertex.cex=.65, vertex.col="sex") ## ------------------------------------------------------------------------ rbind(sim=summary(ego.sim1 ~ degree(c(0:3))), obs=ego.deg) mixingmatrix(ego.sim1, "sex") ego.mixmat ## ------------------------------------------------------------------------ ego.sim100 <- simulate(ego.fit, nsim=100) ego.sim100 ## ------------------------------------------------------------------------ sim.stats <- attr(ego.sim100,"stats") rbind(sim=colMeans(sim.stats), obs=ego.target.stats) ## ------------------------------------------------------------------------ matplot(1:nrow(sim.stats), sim.stats, pch=c("e","m","0","+"), cex=.65, main="100 simulations from ego.fit model", sub="(default settings)", xlab="Replicate", ylab="frequency") abline(h=ego.target.stats, col=c(1:4)) ## ------------------------------------------------------------------------ sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10))) colnames(sim.fulldeg) <- paste("deg",0:10, sep='') sim.fulldeg[1:5,] ## ------------------------------------------------------------------------ sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum)) colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+") rbind(sim=colMeans(sim.deg), obs=ego.deg) ## ------------------------------------------------------------------------ matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5, main="Comparing ego.sims to non-targeted degree frequencies", sub = "(only total edges targeted)", xlab = "Replicate", ylab = "Frequencies") abline(h=c(180, 245, 60, 15), col=c(1:4)) ## ------------------------------------------------------------------------ ego.isolates <- ego.deg[1] ego.target.stats <- c(ego.edges, ego.sexmatch, ego.isolates) ego.fit <- ergm(ego.net ~ edges + nodematch('sex') + degree(0), target.stats = ego.target.stats) summary(ego.fit) ## ------------------------------------------------------------------------ ego.sim100 <- simulate(ego.fit, nsim=100, control=control.simulate.ergm(MCMC.interval=10000)) sim.stats <- attr(ego.sim100,"stats") rbind(sim=colMeans(sim.stats), obs=ego.target.stats) ## ------------------------------------------------------------------------ sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10))) sim.deg <- cbind(sim.fulldeg[,1:3], apply(sim.fulldeg[,4:11],1,sum)) colnames(sim.deg) <- c(colnames(sim.fulldeg)[1:3],"degree3+") rbind(sim=colMeans(sim.deg), obs=ego.deg) ## ------------------------------------------------------------------------ matplot(1:nrow(sim.deg), sim.deg, pch=as.character(0:3), cex=.5, main="Comparing ego.sims to non-targeted degree frequencies", sub = "(only 0, 2+ and total edges targeted)", xlab = "Replicate", ylab = "Frequencies") abline(h=c(180, 245, 60, 15), col=c(1:4))