
\documentclass{article}
\usepackage{natbib} %for easy biblo
\usepackage{hyperref} %for url links
\usepackage{comment} 
\usepackage{color} 
\usepackage{geometry}
\geometry{verbose,tmargin=3cm,bmargin=3cm,lmargin=3cm,rmargin=3cm}
\usepackage{multirow}
\setlength{\parskip}{\smallskipamount}

%\VignetteIndexEntry{ergm Example}

\begin{document}

\definecolor{Sinput}{rgb}{0.19,0.19,0.75}
\definecolor{Soutput}{rgb}{0.2,0.3,0.2}
\definecolor{Scode}{rgb}{0.75,0.19,0.19}
\DefineVerbatimEnvironment{Sinput}{Verbatim}{formatcom = {\color{Sinput}}} 
\DefineVerbatimEnvironment{Soutput}{Verbatim}{formatcom = {\color{Soutput}}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{formatcom = {\color{Scode}}} 
\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}

\SweaveOpts{concordance=TRUE}
\SweaveOpts{width=6,height=6}

\setkeys{Gin}{width=0.7\textwidth} 

<<foo,include=FALSE,echo=FALSE>>=

foo <- packageDescription("ergm")
@

<<echo=false>>=
options(width=75)
@

\title{Introduction to Exponential-family Random Graph\\
(ERG or p{*}) modeling with ergm\\
\small{Version \Sexpr{foo$Version}}  }

\author{The Statnet Development Team}

\maketitle

\tableofcontents

\section{Getting Started}
This vignette is based on the \verb@ergm@ tutorial presented at INSNA Sunbelt - St. Pete Beach, Florida, Feb 2011.

Open an R session, and set your working directory to the location
where you would like to save this work. You can do this with the pull-down
menus (File>Change Dir) or with the command:
<<eval=FALSE>>=
setwd('full.path.for.the.folder')
@

To install all of the packages in the statnet suite: 

<<eval=FALSE>>=
install.packages('statnet')
library(statnet)
@

Or, to only install the specific statnet packages needed for this
tutorial:

<<eval=FALSE>>=
install.packages('network')
install.packages('ergm')
install.packages('sna')
library(network)
library(ergm)
library(sna)
@


After the first time, to update the packages one can either repeat
the commands above, or use:

<<eval=FALSE>>=
update.packages('name.of.package')
@

For this tutorial, we will need one additional package (coda), which
is recommended (but not required) by ergm:

<<eval=FALSE>>=
install.packages('coda')
library(coda)
@

\section{Statistical network modeling; the \textit{ergm} command and \textit{ergm}
object}

Make sure the statnet package is attached:

<<eval=FALSE>>=
library(statnet)
@

or 

<<results=hide>>=
library(ergm)
library(sna)
set.seed(1)
@


The ergm package contains several network data sets that you can use
for practice examples.

<<fig=TRUE>>=
data(package='ergm') # tells us the datasets in our packages
data(florentine) # loads flomarriage and flobusiness data
flomarriage # Let's look at the flomarriage data
plot(flomarriage) # Let's view the flomarriage network
@

Remember the general ergm representation of the probability of the
observed network, and the conditional log-odds of a tie:

