ergm.Rmd

---
title: "Fitting Exponential Random Graph Models with `statnet`"
bibliography: references.bib
---

```{r setup, echo=FALSE, results="hide", cache=FALSE}
suppressMessages({
  library(knitr)
  library(methods)
  } )

pdf.options(
  family="Palatino"
)

opts_chunk$set(
  fig.retina=NULL,
  fig.align="center",
  cache=TRUE
  )

.small_mar <- function(before, options, envir) {
  # save current
  op <- par(no.readonly=TRUE)
  if(before) {
    par(mar=c(1, 1, 1, 3))
  } else {
    # restore
    par(op)
  }
}

.no_mar <- function(before, options, envir) {
  # save current
  op <- par(no.readonly=TRUE)
  if(before) {
    par(mar=rep(0.5, 4))
  } else {
    # restore
    par(op)
  }
}

.use_igraph <- function(before, options, envir) {
  if(before) {
    detach("package:network", force=TRUE)
    library(igraph)
    igraph_options(
      vertex.color="lightskyblue",
      edge.arrow.size=0.5,
      vertex.size=25,
      vertex.label.cex=0.7,
      vertex.label.family="Palatino"
    )
  } else {
    detach("package:igraph", force=TRUE)
    library(network)
  }
}

knit_hooks$set( 
  small.mar = .small_mar, 
  no.mar=.no_mar, 
  use.igraph=.use_igraph
  )

set.seed(123)
```




# Introduction

