---
title: "Diversity and Regression"
author: "by James Hollway"
output: 
  learnr::tutorial:
    theme: journal
runtime: shiny_prerendered
description: >
  This tutorial aims to teach you how to measure and test network diversity,
  moving from univariate to multivariate tests, including 
  network linear models (multiple regression quadratic assignment procedures).
---

```{r pkgs, include = FALSE}
library(manynet)
library(autograph)
library(migraph)
```

```{r setup, include = FALSE, purl=FALSE, eval=TRUE}
library(manynet)
library(autograph)
library(migraph)
library(learnr)
knitr::opts_chunk$set(echo = FALSE)
manynet::clear_glossary()
learnr::random_phrases_add(language = "fr",
                           praise = c("C'est génial!",
                                      "Beau travail",
                                      "Excellent travail!",
                                      "Bravo!",
                                      "Super!",
                                      "Bien fait",
                                      "Bien joué",
                                      "Tu l'as fait!",
                                      "Je savais que tu pouvais le faire.",
                                      "Ça a l'air facile!",
                                      "C'était un travail de première classe.",
                                      "C'est ce que j'appelle un bon travail!"),
                          encouragement = c("Bon effort",
                                            "Vous l'avez presque maîtrisé!",
                                            "Ça avance bien.",
                                            "Continuez comme ça.",
                                            "Continuez à travailler dur!",
                                            "Vous apprenez vite!",
                                            "Vous faites un excellent travail aujourd'hui."))
learnr::random_phrases_add(language = "en", 
                           praise = c("C'est génial!",
                                      "Beau travail!",
                                      "Bravo!",
                                      "Super!"),
                           encouragement = c("Bon effort"))
marvel_friends <- to_unsigned(to_uniplex(fict_marvel, "relationship"), keep = "positive")
marvel_friends <- to_giant(marvel_friends)
marvel_friends <- marvel_friends %>% to_subgraph(Appearances >= mean(Appearances))
```

## This tutorial

<img src="https://i.giphy.com/media/v1.Y2lkPTc5MGI3NjExeDFidmtrcTJhenU1NGRtdXA5N292Ynp3Y3BxcHRhYjR2MnowOGJ0ZiZlcD12MV9pbnRlcm5hbF9naWZfYnlfaWQmY3Q9Zw/djRJNZqj508sE/giphy.gif" alt="gif of deadpool greeting" width = "900"/>

Within a network, nodes and ties vary.
But how much do they vary and are there associations in this variance?
By the end of this tutorial, you will be able to:

- [ ] &nbsp; measure the richness of a network in terms of attribute types
- [ ] &nbsp; measure how diverse a network is
- [ ] &nbsp; measure how heterophilous the ties in a network are
- [ ] &nbsp; test how heterophilous the ties in a network are compared to random networks
- [ ] &nbsp; test how heterophilous the ties in a network are compared to permutations of the original network
- [ ] &nbsp; regress a network's ties on other ties and/or nodal or dyadic attributes

## Initial visualisation

For this session, we'll explore a couple of different datasets.
First, let's examine homogeneity/`r gloss("heterogeneity")` in the Marvel dataset from `{manynet}`,
`fict_marvel`.
The dataset is quite complicated,
so to make this simpler, let's concentrate on:

- just the relationships part of the network
- just the positive (friendship) ties and not the negative (enmity) ties
- the main (giant) component without any isolates
- just those characters that appear in the comics more than average

Fortunately, all these data cleaning moves are easy to do in `{manynet}`,
and can be seen in the following chunk in order:

```{r friends, exercise=TRUE, purl = FALSE}

```

```{r friends-hint-1, purl = FALSE}
to_uniplex(____, "relationship")
```

```{r friends-hint-2, purl = FALSE}
# since the dataset is a 'signed' graph, we want to get just the
# positively signed ties to get the friendship graph 
# (and lose the enmity relations)
to_unsigned(____, keep = "positive")
```

```{r friends-hint-3, purl = FALSE}
# to_giant() is a quick easy way to get the giant/main component
to_giant(____)
```

```{r friends-hint-4, purl = FALSE}
to_subgraph(____, Appearances >= mean(Appearances))
```