\[
\Pr(Y=y)=\exp[\theta'g(y)]/k(\theta)
\]

$Y$ is a network; $g(y)$ is a vector of network stats; $\theta$
is the vector of coefficients; $k(\theta)$ is a normalizing constant.

\[
\mbox{logit}\left(\Pr(Y_{ij}=1|Y^{c})\right)=\theta'\Delta\left(g(y)\right)_{ij}
\]

$Y_{ij}$ is an actor pair in $Y$; $Y^{c}$ is the rest of the network;
$\Delta\left(g(y)\right)_{ij}$ is the change in $g(y)$ when the
value of $Y_{ij}$ is toggled on.

We begin with the simplest possible model, the Bernoulli or Erd?s-R?nyi
model, which contains only an edge term.

<<>>=
flomodel.01 <- ergm(flomarriage~edges) # fit model
flomodel.01 

summary(flomodel.01) # look in more depth
@

How to interpret this model? The log-odds of any tie occurring is:
\begin{eqnarray*}
 &  & -1.609\times\mbox{change in the number of ties}\\
 & = & -1.609\times1
\end{eqnarray*}
for all ties, since the addition of any tie to the network changes
the number of ties by 1!

Corresponding probability is:
\begin{eqnarray*}
 &  & \exp(-1.609)/(1+\exp(-1.609))\\
 & = & 0.1667
\end{eqnarray*}
which is what you would expect, since there are 20/120 ties.

Let's add a term often thought to be a measure of
\textquotedblleft{}clustering\textquotedblright{}: the number of completed
triangles. Note we're in stochastic simulation now
\textendash{} your output will differ

<<>>=
flomodel.02 <- ergm(flomarriage~edges+triangle) 
summary(flomodel.02)
@

<<results=hide>>=
coef1 = flomodel.02$coef[1]
coef2 = flomodel.02$coef[2]
logodds = coef1 + c(0,1,2) * coef2
expit = function(x) 1/(1+exp(-x))
ps = expit(logodds)
coef1 = round(coef1, 3)
coef2 = round(coef2, 3)
logodds = round(logodds, 3)
ps = round(ps, 3)
@

Again, how to interpret coefficients?

Conditional log-odds of two actors forming a tie is:
\[
\Sexpr{coef1} \times \mbox{change in the number of ties} + \Sexpr{coef2} \times \mbox{change in number of triangles}
\]

\begin{itemize}
\item if the tie will not add any triangles to the network, its log-odds
is: $\Sexpr{coef1}$.
\item if it will add one triangle to the network, its log-odds is: $\Sexpr{coef1}+\Sexpr{coef2}=\Sexpr{logodds[2]}$
\item if it will add two triangles to the network, its log-odds is: $\Sexpr{coef1}+\Sexpr{coef2}\times2=\Sexpr{logodds[3]}$
\item the corresponding probabilities are \Sexpr{ps[1]}, \Sexpr{ps[2]}, and \Sexpr{ps[3]}.
\end{itemize}
Let's take a closer look at the ergm object itself:

<<>>=
class(flomodel.02) # this has the class ergm

names(flomodel.02) # let's look straight at the ERGM obj.
@

<<fig=TRUE>>=
flomodel.02$coef 
flomodel.02$formula 
flomodel.02$mle.lik
wealth <- flomarriage %v% 'wealth' # the %v% extracts vertex
wealth # attributes from a network
plot(flomarriage, vertex.cex=wealth/25) # network plot with vertex size 
                                        # proportional to wealth
@

We can test whether edge probabilities are a function of wealth:

<<>>=
flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth'))
summary(flomodel.03)
@

Yes, there is a significant positive wealth effect on the probability
of a tie.

\begin{center}
\rule{10cm}{0.4pt}
\end{center}

Let's try a model or two on:

Is there a statistically significant tendency for ties to be reciprocated ('mutuality')?

<<eval=FALSE>>=
data(samplk) 
ls() # directed data: Sampson's Monks
samplk3
plot(samplk3)
sampmodel.01 <- ergm(samplk3~edges+mutual)
summary(sampmodel.01)
@

Let's try a larger network

<<>>=
data(faux.mesa.high) 
mesa <- faux.mesa.high
@

<<eval=FALSE>>=
plot(mesa) 
mesa
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)
@

Note that two of the coefficients are estimated as -Inf
(the nodematch coefficients for race Black and Other). Why is this?

<<>>=
table(mesa %v% 'Race') # Frequencies of race
mixingmatrix(mesa, "Race")
@

So the problem is that there are very few students in the Black and
Other race categories, and these students form no homophilous (within-group)
ties. The empty cells are what produce the -Inf estimates.

Time to consider some missing data:

<<eval=FALSE>>=
missnet <- network.initialize(10,directed=F)
missnet[1,2] <- missnet[2,7] <- missnet[3,6] <- 1
missnet[4,6] <- missnet[4,9] <- NA
missnet
plot(missnet)
ergm(missnet~edges)
@

The coefficient equals -2.590. This is the log-odds of the probability
.0698. Our network has 3 ties, out of the 43 nodal pairs (10 choose
2 minus 2) whose dyad status we have observed. 3/43 = 0.0698.

<<eval=FALSE>>=
ergm(missnet~edges+degree(2))
missnet[4,6] <- missnet[4,9] <- 0
ergm(missnet~edges+degree(2))
@


The two estimates for the degree 2 coefficient differ considerably.
In the first case, there is one node we know for sure has degree 2,
two that may or may not, and seven that we know for sure do not. In
the latter, there is one node that has degree 2, and nine that do
not.


\section{Model terms available for \textit{ergm} estimation and simulation }

Model terms are the expressions (e.g. \textquotedblleft{}triangle\textquotedblright{})
used to represent predictors on the right-hand size of equations used
in:
\begin{itemize}
\item calls to \verb@ergm@ (to estimate an ergm model)
\item calls to \verb@simulate@ (to simulate networks from an ergm model
fit)
\item calls to \verb@summary@ (to obtain measurements of network statistics
on a dataset)
\end{itemize}

\subsection{Terms provided with ergm}

For a list of available terms that can be used to specify an ERGM,
see Appendix B, or type:

\verb@help('ergm-terms')@