This tutorial supplements the topics covered on WSAD Summer School workshop
"Introduction Statistical Social Network Analysis". The workshop presented how
Exponential Random Graph Models (ERGMs) can be fit using program PNet
[@wang_etal_2009]. Here we provide examples how ERGMs can be simulated and
fitted using [`statnet` suite](http://www.statnet.org) of R packages for network analysis
[@handcock_etal_2003].

In this tutorial we will use network data in the form of objects of class
`network` instead of objects of class `igraph`.  In principle, you can use both
packages at the same time. However, we recommend to use only one at a time to
avoid function name conflicts [^1]. You can easily convert `igraph` objects to
`network` objects (or vice versa) using `asIgraph` and `asNetwork` functions
from `intergraph` package.

[^1]: For more details see for example
http://bc.bojanorama.pl/2010/08/namespaces-and-name-conflicts/





# Exponential Random Graph Models in a nutshell

ERG models, as they are currently known, stem from the work by
@holland_leinhardt_1976 and @frank_strauss_1986, and were further developed by,
among others, @frank_1991, @wasserman_pattison_1996, @pattison_robins_2002,
@snijders_etal_2006. ERGMs have established their position as one of the most
important tools in social networks analysis. Comprehensive description of ERG
models is provided in the book by @lusher_etal_2012.

Exponential Random Graph Models are a family of statistical models for
understanding processes that shape the global structure of a network. This goal
is achieved by assigning probability to networks according to network
statistics -- summary measures describing selected features of the network.

Formally, let $A$ represent a network in the form of an adjacency matrix of
given size $n$. Let $\mathcal{A}$ be a collection of all possible networks of
size $n$ that can be constructed. An ERG model postulates that a probability
distribution defined over all networks in $\mathcal{A}$ can be represented as:

$$
P_\theta(A=a)=\frac{\exp{(\theta^Tg(a))}}{\kappa(\theta, \mathcal{A})}, \qquad a \in \mathcal{A},
$$

where

- $g(a)$ is a vector of network statistics, 
- $\theta$ is a vector of weights associated with these statistics, 
- $\kappa(\theta, \mathcal{A})$ is a normalizing constant ensuring that all the probabilities sum-up to 1.

$$
\kappa(\theta, \mathcal{A}) = \sum_{a \in \mathcal{A}}\exp{(\theta^Tg(a))}.
$$

In words, the probability that we observe a particular network $a$ out of the
set $\mathcal{A}$ of all possible networks of size $n$ is a function of chosen
network statistics $g(a)$, associated parameters $\theta$, and the normalizing
constant $\kappa$.

The model can also be written as a model for conditional log-odds of tie
existence between a pair of nodes.  This makes it a little bit similar to a
logistic regression model (see section [Dyadic independence
ERGMs](#dyadic_independence) for more details):

$$
\operatorname{logit}{ (Y_{ij} = 1 | y^c_{ij} )  = \theta \delta( y_{ij} ) }
$$

where:

- $Y_{ij}$ is a binary variable capturing the state of a pair of actors $i$ and $j$, whether there is a tie ($Y_{ij}=1$) or not ($Y_{ij}=0$).
- $y_{ij}$ is a particular realization of variable $Y_{ij}$.
- $y^c_{ij}$ is the complement of $y_{ij}$, i.e. all other actor pairs in the network apart from $(i, j)$.
- $\delta(y_{ij})$ is a vector of *change statistics* for each term in the model.
- $\theta$ is a vector of parameters, as before.

In words, the log-odds that actors $i$ and $j$ are connected with a tie, given
the particular configuration of all the remaining ties, is a function of change
statistics for each term in the model.

The change statistics capture how a particular term $g(y)$ changes if the
$y_{ij}$ tie is formed or not:

$$
\delta( y_{ij} ) = g( y^{+}_{ij} ) - g( y^{-}_{ij} )
$$

where $y^{+}_{ij}$ represents $y^c$ and $y_{ij}$ set to 1, while $y^{-}_{ij}$
represents $y^c$ and $y_{ij}$ set to 0.

As a consequence, model parameters $\theta$ can be interpreted as log-odds of a
presence of a network tie conditional on the state of all other dyads in the
network.






## Network statistics

Different ERG models can be specified by choosing different network statistics.
A set of network statistics is chosen on the basis of theoretical premises
concerning a particular research problem at hand. However, there are a few
basic configurations often included in the model. Example of network statistics
include:

1. Number of edges in the network (`edges`)
$$
S_1(a) = \sum_{1\leq i\leq j\leq n}a_{ij},
$$

2. Number of _k_-stars (for $k\geq2$). A _k_-star consists of a central node and $k$ neighbors. (`kstar`).
$$
S_k(a) = \sum_{1\leq i\leq n}\binom{k_i}{k}.
$$

3. Number of triangles (`triangles`).
$$
T(a) = \sum_{1\leq i < j < h\leq n}a_{ij}a_{ih}a_{jh}.
$$

ERGMs could easily incorporate additional information about actors attributes.
To indicate extra information function $g(a)$ could be replaced with $g(a,X)$,
where $X$ is a matrix containing attributes. This leads to another statistics,
which measure the effect of actor's attributes. $y_i$ denotes value of an
attribute $y$ of the $i$-th node:

4. Attribute-based activity (`nodeofactor`)
$$
\sum_{1\leq i\leq j\leq n}a_{ij}(y_i + y_j),
$$

5. Homophily (binary attribute) (`nodematch`).
$$
\sum_{1\leq i\leq j\leq n}a_{ij}y_i  y_j,
$$

6. Homophily (continuous attribute) (`absdif`)
$$
\sum_{1\leq i\leq j\leq n}a_{ij}|y_i-  y_j|.
$$

See `?ergm.terms` for a complete list of network statistics implemented in
`ergm` or view it on the www.rdocumentation.org [here](http://www.rdocumentation.org/packages/ergm/functions/ergm-terms).






## Calculating network statistics

The workhorse for estimating ERG models is the `ergm` function from the `ergm`
package. We will describe it in detail later. Function `ergm` uses a formula
interface to specify the network statistics, or `terms`, to be included in the model.
If you simply want to calculate the values of network statistics in a given
network you can use the function `summary` with a single argument -- an ERGM
formula. Such formula should have an object of class `network` on the left hand
side, and a sum of network statistics on the right hand side.

As an example, let us take the classroom network we used in other tutorials. We
need it now as a `network` object:

```{r, classroom_as_network, no.mar=TRUE}
library(network)
# load data
data(IBE121, package="isnar")
# select the "would like to play with" network
# by dropping other edges
playnet <- igraph::delete.edges(IBE121, igraph::E(IBE121)[question != "play"])
# convert to 'network' object
ibe <- intergraph::asNetwork(playnet)
ibe
plot(ibe, vertex.col=ifelse(ibe %v% "female", "pink", "lightskyblue"))
legend("topright", pch=21, pt.bg=c("pink", "lightskyblue"), col="black",
       legend=c("Girl", "Boy"), title="Gender", bty="n")
```

We can now use `summary` to calculate some network statistics of interest. For
example, let's calculate

- number of ties
- number of reciprocated ties
- Gender-based activity (out-degree)


```{r}
library(ergm)
summary(ibe ~ edges + mutual + nodeofactor("female"))
```

```{r, silent_statistics, results="hide", echo=FALSE}
s <- summary(ibe ~ edges + mutual + nodeofactor("female"))
```

Consequently, this network has `r s[1]` ties, `r s[2]` of which are reciprocated, and girls
send ties `r s[3]` times.






## Fitting ERGMs

Before we go into more complex examples let us illustrate how ERG models are
fit with `statnet` using two very simple models.

### Homogenous Bernoulli model       {#model_bernoulli}

First let's estimate the simplest possible model. The model assumes that the
network is purely random: every tie in the network appears with the same
unknown probability $p$. It is also known as *homogeneous Bernoulli model* as
each tie is sampled independently with a given probability $p$. Appropriate
ERGM formulation is:

```{r, bernoulli}
model.ibe0 <- ergm(ibe ~ edges)
summary(model.ibe0)
```

The model consists of a single term `edges`: number of edges in the network.

The `ibe` network has 26 nodes. This means that the number of all possible ties
in this network is equal to $26 \times 25 = 650$. Odds for tie existence are
then equal to $88 / (650 - 88) = `r (o <- 88/(650-88))`$.  If we calculate log-odds we will
get $\log (0.156583) = `r log(o)`$ which is exactly the value of the `edges` term
in the model above. This corresponds to a constant tie probability of $88/650 = `r 88/650`$.
This, in turn, is obviously equal to network density:

```{r, ibe_density}
network.density(ibe)
```

As network density is usually not of our primary research interest, the `edges`
term is usually included in models and interpreted in a similar manner like the
intercept in linear regression.



### Are ties reciprocated?       {#model_reciprocity}

As a little bit more complex example, let's assess the extent, to which
children nominations are reciprocated, i.e., if A nominated B is it more likely
for B to nominate A back, instead of nominating someone else.  We need a model
with two terms:

- `edges` which is the number of edges in the network. With this term we model overall network density.
- `mutual` which is the number of symmetrical (reciprocated) dyads in the network.

To fit the model to our classroom network we use the function `ergm` and
provide a model formula as an argument. The formula has a network object on the
left hand side and a sum of the two terms on the right hand side:

```{r, reciprocity}
model.ibe1 <- ergm(ibe ~ edges + mutual, control=control.ergm(seed=666))
summary(model.ibe1)
```


From the results we see that the `edges` effect is negative and significant
while the `mutual` term is positive and significant too. Both effects are
roughly equal in absolute size, so cancel each other out, which means that odds
for a mutual tie are about 1 and the probability about 0.5. Conditional odds
for a non-mutual tie are $\exp(-2.466) = 0.0849$ so the probability about 0.08.
Indeed, there is a strong tendency to reciprocate ties.



## Goodness of fit 

### Model comparison

ERG models do not have associated goodness of fit  measures like $R^2$ in
linear regression.  To compare different models a Akaike Information Criterion
(AIC) measure [@akaike_1998] is used. Models with smaller AIC should be
preferred.  AIC can be used to compare non-nested models. Its design penalizes
models with more parameters, so adding more terms to the ERG model does not
necessarily lead to better fitting models according to AIC. 

AIC can be calculated with function `AIC`.  For the two models estimated in
previous sections we obtain:

```{r, aic}
AIC(model.ibe0, model.ibe1)
```

So the "reciprocity" model seems to fit the data better than homogeneous
Bernoulli model.



### Goodness of fit through simulation

To examine how well a given ERG model fits the data a simulation-based methods
are used.  Given an estimated model we can simulate networks consistent with
it. Once we simulate a lot of networks from the model, we can check, whether
these simulated networks posses the same kind of global network properties as
the observed network.

Global properties usually analyzed in this way include:

- (in-/out-)degree distributions.
- distribution of shortest path lengths (minimum geodesic distances).
- distribution of the number of edge-wise shared partners.

This can be done with the function `gof`. Let us examine the homogeneous
Bernoulli model fitted earlier.

```{r, gof_bernoulli}
fit0 <- gof(model.ibe0)
layout(matrix(1:4, 2, 2, byrow=TRUE))
plot(fit0)
```

We are presented with four charts. Each chart presents two types of data:

First, solid thick black lines represent properties (distributions) of our
observed data (the `ibe` network). We have in-degree distribution, out-degree
distribution, distribution of the number of edge-wise shared partners, and the
distribution of shortest path lengths.

Second, boxplots summarize the properties of the model-based simulated
networks. By default 100 networks are simulated. Due to the random nature of
the simulation procedure each simulated network will be a little bit different
from the others. That's why boxplots are used to visualize how the simulated
distributions vary.

To examine the fit, we have to compare visually empirical data (solid thick
black lines) to simulated distributions (boxplots). Median values of simulated
distributions are shown with black horizontal lines within the boxes. If our
model would fit the data perfectly, the empirical lines should follow the
medians of all the boxes. If, on the other hand, the empirical line strays far
away, even beyond the gray lines symbolizing 95% confidence intervals, it is an
indication of a very poor fit.

From the charts above we see that the Bernoulli model fits the data very poorly
indeed. While the out-degree distribution the model is able to recover, the
remaining three statistics are not modeled in a satisfactory manner.


Let us examine the goodness of fit of the "reciprocity" model.

```{r, gof_reciprocity}
fit1 <- gof(model.ibe1)
layout(matrix(1:4, 2, 2, byrow=TRUE))
plot(fit1)
```

Again, we see that this model also does not fit the data very well.  This means
that there must be other social process, apart from reciprocity, that govern
the structure of the network in this classroom.  We will explore other
possibilities in the sections below.











# (In)dependence assumptions        {#independence_assumptions}

You are probably familiar with Ordinary Least Squares (OLS) regression:

$$
Y_i = \alpha + \beta X_i + \varepsilon_i
$$

If you recall, on of the assumption of OLS is that observations are
independent, which translates to the requirement
$\operatorname{Cor}{(\varepsilon_i, \varepsilon_j)} = 0$.  In other words, that
the correlation between errors for any two observations is 0. Some
generalizations of OLS regression go in the direction of relaxing this
assumption, by allowing some observations to be correlated. In particular,
multilevel models [e.g. @snijders_bosker_2011] can be considered such a
generalization as, while we assume zero correlation between observations
belonging to different groups, observations belonging to the same group can be
correlated.

Part of the challenge in modeling the structure of social networks is in
dependence.  Social networks represent relationships between actors, which make
these actors and their actions dependent in some respect. In particular, the
process through which they form relations with others may be dependent.
Consequently, in social networks, in principle, every tie may depend on other
ties in the network. To approach this "Hell of Dependence" in a tractable way,
proposed models can be characterized with *independence assumptions* that they
pose. These statements specify what kind of independence, and between which
ties, is assumed behind a given model. Much like in OLS regression, a lot of
new developments in statistical analysis of social networks have been made by
postulating new assumptions that relax some independence assumptions
made by previous ones.

In this section we review the following independence assumptions.

1. The Bernoulli assumption.
2. Dyadic independence assumption.
3. Markov dependence assumption.
4. Other forms of dependence.

Practical examples showing ERG models that conform to particular assumptions are
presented in the [section with practical examples](#practical_examples).



## Bernoulli assumption        {#ass_bernoulli}

The Bernoulli assumption is the simplest and the strongest. We have already
discussed it in the section with the [homogeneous Bernoulli model example](#model_bernoulli) above. 
According to this model every *tie* is assumed to be a realization of an
independent [Bernoulli Trial](http://en.wikipedia.org/wiki/Bernoulli_trial), a
flip of a coin. The probability that an arc from $i$ to $j$ exist does not
depend in any way on the structure of the remaining ties of the network. In
particular, it does not depend on how many ties $i$ and $j$ already have,
whether they have any network neighbors in common, and so on. It also does not
depend on whether or not there is an arc from $j$ to $i$.



## Dyadic independence assumption      {#ass_dyadic}

Recall that a dyad is a pair $(Y_{ij}, Y_{ji})$ that is the state of ties
between actors $i$ and $j$. According to the Bernoulli assumption all
$Y_{ij}$, for all $i$ and $j$, are independent from one another. In particular,
$Y_{ij}$ is assumed independent from $Y_{ji}$, as mentioned in the previous
section.

The assumption of *dyadic independence* relaxes the assumption of independence
*within dyads*. Consequently, $Y_{ij}$ and $Y_{ji}$ can be dependent.

The reciprocity model [discussed previously](#reciprocity) is a
dyad-independent model.  The presence or absence of a tie from $i$ to $j$
affects the probability of tie from $j$ to $i$ because of the presence of the
`mutual` term.  However, it does not affect the probabilities of any other ties
in the network.

There are also other important dyad-independent models:

- Models that include node covariates. Sending or receiving ties might be
  affected by a value of a nodal attribute, e.g., gender, age, and so on.
- Models that include dyadic covariates. Ties may be more likely between nodes
  that are similar to one another (e.g. same gender, smaller age difference) or
  are in some way closer to one another, e.g., are separated by a smaller
  geographical distance.

We will explore both types of models when analyzing homophily in the subsequent
sections.

Some dyad-independent ERGMs simplify to Generalized Linear Models, which we will
show in section [Dyad-independent ERGMs](#dyadic_independence).




## Markov dependence              {#ass_markov}

Assumption of *Markov dependence* was proposed by @frank_strauss_1986 and
further relaxes the dyadic independence assumption.

Dyadic independence implies that all dyads are independent, including those
involving the same actor, i.e.: 

- $Y_{ij}$ is independent from $Y_{ik}$ (the same tie sender)
- $Y_{ij}$ is independent from $Y_{kj}$ (the same tie receiver)

Markov dependence assumes that dyads are independent *unless* share an actor on one of the endpoints. See left on the figure below. The probability of tie $ij$ might not be independent from the probability of tie $jk$ because these are both ties of actor $j$.

```{r illustration_dependence, echo=FALSE, out.width="0.4\\textwidth", small_mar=TRUE}
layout(matrix(1:2, 1, 2))
g <- igraph::graph.formula(i -- j, j -- k)
l <- matrix(c(1,2, 2,1, 3,2), 3, 2, byrow=TRUE)
igraph::plot.igraph(g, layout=l, vertex.color="white", vertex.size=40,
                    edge.width=2, vertex.label.color="black", vertex.label.font=3)
title(main="Markov dependence")

g <- igraph::graph.formula(i -- j, j -- u, u -- v, v -- i)
l <- matrix(c(1,1, 1,2, 2,2, 2,1), 4, 2, byrow=TRUE)
igraph::plot.igraph(g, layout=l, vertex.color="white", edge.lty=c(1,2,2,1),
                    vertex.size=40, vertex.label.color="black",
                    edge.width=2, vertex.label.font=3)
title(main="Partial conditional dependence")
```

ERG models conforming to the Markov dependence assumption allow to
model the shape of degree distribution and processes like preferential
attachment or transitivity (e.g. meeting friends through friends, etc.).

Further relaxations of independence assumptions lead to the so-called "Socal Circuit Model" in which two ties $ij$ and $vu$ are assumed independent *unless* there is a link $iu$ or $jv$.
See @koskinen_daraganova_2012 for further details.








# ERG models and logistic regression        {#ergm_and_logit}

Dyadic independence is a very strong assumption which rarely (if ever) holds in
real social networks. Dyadic independence ERGMs treat every dyad as a single
object, independent from its surrounding.  Therefore, in such models we could
use only such measures, that are dyadic independent. That means that we could
calculate change statistic knowing only the state of a given dyad and maybe
some attributes of nodes (in this dyad). For instance number of edges is dyadic
independent -- if we toggle an edge from 0 to 1 we are sure that number of
edges will increase by 1, no matter how the rest of the network looks like.
Another example is the homophily effect -- we could compute change in the
number of edges between nodes of the same type knowing only attributes of these
two specific nodes in a dyad. On the other hand, number of 2-stars is dyadic
dependent - we need two know structure of the neighborhood if we want to
calculate the change in the number of twopaths.

Assuming that our statistics are dyadic independent we could rewrite edge
probability to logit form:

$$
\operatorname{logit}{(P_\theta(A_{ij}=1))}  = \theta^T \delta[g(a)]_{ij}
$$

That looks like a logistic regression model. Indeed, dyad-independent ERGMs
simplify to a logistic regression for data in which dyads are observations, the
dependent variable is binary and equals to 1 if an edge exists and 0 otherwise.
Independent variables are the change statistics (one for each term).






# Practical examples of ERG models       {#practical_examples}

In the following subsections we elaborate some example ERG models fit to the
classroom data introduced above. The primary goal of these examples is to
illustrate

- how to specify different ERG models matching different research questions
- show ERG models pretending to different independence assumptions (see [section on (in)dependence assumptions](#independence_assumptions)).

Below we present brief examples of modeling homophily and triadic closure with ERGMs.



## Actor relation effects: gender homophily in a school class

The first example shows how ERG models can be used to model homophily.  Topics
of homophily and segregation in social networks are covered in a separate
tutorial.

Recall that homophily is a tendency for ties to form between similar actors.
Similarity need to be specified according to a specific node attribute, which
can be continuous or categorical. Let us focus on categorical attribute first.

The key ERGM term in modeling homophily on a categorical attribute is `nodematch`.







Consider the following example of fitting a dyad-independent ERGM to classroom
data.

```{r}
model.ibe1 <- ergm(ibe ~ edges + nodematch("female"))
summary(model.ibe1)
```

Let us now try to fit the same model using logistic regression. We need to
prepare the data first.

```{r, mm}
mm <- isnar::mixingm(playnet, "female", full=TRUE)
mm
```

Object `mm` is a so-called *mixing matrix* that cross-classifies all dyads in
the network according to three characteristics:

1. Gender of first node.
2. Gender of second node.
3. Whether there is a tie from first node to second node or not.

For example, there are `r mm[2,2,1]` girl-girl pairs that are connected with a
tie and only `r mm[1,2,2]` ties that are send by a boy towards a girl, and so
on.

We now transform the mixing matrix to a data frame:

```{r, mm_df}
d <- as.data.frame(as.table(mm))
d
```

Columns `ego` and `alter` mark whether, respectively, tie sender or tie
receiver is a female.  Variable `tie` is `TRUE` for connected pairs. Finally,
variable `Freq` contains the frequency.

We now add a variable `match` differentiating same-gender pairs: `match` is
`TRUE` whenever `ego` and `alter` are equal:

```{r, add_match}
d$match <- with(d, ego == alter)
d
```

We are now ready to fit our logit model. Our binary dependent variable is `tie`
(is there a tie or not). Our only independent variable is `match`, also binary,
representing whether dyad involves actors of the same sex. We weight the cases
with variable `Freq`.

```{r, logit}
ibe.logit <- glm( tie ~ match, data=d, weight=Freq, family=binomial("logit"))
ibe.logit
```


As you can see, the coefficients are identical:

```{r}
cbind( ERGM=coef(model.ibe1), Logit=coef(ibe.logit))
```

In the next section we show how homophily can be modeled in some more detail.






## Gender homophily example

In our attempt to model the `ibe` network let's pursue the topic of gender
*homophily*. Network homophily is a pattern in which ties are more likely to
exist between nodes similar to each other according to some attribute, e.g.
gender [@mcpherson_etal_2001; bojanowski_corten_2014].

ERG models for homophily are also dyad-independent models, similar to the
one from the previous section. Conditional probability of a tie exist between
two nodes can be represented as a function of attributes of the first node,
attributes if the second node, and possible interaction effects.

```{r, homophily_main}
model.ibe2a <- ergm(ibe ~ edges + mutual + nodeofactor("female") +
                    nodeifactor("female"),
                  control=control.ergm(seed=666))
summary(model.ibe2a)
```




```{r, homophily_mixing}
model.ibe2b <- ergm(ibe ~ edges + mutual +  nodeofactor("female") + nodeifactor("female")
                    + nodematch("female"),
                  control=control.ergm(seed=666))
summary(model.ibe2b)
```

Compare Bernoulli model, reciprocity model, and the two homophily models

```{r, homophily_AIC}
AIC(model.ibe0, model.ibe1, model.ibe2a, model.ibe2b)
```

Homophily model seems to fit best. Let's examine the (lack of) fit in some more
detail with the `gof` function.

```{r gender_homophily_gof}
fit2b <- gof(model.ibe2b)
layout(matrix(1:4, 2, 2, byrow=TRUE))
plot(fit2b)
```















## Social circuit models

Here we present an example of an ERGM which contains a statistic which subsumes to the "Social Circuit" dependence assumptions presented earlier. It is Geometrically Weighted Edgwise Shared Partners (GWESP). With this statistic we model the distribution of the number of network partners of nodes, that are themselves connected (hence "edge" in the name of the statistic). Positive effect indicates the tendency for creating transitive triplets.

Let us fit an ERGM with the GWESP statistic together with gender homophily effects investigated earlier.

```{r gwesp2}
model.ibe4 <- ergm(ibe ~ edges + mutual + 
                   nodeofactor("female") + nodeifactor("female") + nodematch("female") +
                   gwesp(alpha=0.2, fixed=TRUE), 
                 control=control.ergm(seed=666))
fit4 <- gof(model.ibe4)
layout(matrix(1:4, 2, 2, byrow=TRUE))
plot(fit4)
```

Model fit still not ideal. In particular, the models does not explain the spike in the in-degree distribution for degree 2. Nevertheless, the fit is much better than in models estimated earlier. Let us have a look at the estimates:

```{r gwesp2_model_summary}
summary(model.ibe4)
```

From the results we can observe:

- Friendship nominations are reciprocated
- There is a strong gender homophily effect
- GWESP effect is positive, hence friendship nominations are transitive. In other words, ceteris paribus, friend of a friend is likely to be a friend too.















# More technical details for interested and technically inclined readers

In the following two sections we present a somewhat more detailed discussion of
how an ERG model can be interpreted as:

1. A model for a probability distribution over all networks of a given size.
2. A model for conditional probabilities of tie existence given the structure of the remainder of the network.

Interpretation (1) can be called *global* while interpretation (2) can be
called *local*.





## Model-based network probabilities

This is a more advanced topic.

As mentioned, ERGMs assign probabilities to all networks with given size
(number of nodes) according to some statistic (measure/attribute). To see how
it works let's consider undirected network of size 4. There are 64 such
networks in total, but some of them differ only in node permutation, so there
are only 11 topologically different networks.

Assume we have three models: null model, model with one parameter (number of
edges) and with two parameters (number of edges and 2-stars/twopaths).
Parameters are equal $\theta_1=-0.5$ for number of edges (both models) and
$\theta_2 = 0.2$ for number of twopaths. We will calculate probability of each
type of network under all three models.

First we need to generate all possible networks of size 4 and select one
representative for each canonical form.

```{r, fig.width=5, fig.height=5, fig.align='center'}
# all 4-node networks
adj <- as.matrix(expand.grid(0:1, 0:1, 0:1, 0:1, 0:1, 0:1))
full_edgelist <- subset(expand.grid(1:4, 1:4), Var1 > Var2)
nets <- lapply(seq(nrow(adj)), function(i) {
  network.edgelist(full_edgelist[as.logical(adj[i,]), ],
                   network.initialize(4, directed = FALSE))
})
nets <- nets[order(sapply(nets, network.edgecount))]

# unique canonical permutation (based on degree sequence, works for 4 nodes)
degrees <- as.data.frame(t(sapply(nets, function(net) summary(net ~ degree(0:3)))))
degrees <- data.frame(degrees, id = apply(degrees, 1, paste, collapse = ""))
canonical_count <- table(degrees$id)
nets_unique <- nets[match(unique(degrees$id), degrees$id)]
```

Now we calculate the numerator in equation for probability for every canonical
form and each model. Under null model every network is equally probable, so we
don't have to compute anything. For other functions small function will come in
handy.

```{r}
model0 <- rep(1, 11)

prob <- function(net, coeff) {
  form <- as.formula(paste("net ~", paste(names(coeff), collapse = "+"), collapse =" "))
  z <- summary(form)
  exp(sum(z * coeff))
}
model1 <- sapply(nets_unique, prob, coeff = c(edges = -0.5))
model2 <- sapply(nets_unique, prob, coeff = c(edges = -0.5, twopath = 0.2))
```

Next step is calculating proper probabilities. To do so we sum numerators over
all canonical forms taking the size of each group into account.

```{r}
model0 <- model0 / sum(model0 * canonical_count)
model1 <- model1 / sum(model1 * canonical_count)
model2 <- model2 / sum(model2 * canonical_count)
```

And plot them.

```{r, echo=FALSE}
library("RColorBrewer")
pal <- brewer.pal(3, "Set2")

par(mfrow = c(6,4), mar = c(2, 1, 1, 1))
for (i in c(1,7,2,8,3,9,4,10,5,11,6)) {
  plot(nets_unique[[i]], coord = matrix(c(0,0,1,1,0,1,0,1), ncol = 2), pad = 0, 
       jitter = FALSE)
  barplot(rev(c(model0[i], model1[i], model2[i])), horiz = TRUE, col = rev(pal),
          xlim = c(0, 0.06))
}
plot.new()
plot.new()
legend("center", legend = c("Null model", "Edges only", "Edges and twopaths"),
       fill = pal)
par(mfrow = c(1, 1))
```

All networks has the same probabilities under null model obviously. It is also
clear that under model with edges only probability is decreasing as density (or
edge count) is increasing. This is caused by negative edge parameter - less
dense networks are preferable. Last model is the most interesting one - we
could notice that probability decreases in the beginning as density increases,
but afterwards it starts to increase again. Edge parameter is still negative so
dense networks are penalised. However, second parameter is positive and number
of twopath (2-stars) increases more or less faster than number of edges.
Therefore probabilities of dense networks are increasing. 




## Conditional edge probabilities

This a more advanced topic.

The model could be equivalently described by conditional probabilities that an
edge exists given the rest of the network. Conditional log-odds could be easily
derived from the general form of the model

$$\frac{ P_\theta ( A_{ij}=1 | A_{ij}^c = a_{ij}^c ) }{ P_\theta ( A_{ij}=0 | A_{ij}^c = a_{ij}^c ) }=
\exp\{\theta^T \delta[g(a)]_{ij}\}$$

where $\delta[g(a)]_{ij}$ is the change of $g(a)$ when $a_{ij}$ is changed from 0 to 1. So the probability of an edge is equal to

$$
\begin{split} 
 P_\theta(A_{ij}=1|A_{ij}^c=a_{ij}^c) &=  \frac{\exp\{\theta^T \delta[g(a)]_{ij}\}}{1+\exp\{\theta^T \delta[g(a)]_{ij}\}} \nonumber\\
 & =  (1 + \exp\{-\theta^T \delta[g(a)]_{ij}\})^{-1} ,
\end{split}
$$

Coming back to 4-nodes networks. Let's say we have a model with two statistics
defined above: number of edges and number of 2-stars (twopaths) with
corresponding parameters $\theta_1=-0.5$ and $\theta_2 = 0.2$.

```{r, echo=FALSE, fig.width=8, fig.height=5, fig.align='center'}
source("function_compute_prob.R")
coeff <- c(edges = -0.5, twopath = 0.2)
probs <- lapply(nets_unique, function(net) {
  p <- compute_prob(net, coeff)
  p[upper.tri(p)]
})

par(mfrow = c(3, 4), mar = rep(1,4))
for (i in seq_along(nets_unique)) {
  plot(nets_unique[[i]], coord = matrix(c(0,0,1,1,0,1,0,1), ncol = 2), pad=0.3,
       jitter = FALSE, main=paste0("(", i, ")"), vertex.cex=3, edge.col="grey")
  text(x = c(-0.2, 0.5, 0.3, 0.3, 0.5, 1.2), 
       y = c(0.5, -0.1, 0.15, 0.85, 1.1, 0.5),
       zapsmall(probs[[i]], digits = 2))
}
par(mfrow = c(1, 1))
```


```{r, results="hide", echo=FALSE}
p <- exp(-0.5+0.6)/(1+exp(-0.5+0.6))
```

Look close on 7th network. Where has probability `r round(p, 2)` come from? So assume
that we create an edge between top two nodes. Then number of edges is obviously
increased by 1, but number of twopaths is increased by 3 -- new edge creates
twopath with every other edges. Now we have change statictics so we could
compute probability from equation

$$
 \frac{\exp\{-0.5 \cdot 1 + 0.2 \cdot 3\}}{1+\exp\{-0.5 \cdot 1 + 0.2 \cdot 3\}} = `r exp(-0.5+0.6)/(1+exp(-0.5+0.6))`
$$






```{r warns, echo=FALSE, results="hide"}
stopifnot(is.null(warnings()))
```



# References

<!--
vim: wrap:linebreak:spell:spelllang=en_us
-->