```{r friends-hint-5, purl = FALSE}
marvel_friends <- to_unsigned(ison_marvel_relationships, keep = "positive")
marvel_friends <- to_giant(marvel_friends)
marvel_friends <- marvel_friends %>% to_subgraph(Appearances >= mean(Appearances))
marvel_friends
```

```{r friends-solution}
marvel_friends <- to_unsigned(ison_marvel_relationships, keep = "positive")
marvel_friends <- to_giant(marvel_friends)
marvel_friends <- marvel_friends %>% to_subgraph(Appearances >= mean(Appearances))
marvel_friends
```

This gives us an `r gloss("undirected")` network of nearly twenty characters and a little more than 100 edges.
Recall that this data has several nodal attributes.
Let's explore a couple of these attributes, "Gender" and "PowerOrigin", visually
using `graphr()`.

```{r plotfriends, exercise=TRUE, purl = FALSE, fig.width=9}

```

```{r plotfriends-hint, purl = FALSE}
# Since both Gender and PowerOrigin are categorical variables
# you will need to use two different aesthetic dimensions to
# represent them together.
# Which will you present as shapes and which as colors?
graphr(____, 
       node_shape = ____,
       node_color = ____)
```

```{r plotfriends-solution}
graphr(marvel_friends, 
       node_shape = "Gender",
       node_color = "PowerOrigin")
```

These variables seem to be distributed unevenly across the network.
There seems to be some `r gloss("homophily")` -- or like choosing like --
operating here, but it is difficult to say conclusively because there are
clearly more male than female superheros,
as well as clearly more superheros of mutant origin than others.
So what might seem like homophily could just be a result of there being
many more opportunities for ties between nodes of some categories.
We therefore need to establish how diverse this network really is.

## Measuring richness

Some categorical variables have more categories than others,
and this can be thought of as indicating greater diversity than e.g. a binary
categorical variable.
A measure of this 'richness' is inherited from the study of biodiversity.
`r gloss("Richness","richness")` calculates the number of different categories that appear in
the data for a given variable.
Let's calculate the richness of all of the variables in our Marvel data.

```{r rich, exercise=TRUE, purl = FALSE}
net_node_attributes(marvel_friends)
```

```{r rich-hint, purl = FALSE}
net_richness(____, ____)
```

```{r rich-solution}
net_richness(marvel_friends, "Gender")
net_richness(marvel_friends, "PowerOrigin")
net_richness(marvel_friends, "Attractive")
net_richness(marvel_friends, "Rich")
net_richness(marvel_friends, "Intellect")
```

```{r richness-question, echo=FALSE, purl = FALSE}
question("Which variable is the most 'diverse' according to this richness measure?",
  answer("Gender"),
  answer("PowerOrigin", 
         correct = TRUE,
         message = "There are four categories available in this data for power origin, while the other variables include only two categories each."),
  answer("Attractive"),
  answer("Rich"),
  answer("Intellect"),
  random_answer_order = TRUE,
  allow_retry = TRUE
)
```

Note though that 'richness' as a network measure does not include 
any sense of how distributed these categories are around the network.
There is a measure of nodal richness available,
which counts the number of different categories to which each node is tied,
but this does not offer a summary of how evenly distributed the appearance
of categories are.
For that we would need to move on to something like the Blau index.

## Measuring diversity

Another measure that reflects the diversity or `r gloss("heterogeneity")` in the 
network for each attribute is the Blau Index.
Recall that the Blau index for any given variable is:

$$1 - \sum p_i^2$$

where $p_i$ represents the proportion of observations holding $i$, 
which indexes each of the given categories.
A perfectly homogeneous group (all the same category) would receive a score of 0,
while a perfectly `r gloss("heterogeneous","heterogeneity")` group (with members spread evenly over the maximum categories)
would receive a score of 1.
Obtain the network diversity scores for the five attributes used above.

```{r blau, exercise=TRUE, purl = FALSE}

```

```{r blau-hint, purl = FALSE}
net_diversity(____, ____)
```

```{r blau-solution}
net_diversity(marvel_friends, "Gender")
net_diversity(marvel_friends, "PowerOrigin")
net_diversity(marvel_friends, "Attractive")
net_diversity(marvel_friends, "Rich")
net_diversity(marvel_friends, "Intellect")
```

Looks like there is more diversity in terms of where these characters got
their powers, whether they have significant intellectual powers,
and their gender, than their attractiveness or their wealth.