For a more complete discussion of these terms see the 'Specifications'
paper in J Stat Software v. 24. (link is available online at \url{www.statnet.org})


\subsection{Coding new terms}

We have recently released a new package (\verb@ergm.userterms@)
and tutorial aimed at making it much easier than before to write one's
own terms. The package is available on CRAN, and installing it will
also download the tutorial (ergmuserterms.pdf). We teach a workshop
at the Sunbelt meetings, and are also hoping for the tutorial to appear
soon in the \textit{Journal of Statistical Software}. Note that writing
up new \verb@ergm@ terms requires some knowledge of C and the ability
to build R from source (although the latter is covered in the tutorial). 

\section{Network simulation: the \textit{simulate} command and \textit{network.list}
objects}

Once we have estimated the coefficients of an ERGM, the model is completely
specified. It defines a probability distribution across all networks
of this size. If the model is a good fit to the observed data, then
networks drawn from this distribution will be more likely to "resemble"
the observed data. To see examples of networks drawn from this distribution
we use the \verb@simulate@ command:

<<fig=TRUE>>=
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")
@

Voila. Of course, yours will look somewhat different.


\section{Examining the quality of model fit \textendash{} \textit{GOF} }

ERGMs are generative models \textendash{} that is, they represent
the process that governs tie formation at a local level. These local
processes in turn aggregate up to produce characteristic global network
properties, even though these global properties are not explicit terms
in the model. One test of whether a model \textquotedbl{}fits the
data\textquotedbl{} is therefore how well it reproduces these global
properties. We do this by choosing a network statistic that is not
in the model, and comparing the value of this statistic observed in
the original network to the distribution of values we get in simulated
networks from our model. 


<<fig=TRUE>>=
flomodel.03.gof <- gof(flomodel.03~degree)

flomodel.03.gof
plot(flomodel.03.gof)
@

<<fig=TRUE>>=
mesamodel.02 <- ergm(mesa~edges)
mesamodel.02.gof <- gof(mesamodel.02~distance,nsim=10)
plot(mesamodel.02.gof)

@

For a good example of model exploration and fitting for the Add Health
Friendship networks, see Goodreau, Kitts \& Morris, \textit{Demography}
2009.


\section{Diagnostics: troubleshooting and checking for model degeneracy}

The computational algorithms in \verb@ergm@ use MCMC to estimate
the likelihood function. Part of this process involves simulating
a set of networks to approximate unknown components of the likelihood. 

When a model is not a good representation of the observed network
the estimation process may be affected. In the worst case scenario,
the simulated networks will be so different from the observed network
that the algorithm fails altogether. This can occur for two general
reasons. First, the simulation algorithm may fail to converge, and
the sampled networks are thus not from the specified distribution.
Second, the model parameters used to simulate the networks are too
different from the MLE, so even though the simulation algorithm is
producing a representative sample of networks, this is not the sample
that would be produced under the MLE.

For more detailed discussions of model degeneracy in the ERGM context,
see the papers in \textit{J Stat Software} v. 24. (link is available
online at \url{www.statnet.org})

We can use diagnostics to see what is happening with the simulation
algorithm, and these can lead us to ways to improve it.

We will first consider a simulation where the algorithm works. To
understand the algorithm, consider

<<eval=FALSE>>=
fit <- ergm(flobusiness~edges+degree(1), 
  control=control.ergm(MCMC.interval=1, MCMC.burnin=1000, seed=1))
@

This runs a version with every network returned. Let us look at the
diagnostics produced:

<<eval=FALSE>>=
mcmc.diagnostics(fit, center=F)
@

Let's look more carefully at a default model fit:

<<eval=FALSE>>=
fit <- ergm(flobusiness~edges+degree(1))
mcmc.diagnostics(fit, center=F)
@

Now let us look at a more interesting case, using a larger network: 

<<>>=
data('faux.magnolia.high')
magnolia <- faux.magnolia.high
plot(magnolia, vertex.cex=.5)
@


<<eval=FALSE>>=
fit <- ergm(magnolia~edges+triangle, control=control.ergm(seed=1))
mcmc.diagnostics(fit, center=F)
@
<<echo=FALSE,eval=TRUE,results=verbatim,width=75>>= 
try_out <- try(ergm(magnolia~edges+triangle, control=control.ergm(seed=1))) 
cat(try_out) 
@

Very interesting. This model produced degenerate networks.
You could have gotten some more feedback about this
during the fitting, by using:

<<eval=FALSE>>=
fit <- ergm(magnolia~edges+triangle, control=control.ergm(seed=1), verbose=T)
@

You might try to increase the MCMC sample size:

<<eval=FALSE>>=
fit <- ergm(magnolia~edges+triangle,seed=1,
 control = control.ergm(seed=1, MCMC.samplesize=20000),
 verbose=T)
mcmc.diagnostics(fit, center=F)
@

How about trying the more robust version of modeling triangles:
GWESP? (For a technical introduction to GWESP see Hunter and Handcock;
for a more intuitive description and empirical application see Goodreau
Kitts and Morris 2009)

<<eval=FALSE>>=
fit <- ergm(magnolia~edges+gwesp(0.5,fixed=T),
 control =  control.ergm(seed=1))
mcmc.diagnostics(fit)
@

Still degenerate, but maybe getting closer?

<<eval=FALSE>>=
fit <- ergm(magnolia~edges+gwesp(0.5,fixed=T)+nodematch('Grade')+nodematch('Race')+
  nodematch('Sex'),
 control = control.ergm(seed=1),
 verbose=T)

pdf('diagnostics1.pdf')	#Use the recording function if possible, otherwise send to pdf
mcmc.diagnostics(fit)
dev.off()			#If you saved to pdf, look at the file

fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T)+nodematch('Grade')+nodematch('Race')+
	nodematch('Sex'),
 control = control.ergm(seed=1))
mcmc.diagnostics(fit)

@

One more try...
<<results=hide>>=
fit <- ergm(magnolia~edges+gwesp(0.25,fixed=T)+nodematch('Grade')+nodematch('Race')+
  nodematch('Sex'),
 control = control.ergm(seed=1,MCMC.samplesize=4096,MCMC.interval=8192),
 verbose=T)
@

<<fig1,fig=TRUE>>=
mcmc.diagnostics(fit)
@

Success! Of course, in real life one might have a lot more trial and
error.

Changes in version 3.2 of the ergm estimation algorithm mean that the MCMC diagnostic plots can no longer be used to ensure that the mean statistics from the model match the observed network statistics. For that functionality, please use the GOF command: \verb@gof(fit, GOF=~model)@. The plots can still be used to assess MCMC mixing and convergence.


\section{Working with egocentrically sampled network data}

In many empirical contexts, it is not feasible to collect a network
census or even an adaptive (link-traced) sample. Even when one of
these may be possible in practice, egocentrically sampled data are
typically cheaper and easier to collect. 

Long regarded as the poor country cousin in the network data family,
egocentric data contain a remarkable amount of information. With the
right statistical methods, such data can be used to explore the properties
of the complete networks in which they are embedded. The basic idea
here is to combine what is observed, with assumptions, to define a
class of models that describe the distribution of networks that are
centered on the observed properties. The variation in these networks
quantifies some of the uncertainty introduced by the assumptions.

Let's start with a simple fictional example: You have
a sample of persons who were asked about the friends they had seen
face-to-face more than once in the last week. The respondent was asked
their own sex, and the sex of each friend (for up to three friends).
Summary statistics from these data thus include the sex distribution,
the degree distribution (it could be broken down by sex, but we will
just examine the marginal distribution here), and the joint distribution
of the respondent and friend's sex (the sex mixing
matrix). Let's assume there are equal numbers of men
and women in the sampled respondents. The other distributions are
shown below:

\noindent \begin{center}
\begin{tabular}{cccc}
\multicolumn{4}{c}{\textbf{Degree distribution}}\tabularnewline
\textbf{Degree} & \textbf{Frequency} & \textbf{Fraction} & \textbf{Ties}\tabularnewline
\hline 
0 & 180 & 0.36 & 0\tabularnewline
1 & 245 & 0.49 & 245\tabularnewline
2 & 60 & 0.12 & 120\tabularnewline
3 & 15 & 0.03 & 45\tabularnewline
\hline 
Total & 500 & 1.00 & 410\tabularnewline
\end{tabular}\qquad{}%
\begin{tabular}{cc|c|c}
\multicolumn{4}{c}{\textbf{Sex mixing matrix }}\tabularnewline
(410 total) & \multicolumn{1}{c}{} & \multicolumn{2}{c}{\textbf{Friend}}\tabularnewline
\multirow{3}{*}{\textbf{Respondent}} & \multicolumn{1}{c}{} & \multicolumn{1}{c}{M} & F\tabularnewline
\cline{3-4} 
 & M & 164 & \multicolumn{1}{c|}{44}\tabularnewline
\cline{3-4} 
 & F & 26 & \multicolumn{1}{c|}{176}\tabularnewline
\cline{3-4} 
\end{tabular}
\par\end{center}

So, total N respondents = 500, total N friends reported = 410. 