We can also cross-reference this diversity.
For example, we might be interested in whether our comic book heroes
are equally gender diverse across their (power) origin stories,
or equally intellectually diverse across gender.^[Note that this works for calculated categorical variables too, such as cluster/group assignment from community detection or equivalence classes.]

```{r crossref, exercise=TRUE, purl = FALSE}

```

```{r crossref-hint, purl = FALSE}
over_membership(.data, FUN, ..., membership)
```

```{r crossref-solution}
over_membership(marvel_friends, net_diversity, "Gender", 
                membership = "PowerOrigin")
over_membership(marvel_friends, net_diversity, "Intellect", 
                membership = "Gender")
```

Note that the length of the vector returned as a result
is the number of categories in the second category listed.
It looks like some origin stories are more gender diverse than others.
Gods (just Thor here) and humans are all men,
whereas those with mutant or radiation origin stories are more gender diverse.
There doesn't appear to be much difference in intellect 
across gender categories in this data.

Ok, this tells us about how (un)even the distribution of these variables is in this network,
but it doesn't necessarily tell us whether ties are appearing more frequently
_between_ nodes of similar (or different) categories.
For that we need to look at `r gloss("homophily")`/`r gloss("heterophily.", "heterophily")`

## Measuring heterophily

The EI (or E-I) index offers a way to measure the degree to which ties
appear between rather than within groups of nodes of the same category.
Recall that the EI index is calculated as:

$$\frac{E-I}{E+I}$$

where $E$ is the number of ties present between a variable's categories (i.e. external),
and $I$ is the number of ties present within a variable's categories (i.e. internal).
As such, an EI index of -1 suggests perfect homophily, whereas an EI index of +1 suggests perfect `r gloss("heterophily.","heterophily")`
(This is why the function is called `net_heterophily()`).

Check how homophilic three variables in the network are,
"Gender", "PowerOrigin", and "Attractive".
Please assign the results so that we can use them later.

```{r ei, exercise=TRUE, purl = FALSE}

```

```{r ei-hint, purl = FALSE}
net_heterophily(____, ____)
```

```{r ei-solution}
(obs.gender <- net_heterophily(marvel_friends, "Gender"))
(obs.powers <- net_heterophily(marvel_friends, "PowerOrigin")) 
(obs.attract <- net_heterophily(marvel_friends, "Attractive")) 
```

```{r homophily-present, echo=FALSE, purl = FALSE}
question("For which variables is there a sign of some homophily effect according to the EI index? (Choose all that apply)",
  answer("Gender", 
         correct = TRUE, 
         message = "Yes, looks like there might be some gender homophily present."),
  answer("PowerOrigin", 
         message = "The score for power origin homophily is so close to 0 that it does not seem to signal much."),
  answer("Attractive", 
         correct = TRUE, 
         message = "And looks like a fairly large effect for homophily on the basis of looks..."),
  allow_retry = TRUE
)
```

Ultimately though, these are just scores,
and doesn't tell us whether this is any more or less than
what we might expect the score to be by chance for a network
of this size and density and distribution of that attribute.

## CUG tests

To see whether we should be surprised by scores this high/low, 
we compare these scores with those from a series of `r gloss("random graphs","random")`
(Erdös-Renyi/Bernoulli) of the same dimensions and
distribution of the attribute.
This can help us establish whether there is
more homophily or heterophily than expected by chance.

This is often called a `r gloss("conditional uniform graph","cug")` or CUG test,
but `{migraph}` uses more descriptive function names,
such as `test_random()`.
Plot the results of running this function with respect to the EI index on each of the three variables.
These plots can be made if you have `{autograph}` installed and loaded.
You can specify that one thousand simulations should be used using `times = 1000`.

```{r rando, exercise=TRUE, purl = FALSE, exercise.timelimit = 120, fig.width=9}

```

```{r rando-hint-1, purl = FALSE}
rand.____ <- test_random(____, FUN = ____, attribute = ____, times = ___)
```

```{r rando-hint-2, purl = FALSE}
plot(rand.____)
```

```{r rando-hint-3, purl = FALSE}
rand.gender <- test_random(marvel_friends, 
                           net_heterophily, attribute = "Gender", 
                           times = 1000)
rand.power <- test_random(marvel_friends, 
                          net_heterophily, attribute = "PowerOrigin", 
                          times = 1000)
rand.attract <- test_random(marvel_friends, 
                            net_heterophily, attribute = "Attractive", 
                            times = 1000)
plot(rand.gender) + ggtitle("CUG test results for 'Gender' attribute")
plot(rand.power) + ggtitle("CUG test results for 'PowerOrigin' attribute")
plot(rand.attract) + ggtitle("CUG test results for 'Attractive' attribute")
```

```{r rando-solution}
rand.gender <- test_random(marvel_friends, 
                           net_heterophily, attribute = "Gender", 
                           times = 1000)
rand.power <- test_random(marvel_friends, 
                          net_heterophily, attribute = "PowerOrigin", 
                          times = 1000)
rand.attract <- test_random(marvel_friends, 
                            net_heterophily, attribute = "Attractive", 
                            times = 1000)
plot(rand.gender) + ggtitle("CUG test results for 'Gender' attribute")
plot(rand.power) + ggtitle("CUG test results for 'PowerOrigin' attribute")
plot(rand.attract) + ggtitle("CUG test results for 'Attractive' attribute")
```

The plots of these results use a dotted vertical line for 0 where this is in bounds,
a red vertical line for the observed score,
and a density plot of the scores from the randomly generated networks.
The grey tails of the distribution are a visual aid indicating the most extreme 5% of scores
from the distribution.
The results are _really_ interesting.

<img src="https://i.giphy.com/media/v1.Y2lkPTc5MGI3NjExNzlsM3J1dm9tb2tlYjB4bGw0cGxuc3RodXpvdjNhcGJ5cXV2MndjbiZlcD12MV9pbnRlcm5hbF9naWZfYnlfaWQmY3Q9Zw/AbYxDs20DECQw/giphy.gif" alt="gif of ironman explosion" width = "900"/>

Despite being the larger coefficients (in absolute terms),
it looks like we cannot reject the null hypothesis that there is no homophily
for gender nor for attractiveness.
Both observed scores fall within the range of scores we would expect from
randomly generated networks with the same distribution of that variable.

However, we can reject the null hypothesis with respect to their power origin story.
While the coefficient itself is close to 0 (neither strong homophily nor heterophily),
all the random networks generated returned larger EI scores, between .1 and .4.

```{r cuginterp-qn, echo=FALSE, purl = FALSE}
question("Which of the following is true?",
  answer("Attractive superheros only hang out with other attractive superheros.",
         message = paste(learnr::random_encouragement(), '<img src="https://i.giphy.com/media/v1.Y2lkPTc5MGI3NjExeW9jbXg5dDBndjRpcWc2enhyOTZvbjloeGFpbmJndXNrMm05MHdtZiZlcD12MV9pbnRlcm5hbF9naWZfYnlfaWQmY3Q9Zw/ziLadIVnOGCKk/giphy.gif" alt="gif of thor frustration"/>')),
  answer("There is significantly less power origin heterophily than expected.",
         correct = TRUE,
         message = learnr::random_praise()),
  answer("There is no evidence of gender homophily.",
         correct = TRUE,
         message = learnr::random_praise()),
  answer("There is significant power origin homophily.", 
         message = "Remember that the coefficient does not suggest many homophilous ties."),
  random_answer_order = TRUE,
  allow_retry = TRUE
)
```

## QAP tests

Ah, but perhaps random graphs are not the best reference group 
for establishing whether there is a significant homophily effect here.
After all, social networks are not completely random; 
they are _structured_ in particular ways, such as some nodes having higher degrees
than others, or there being a `r gloss("core and periphery","core")` or `r gloss("community")` topology.

Another approach to establishing a baseline for whether we should be surprised
by a given score or not is to use `r gloss("permutations","permutation")` of the underlying network
instead of random graphs.
Permuting the network retains the structure of the network
because the ties are kept and only the labels (variables) are reassigned randomly.
Let's first plot the observed data and some permuted data next to each other.

```{r perm, exercise=TRUE, purl = FALSE, fig.width=9}

```

```{r perm-hint, purl = FALSE}
graphr(generate_permutation(____, with_attr = TRUE), ____)
```