We can use an ERGM to fit the parameters associated with these observed
statistics, then use the fitted model to simulate complete networks
that are drawn from the distribution of networks that is centered
around these statistics. As a theoretical exercise, this provides
a method for investigating the complete network implications of these
observed summary statistics. As an empirical exercise, the primary
assumption needed for inference is that the data we have is sampled
from a population in equilibrium (and, as in all statistical inference,
that our model is correct). The theory for this is developed in Krivitsky,
2009.

We need to make assumptions about size, directedness and bipartite-ness
when we model and simulate the complete network.
\begin{itemize}
\item \textbf{Size}: any size can be simulated, but if the model is fit
using the observed frequencies, it should be used to simulate a population
of that size, unless a size adjustment is made in the simulation (see
Krivitsky, Handcock and Morris 2011). We are going to work with a
population size 500 here, equal to the number of respondents.
\item \textbf{Directedness}: Ego data are in one sense inherently directed
(respondents nominate alters, alters are not observed), but the relationship
may be either. In this case (\textquotedblleft{}seen more than once\textquotedblright{})
it is undirected, so we will fit and simulate an undirected network.
\item \textbf{Bipartite}: Ego data can be bipartite (if no alters are also
respondents, or data are collected on 2-mode networks) or not (if
respondents are also alters). But again, the relationship may be either.
\textquotedblleft{}Seen\textquotedblright{} is undirected, and we
will fit and simulated and undirected network.
\end{itemize}
In sum, we will simulate a one-mode, undirected network of size 500,
assuming the ego statistics we observed contain the information we
need to calculate the statistics that would have been observed in
a self-contained population of that size, noting that other assumptions
are possible.

To ensure consistency between the degree distribution (which is a
tabulation of nodes) and the mixing matrix (which is a cross-tabulation
of ties) in our simulated \textquotedblleft{}complete network,\textquotedblright{}
it is important to recognize that in a complete network, the degree
distribution should imply twice the number of ties observed in the
mixing matrix, because every tie is being reported by both nodes in
the degree distribution. If we are fixing the population size at 500
in this simulation, then our observed mixing matrix data needs to
be divided by 2.

Start by initializing an empty network of the desired size and assign
the \textquotedblleft{}sex\textquotedblright{} attribute to the nodes:

<<>>=
ego.net <- network.initialize(500, directed=F)
ego.net %v% 'sex' <- c(rep(0,250),rep(1,250))
@

Set up the observed statistics (adjusted for this \textquotedblleft{}complete\textquotedblright{}
network) as we will use them to assess the accuracy of the simulation
later:

<<>>=
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
@

Then, pick the observed statistics you want to target \textendash{}
we will start with a simple model here: the total number of ties (edges),
and the number of sex-matched ties (homophily). These are both functions
of the observed statistics:

<<>>=
ego.edges <- sum(ego.mixmat)
ego.sexmatch <- ego.mixmat[1,1]+ego.mixmat[2,2]
@

And combine these into a vector 

<<>>=
ego.target.stats <- c(ego.edges, ego.sexmatch)
ego.target.stats
@

Now, fit an ERGM to this \textquotedblleft{}network\textquotedblright{},
with terms for the statistics you want to match, and the \textquotedblleft{}\texttt{target.stats}\textquotedblright{}
argument for \textbf{ergm} that specifies the target values for those
statistics:

<<>>=
ego.fit <- ergm(ego.net ~ edges + nodematch('sex'),
 target.stats = ego.target.stats)
@

Take a look at the fitted model:

<<>>=
summary(ego.fit) 
@

Now that you have a fitted model, you can simulate a complete network
from it, and look at the results:

<<fig=TRUE>>=
ego.sim1 <- simulate(ego.fit)
plot(ego.sim1, vertex.cex=.65, vertex.col="sex")
@

Does it reproduce the observed degree and mixing frequencies? We only
targeted the total number of edges, not the full degree distribution.

<<>>=
rbind(summary(ego.sim1 ~ degree(c(0:3))), ego.deg)
mixingmatrix(ego.sim1, "sex")
ego.mixmat
@

We only targeted the number of same-sex ties, not the full mixing
matrix.

The simulation stats seem quite different than the observed stats,
and there are two possible reasons: either we mis-specified the original
model (bias), or this one random draw may be unrepresentative of the
distribution described by the model (variance). The latter is easily
examined by simulating 100 networks, to see where the observed data
fall in the distribution of networks produced by the model:

<<>>=
ego.sim100 <- simulate(ego.fit, nsim=100)
ego.sim100
@

More information can be obtained with 
<<eval=FALSE>>=
summary(ego.sim100)
@

First, we\textquoteright{}ll look at how well the simulations reproduced
the target statistics that were included in the model (note, not the
observed full distributions):

<<>>=
sim.stats <- attr(ego.sim100,"stats")
rbind(colMeans(sim.stats), ego.target.stats)
@