```{r perm-solution}
old <- graphr(marvel_friends,
              labels = FALSE, node_size = 6,
              node_color = "PowerOrigin",
              node_shape = "Gender") + ggtitle("Original network")
new <- graphr(generate_permutation(marvel_friends, with_attr = TRUE),
              labels = FALSE, node_size = 6,
              node_color = "PowerOrigin",
              node_shape = "Gender") + ggtitle("Permuted network")
old + new
```

```{r cupqap-qn, echo=FALSE, purl = FALSE}
question("Which of the following is true?",
  answer("Random networks retain the structure of the original network.",
         message = paste(learnr::random_encouragement(), '<img src="https://i.giphy.com/media/v1.Y2lkPTc5MGI3NjExeW9jbXg5dDBndjRpcWc2enhyOTZvbjloeGFpbmJndXNrMm05MHdtZiZlcD12MV9pbnRlcm5hbF9naWZfYnlfaWQmY3Q9Zw/ziLadIVnOGCKk/giphy.gif" alt="gif of thor frustration"/>')),
  answer("Permuted networks retain the structure of the original network.",
         correct = TRUE,
         message = learnr::random_praise()),
  answer("Both random and permuted networks retain the proportion of attributes from the original network.",
         correct = TRUE,
         message = learnr::random_praise()),
  answer("Permuted networks retain the ties among the same nodes from the original network.", 
         message = "Permuted networks randomly reassign any variables and change positions of nodes by swapping the rows and columns of an adjacency matrix."),
  random_answer_order = TRUE,
  allow_retry = TRUE
)
```

This single permutation suggests a more even mixing of these
attributes is possible, but it is just a single permutation.
Let's try a `r gloss("QAP","qap")` test that runs this over a succession of permutations,
just as we did with random graphs.
Plot the results for gender and power according to the random and permutation baselines.

```{r testperm, exercise=TRUE, exercise.setup = "rando-solution", purl = FALSE, exercise.timelimit = 120, fig.width=9}

```

```{r testperm-hint, purl = FALSE}
test_permutation(____, FUN = ____, attribute = ____,
                                times = ____)
```

```{r testperm-solution}
(perm.gender <- test_permutation(marvel_friends, 
                                net_heterophily, attribute = "Gender",
                                times = 1000))
(perm.power <- test_permutation(marvel_friends, 
                                net_heterophily, attribute = "PowerOrigin",
                                times = 1000))

(plot(rand.gender) + ggtitle("CUG test results for 'Gender' attribute") + 
    theme(plot.title = element_text(size=8)) | 
    plot(rand.power) + ggtitle("CUG test results for 'PowerOrigin' attribute") + 
    theme(plot.title = element_text(size=8))) /
(plot(perm.gender) + ggtitle("QAP test results for 'Gender' attribute") + 
   theme(plot.title = element_text(size=8)) | 
   plot(perm.power) + ggtitle("QAP test results for 'PowerOrigin' attribute") + 
   theme(plot.title = element_text(size=8)))
```

Again, we see that there is perhaps nothing so surprising that we got the homophily score
for gender that we did, but the lack of power origin heterophily is surprising.
Note how the distributions are generally wider when permuting the observed network 
than creating a random distribution (be mindful of the scale of the _x_-axis).
That is, taking into account the structure of the network leads us to
expect a larger spread in the EI index than when the variable is distributed across a random network.

```{r cupqap-results, echo=FALSE, purl = FALSE}
question("What can we say from these results?",
  answer("We can reject the null hypothesis that there is no homophily for gender in the comparison with random networks",
         message = "We cannot reject the null hypothesis since the result is in the range of expected scores across randomised networks, indicated by the red line between the grey shaded zones."),
  answer("We can reject the null hypothesis that there is no homophily for power origin in the comparison with random networks", 
         correct = TRUE,
         message = learnr::random_praise()),
  answer("We can reject the null hypothesis that there is no homophily for power origin in both CUG and QAP tests",
         correct = TRUE,
         message = learnr::random_praise()),
  answer("We can reject the null hypothesis that there is no homophily for gender in both CUG and QAP tests", 
         message = "We cannot reject the null hypothesis since the result is in the range of expected scores across randomised and permuted networks, indicated by the red lines between the grey shaded zones in both plots."),
  allow_retry = TRUE
)
```