These look pretty good \textendash{} the means of the simulated target
stats are within 1\% of the observed. We can plot the 100 replicates
to see check the variation for any problematic patterns:

<<fig=TRUE>>=
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))
@

The lines mark the target statistic frequencies in the observed data.
The points represent the frequencies in the simulated networks.

The simulated network statistics vary stochastically around the target
values, not trending over time.

But, there is clear autocorrelation across the replicates, which suggests
we might want to increase the MCMC interval to draw more independent
realizations.

<<fig=TRUE>>=
ego.sim100 <- simulate(ego.fit, nsim=100,
  control=control.simulate.ergm(MCMC.interval=10000))
sim.stats <- attr(ego.sim100,"stats")
matplot(1:nrow(sim.stats), sim.stats,
  pch=c("e","m"), cex=.65,
  main="100 simulations from ego.fit model", sub="(MCMC.interval=10000)",
  xlab="Replicate", ylab="frequency")
abline(h=ego.target.stats, col=c(1:2))
@

With the larger interval, the autocorrelation is no longer detectable,
and all of the statistics from the simulated networks vary in a symmetric
band around their targets.

The variation (about +/- 10\%) represents the range of target statistics
that are consistent with the fitted coefficients.

If you wanted instead to constrain these statistics to equal a specified
value, then you can use the \textquotedblleft{}constraints=\textquotedblright{}
argument during the \verb@ergm@ fit instead. 

This is good for verifying that the simulation reproduces the target
values we specified. So now let\textquoteright{}s see whether the
simulated complete networks also match statistics that were not set
by the targets. We targeted edges and homophily. How well does this
model reproduce the full degree distribution?

<<>>=
sim.fulldeg <- summary(ego.sim100 ~ degree(c(0:10)))
colnames(sim.fulldeg) <- paste("deg",0:10, sep='')
sim.fulldeg[1:5,]
@

Recall that the degree range in our data was {[}0,3{]} by design,
but we did not constrain the simulations to this range. If our model
correctly captured the processes that led to the aggregate statistics
we observe in our data, the simulated networks would show us what
we missed. Here the simulated networks suggest that the fully observed
network would have a wider range of degrees, which we might have observed,
had we not truncated the data collection at 3 friends per respondent.
About 1\% of nodes have degree 4 or 5, and the max observed is 6.

But did our model did correctly capture the underlying processes?
How well does the simulated degree distribution from this model match
the frequencies we did observe? Aggregating the degrees of 3 or more
in the simulations, we find:

<<>>=
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(colMeans(sim.deg),ego.deg)
@

As with the single simulation above, the discrepancies between the
simulation averages and the observed statistics are quite large. We
can see this more clearly by plotting the degree frequencies for the
100 replicate networks we simulated:

<<fig=TRUE>>=
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))
@

The simulations are producing systematically more isolates than expected
(the \textquotedblleft{}0\textquotedblright{} points vs. the black
line), and systematically more degree 1 nodes. In fact, the two degree
frequencies are essentially reversed in the simulation. 

The fraction of nodes with 2 or 3+ partners is systematically off
but by a much smaller amount.

So our observed network has fewer isolates than expected in a network
of this density, more degree 1 nodes than expected, and fewer degree
2 and 3+ nodes.

This suggests the model is mis-specified. Since the degree 0 vs. degree
1 is the worst fitting aspect, we will try using the number of isolates
as a target statistic in the model.

<<>>=
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)
@

Simulating from this model, the target statistics are again well matched:

<<>>=
ego.sim100 <- simulate(ego.fit, nsim=100,
   control=control.simulate.ergm(MCMC.interval=10000))
sim.stats <- attr(ego.sim100,"stats")
rbind(colMeans(sim.stats), ego.target.stats)
@

And the full degree frequencies look much better:

<<>>=
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(colMeans(sim.deg),ego.deg)
@

and finally, plotting the full degree frequencies

<<fig=TRUE>>=
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))
@


The degree frequencies in these simulated networks are now well centered
on the observed frequencies. So adding the one additional parameter
to capture the lower than expected number of isolates did a good job
of capturing how our observed network deviates from a random network
with this density. 

The fraction of nodes with 3+ partners produced by our new model might
still be a bit low.

\textbf{Moral}: We can use ERGMs to estimate network models using
target statistics from egocentrically sampled data. The fact that
the target statistics are reproduced by this model does not guarantee
that additional features of the network would also be reproduced.
But starting with simple models can help to identify whether and how
the aggregate statistics we observe from an egocentric sample deviate
from those we would expect from the model. If we fit all of the observed
statistics without a saturated model, we cannot reject the hypothesis
that this model produced the network we sampled from.

We can also use this approach to explore network statistics that are
not visible at all from the egocentric data, e.g., the geodesic distribution,
betweenness, etc., but it must always be remembered that the distributions
we will produce are based on our model. They faithfully reproduce
the model, but that does not mean that the model faithfully represents
the population.

In the STERGM workshop, we show how complete dynamic networks also
can be simulated over time on the basis of egocentric data like these,
with the minimal addition of a single estimate of partnership duration.
For a movie of an example dynamic simulation used to explore the impact
of network structure on HIV transmission, 
see \href{http://statnet.org/movies}{statnet.org/movies}.


\section{Additional functionality in the statnet family of packages}


\subsection{In the \texttt{ergm} and \texttt{network} packages:}

\begin{itemize}

\item ERGMs for valued ties \textendash{} see the paper by Pavel Krivitsky (2012)

\item Analysis of bipartite networks \textendash{} \verb@network@ recognizes this as an attribute of the network and \verb@ergm@ provides specific model terms for such networks that begin with b1 or b2 \\
(try: \texttt{search.ergmTerms(categories=c('bipartite'))}).

\end{itemize}

\subsection{In other packages from the statnet suite:}

These are in stand-alone packages that can be downloaded either from CRAN, 
or from the \verb@statnet@ website. Many have online training materials from our
workshops.  For more detailed information, 
please visit the \verb@statnet@ webpage 
(\href{http://statnet.org}{http://statnet.org}).

\subsubsection{Static (cross-sectional) network analysis packages}

\begin{itemize}

\item \href{http://cran.r-project.org/web/packages/sna/index.html}{\texttt{sna}} \textendash{} Traditional SNA methods and summaries.

\item \href{http://cran.r-project.org/web/packages/latentnet/index.html}{\texttt{latentnet}}
\textendash{} Latent space and latent cluster analysis.

\item \href{http://cran.r-project.org/web/packages/netperm/index.html}{\texttt{netperm}} 
\textendash{} Network permutation models.

\item \href{http://cran.r-project.org/web/packages/degreenet/index.html}{\texttt{degreenet}} 
\textendash{} MLE estimation for degree distributions (negative binomial, Poisson, scale-free, etc.) 

\item \href{http://cran.r-project.org/web/packages/networksis/index.html}{\texttt{networksis}}
\textendash{} Simulation of bipartite networks with given degree distributions. 

\end{itemize}

\subsubsection{Dynamic (longitudinal) network analysis packages:}

\begin{itemize}

\item \href{http://cran.r-project.org/web/packages/tergm/index.html}{\texttt{tergm}}
\textendash{} Temporal ERGMs (TERGMs): discrete-time dynamic network models for longitudinal network panel data, and other temporal extensions.

\item \href{http://cran.r-project.org/web/packages/relevent/index.html}{\texttt{relevent}}
\textendash{} Relational event models: continuous-time dynamic network models for longitudinal network data.

\item \href{http://cran.r-project.org/web/packages/networkDynamic/index.html}{\texttt{networkDynamic}}
\textendash{} Dynamic network data storage and manipulation.

\item \href{http://cran.r-project.org/web/packages/ndtv/index.html}{\texttt{ndtv}}
\textendash{} Network movie maker.

\end{itemize}


\subsection{R packages that build on statnet}

There is a growing number of R packages written by other folks that build on or extend the functionality of the statnet suite.  You can get a current list of those packages by looking at the reverse depends/suggests on CRAN.  A partial list includes:

\begin{itemize}

\item \href{http://cran.r-project.org/web/packages/EpiModel/index.html}{\texttt{EpiModel} package}
\textendash{} Mathematical Modeling of Infectious Disease, includes functions for deterministic compartmental modeling, stochastic individual contact modeling, and stochatic network modeling

\item \href{http://cran.r-project.org/web/packages/RDS/index.html}{\texttt{RDS} package}
\textendash{} Estimation with data collected using Respondent-Driven Sampling.

\item \href{http://cran.r-project.org/web/packages/Bergm/index.html}{\texttt{Bergm} package} \textendash{} Bayesian ERGM estimation 

\item \href{http://cran.r-project.org/web/packages/hergm/index.html}{\texttt{hergm} package} \textendash{} Hierarchical Exponential-Family Random Graph Models with Local Dependence (for latent groups).

\item \href{http://cran.r-project.org/web/packages/lvm4net/index.html}{\texttt{lvm4net} package} 
\textendash{} Latent variable models.

\item \href{http://cran.r-project.org/web/packages/VBLPC/index.html}{\texttt{VBLPCM} package}
\textendash{} Variational Bayes Latent Position Cluster Models. 

\item \href{http://cran.r-project.org/web/packages/xergm/index.html}{\texttt{xergm} package}
\textendash{} Temporal exponential random graph models (TERGM) by bootstrapped pseudolikelihood, MCMC MLE and (temporal) network autocorrelation models.

\end{itemize}

\subsection{Additional functionality in development:}
\begin{itemize}

\item \texttt{ergm.ego} package \textendash{} ERGM estimation and inference from egocentrically sampled data (expected May 2015)

\item \texttt{tsna} package \textendash{} Temporally extended (vertex and edge) SNA methods for dynamic longitudinal network data (expected May 2015)

\item \texttt{MLergm} package \textendash{} ERGM estimation and inference for multi-level data (for observed groups) (expected 2016)

\end{itemize}

\section{Statnet Commons: The core development team}

  Mark S. Handcock <\href{mailto:handcock@stat.ucla.edu}{handcock@stat.ucla.edu}>\\
  David R. Hunter <\href{mailto:dhunter@stat.psu.edu}{dhunter@stat.psu.edu}>\\
  Carter T. Butts <\href{mailto:buttsc@uci.edu}{buttsc@uci.edu}>\\
  Steven M. Goodreau <\href{mailto:goodreau@u.washington.edu}{goodreau@u.washington.edu}>\\
  Skye Bender-deMoll <\href{mailto:skyebend@skyeome.net}{skyebend@skyeome.net}>\\
  Martina Morris <\href{mailto:morrism@u.washington.edu}{morrism@u.washington.edu}>\\
  Pavel N. Krivitsky <\href{mailto:pavel@UOW.EDU.AU}{pavel@uow.edu.au}>
  
\section*{Appendix A: Clarifying the terms \textendash{} ergm and network }

You will see the terms ergm and network used in multiple contexts
throughout the documentation. This is common in R, but often confusing
to newcomers. To clarify:

\textbf{ergm} 
\begin{itemize}
\item \textbf{ERGM}: the acronym for an Exponential Random Graph Model;
a statistical model for relational data that takes a generalized exponential
family form. 
\item \textbf{ergm package}: one of the packages within the \verb@statnet@
  suite
\item \textbf{ergm function}: a function within the ergm package; fits an
ERGM to a network object, creating an ergm object in the process. 
\item \textbf{ergm object}: a class of objects produced by a call to the
ergm function, representing the results of an ERGM fit to a network. 
\end{itemize}

\textbf{network} 
\begin{itemize}
\item \textbf{network}: a set of actors and the relations among them. Used
interchangeably with the term graph. 
\item \textbf{network package}: one of the packages within the \verb@statnet@
  suite; used to create, store, modify and plot the information found
in network objects. 
\item \textbf{network object}: a class of object in \verb@R@ used to represent
a network.
\end{itemize}


\begin{thebibliography}{}

\bibitem[Goodreau, S., J. Kitts and M. Morris (2009)]{goodreau}
Goodreau, S., J. Kitts and M. Morris (2009).
\newblock Birds of a Feather, or Friend of a Friend? Using Statistical Network Analysis to Investigate Adolescent Social Networks.
\newblock \emph{Demography 46(1)}:103-125. 

\bibitem[Handcock, M. S., D. R. Hunter, C. T. Butts, S. M. Goodreau and M. Morris (2008)]{handcock}
Handcock, M. S., D. R. Hunter, C. T. Butts, S. M. Goodreau and M. Morris (2008).
\newblock statnet: Software Tools for the Representation, Visualization, Analysis and Simulation of Network Data.
\newblock \emph{Journal of Statistical Software 42(01)}.

\bibitem[Hunter et al.(2008b)]{ergm}
Hunter DR, Handcock MS, Butts CT, Goodreau SM, Morris M (2008b). 
\newblock ergm: A Package to Fit, Simulate and Diagnose Exponential-Family Models for Networks. 
\newblock \emph{Journal of Statistical Software}, 24(3). \url{http://www.jstatsoft.org/v24/i03/}. 

\bibitem[Krivitsky, P.N.(2009)]{pavelthesis}
Krivitsky, P.N.(2009). 
\newblock PhD thesis.
\newblock \emph{University of Washington, Seattle, WA}

\bibitem[Krivitsky, P. N., M. S. Handcock and M. Morris (2011).]{pavelSize}
Krivitsky, P. N., M. S. Handcock and M. Morris (2011). 
\newblock Network Size and Composition Effects in Exponential-Family Random Graph Models.
\newblock \emph{Statistical Methodology 8:}319-339

\bibitem[Krivitsky, P. N. (2012).]{pavelValued}
Krivitsky PN (2012). 
\newblock Exponential-family random graph models for valued networks.
\newblock \emph{Electronic Journal of Statistics 6:}1100-1128

\end{thebibliography}


\end{document}