CUG and QAP tests are useful for establishing whether there is
more or less homophily than expected by chance,
with QAP tests being more conservative and better taking into account the structure of the network.
Alternative reference networks can also be used that condition on more structural
information than completely random graphs, 
but introduce potentially more variation in the network structure than permutations.
These include for example `r gloss("configuration")` networks that retain the degree distribution observed in the original network, but could result in varying clustering,
or networks that are conditioned on the dyad census (i.e. `r gloss("MAN","man")` models).
And yet none of these approaches allow us to examine multiple variables at once
or control for other factors.
For that we need to move on to network regression.

## Network regression

Next let us examine homophily in another network.
The data were collected as part of an early experiment on communication between social network researchers who were using an Electronic Information Exchange System (EIES).
You may recognise some of the names.
The main network consists of 32 scholars with directed ties weighted by the total number of messages sent from $i$ to $j$ over the period of the study.
Nodal attributes collected include the primary discipline and number of citations in the social science citation index at the start of the study.

```{r introeies, exercise=TRUE, purl = FALSE, fig.width=9}
ison_networkers
graphr(ison_networkers, node_color = "Discipline")
```

Let's use both the continuous `Citations` and the categorical `Discipline` variables 
and come up with a couple of key hypotheses:

- $H_1$: researchers send more emails to those who are cited more (popularity hypothesis)
- $H_2$: researchers send more emails to those of the same discipline (disciplinary homophily hypothesis)

Let's start with a pretty maximally specified model 
(note that it doesn't make sense to include both `r gloss("ego")` and `r gloss("alter")` effects because these networks are undirected).
We are using `times = 200` here because of time-out limitations in the `learnr` tutorial system,
but for publication quality results you would want to base your conclusions on 2000
simulations or more.

```{r qapmax, exercise=TRUE, exercise.timelimit = 3600, purl = FALSE}

```

```{r qapmax-hint-1, purl = FALSE}
net_regression(____, ison_networkers, times = 200)
# If the model runs into a timeout error, please reduce the number of 'times' in the function above.
```

```{r qapmax-hint-2, purl = FALSE}
weight ~ ego(Citations) + alter(Citations) + sim(Citations) + 
                      ego(Discipline) + same(Discipline)
```

```{r qapmax-solution}
model1 <- net_regression(weight ~ ego(Citations) + alter(Citations) + sim(Citations) + 
                      ego(Discipline) + same(Discipline), 
                      ison_networkers, times = 200)
```

By default, the printing method for this model includes the key information you might need.
We can also use tidy methods to obtain just the table of coefficients (e.g. `tidy()`) 
or the model summary statistics (e.g. `glance()`).
`{autograph}` also includes a plot method for these results to
facilitate the quick interpretation of these results,
and as long as this is loaded you can simply call `plot()` on the model object.

```{r qapinterp, exercise=TRUE, exercise.setup = "qapmax-solution", purl = FALSE, fig.width=9}

```

```{r qapinterp-solution}
model1
tidy(model1)
glance(model1)
plot(model1)
```

This violin plot presents the distribution of coefficients from permutations of the network,
with the coefficient fitted from the data as a red dot.
Lines are used to indicate 95% thresholds,
but here the distributions are rendered so wide that they are often not visible.

```{r qap-interp, echo=FALSE, purl = FALSE}
question("What can we say from the results from model 1?",
  answer("Researchers send more messages to those who are cited more", 
         message = "Looks like alter Citations is not significant."),
  answer("Researchers send more messages to those who similarly cited", 
         message = "Looks like sim Citations is not significant."),
  answer("Researchers send more messages to mathematicians than anthropologists", 
         message = "Looks like alter Discipline Mathematics is not significant."),
  answer("Researchers send more messages to sociologists than anthropologists", 
         message = "Looks like alter Discipline Sociology is not significant."),
  answer("Not much", 
         correct = TRUE, 
         message = "Yes, the fitted coefficients are, if not typical, at least unsurprising from permutations of the network and so there's no evidence for rejecting the null hypothesis on the basis of this data."),
  allow_retry = TRUE
)
```

## Free play

While these are the conclusions from this 'play' data, 
you may have more and more interesting data at hand.
How would you go about specifying such a model?
Why is such an approach more appropriate for network data than linear
or logistic regression?

```{r freeplay, exercise = TRUE, fig.width=9}

```

## Glossary

Here are some of the terms that we have covered in this module:

`r print_glossary()`
