README.Rmd

---
title: "useR!2017 poster"
author: By David Selby and David Firth, Department of Statistics, University of Warwick.
output:
  github_document:
    df_print: kable
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      cache = TRUE,
                      message = FALSE,
                      fig.path = 'img/')
options(scipen = 999) # disable scientific notation

# Load system fonts (Arial Narrow) for PNG and PDF plots
extrafont::loadfonts(device = 'win', quiet = TRUE)
extrafont::loadfonts(device = 'pdf', quiet = TRUE)
```

A poster for the useR!2017 conference in Brussels. 

### Title

Ranking influential communities in networks

### Abstract

Which scientific fields export the most intellectual influence, through recent research, to other fields?
Citation behaviour varies greatly across disciplines, making inter-field comparisons difficult.
Recent approaches to measuring influence in networks, based on the PageRank algorithm, take the source of citations (or recommendations or links) into account.

By aggregating all publications in each Web of Science field into "super-journals", we model the exchange of citations between fields.
We then fit a Bradley–Terry paired comparisons model—which can be shown to be equivalent to scaled PageRank—to measure the relative influence of academic communities.
Uncertainty intervals are provided to quantify the robustness of the ranking.
All analyses were performed in R.

Which field is top of the table? Visit the poster to find out!

**Keywords**: PageRank, Bradley–Terry model, networks, ranking, bibliometrics

## Design inspiration

Alberto Cairo's [two-page spread on the *Galileo* telescope][galileo] is a good basis for a layout.
This features one huge diagram just left of centre, with an introductory paragraph on the far left, a self-contained series of fun facts along the bottom and a set of more detailed diagrams on the right.
See also [Giants of the Ocean][giants].

[galileo]: http://www.domusweb.it/content/dam/domusweb/en/interviews/2014/01/31/the_art_of_informationvisualization/rbig/Cairo2.jpg
[giants]: http://www.domusweb.it/content/dam/domusweb/en/interviews/2014/01/31/the_art_of_informationvisualization/rbig/Cairo1.jpg

See also the diagram in Figure 3 of *Fast unfolding of communities in large networks* by [Blondel et al (2008)](https://arxiv.org/abs/0803.0476), where we have a large node-link diagram with communities aggregated into nodes, then a 'zoom' on a selected community, represented as its constituent nodes and links.

Since useR!2017 specifies posters must be A0 in portrait orientation, we can't implement landscape layouts like those mentioned above.
While using some of those ideas, we might consider Colin Purrington's [results area design for portrait posters](https://i0.wp.com/colinpurrington.com/wp-content/uploads/2011/09/poster-template-vertical-2-purrington.jpg).

Another useful resource is the [Better Posters blog](http://betterposters.blogspot.co.uk/).

A large dose of inspiration may also be drawn from Dorling Kindersley *Eyewitness* books.

# Method

Unfortunately, some of the raw data is property of Thomson Reuters / Clarivate Analytics so I cannot republish it in this repository.
Nonetheless, hopefully you can get a general idea of what I have done from the source code and results.

We will need the following packages.

```{r packages, cache = FALSE, message = FALSE}
library(igraph)
library(ggraph)
```

## Data wrangling

First we load the 2013 Web of Science citation data into R.

```{r load}
WoS <- readRDS('../thesis/data/thomsonreuters.Rds')
```

We are only interested in citations from 2013 to publications from 2003–2012.
Furthermore, some miscellaneous journals are unhelpfully classified as `ALL OTHERS` and we want to omit any reference to these as well.
Citations of count (`weight`) equal to zero should be deleted to avoid connectivity problems later.

```{r filter}
WoS <- data.frame(from = WoS$Citing, to = WoS$Cited, weight = WoS$AllYears - WoS$Earlier)
WoS <- WoS[- union(grep('ALL OTHERS', WoS$from), grep('ALL OTHERS', WoS$to)), ]
WoS <- subset(WoS, weight > 0)
```

It's important to put the citing journal in the first column, the cited journal second and call the citation counts `weight`, because of how the `igraph` function `graph_from_data_frame` works.
We will now turn our data into an `igraph` object.

```{r igraph}
ig <- graph_from_data_frame(WoS)
```

We now have a weighted, directed graph.

Before we run our community detection algorithm, there is some other housekeeping to do.
A number of journals either gave out or received zero citations during the study period.
These will break the graph up into disconnected components, so we want to remove all these singletons.

Firstly, we calculate the strongly-connected components of the graph.
To be strongly connected to the rest of the network, a journal must both cite and be cited by others.

```{r strong}
strong <- components(ig, mode = 'strong')
core <- which(strong$csize > 1)
if (length(core) > 1) stop('There should be only one core component')
```

The whole graph contains `r scales::comma(vcount(ig))` nodes (journals).
Of these, `r scales::comma(strong$csize[core])` belong to our "core" strongly-connected component and `r scales::comma(sum(strong$csize == 1))` are singletons, either weakly connected or completely disconnected from the rest of the graph.
Let's get rid of these singleton nodes.

```{r singletons}
ig <- induced_subgraph(ig, which(strong$membership == core))
```

We now have `r scales::comma(vcount(ig))` journals in our network.

## Super journals

We will put our directed, weighted graph through the Infomap algorithm, as implemented in the `igraph` package.
Results can be nondeterministic, so we will fix the random seed for reproducible results.
Community detection can take a long time, so you might want to `cache` this chunk!

```{r infomap}
set.seed(2017)
infomap <- cluster_infomap(ig)
```

The algorithm returns `r length(sizes(infomap))` communities; the largest contains `r max(sizes(infomap))` journals and the smallest contains `r min(sizes(infomap))`.
The mean community size is `r round(mean(sizes(infomap)))` journals.

The journals of community will be aggregated into a *super-journal* representing all incoming and outgoing citations for that community.
Edge weights are summed and other edge and vertex attributes are ignored.

```{r aggregate}
sj <- contract.vertices(ig, membership(infomap), 'ignore')
sj <- simplify(sj, remove.multiple = TRUE, remove.loops = TRUE,
               edge.attr.comb = list(weight = 'sum', 'ignore'))
```

For later reference, each super-journal will be assigned a unique ID.

```{r id}
V(sj)$name <- 1:vcount(sj)
```

## Visualisation

A nice way to visualise these graphs is with the `ggraph` package.
Each node will be represented by a point, proportional to its PageRank score,
and each edge will be represented by an arc between these points.

There is a problem, however: `ggraph`, `ggplot2` and our PDF reader won't like it if we try to plot `r scales::comma(sum(E(sj)$weight))` arcs in a single graphic!

(Another approach is to have a single arc for each pair of nodes, with opacity or width proportional to the number of citations, but this doesn't look very good in my opinion.)

Calculating the positions of nodes in our graphic will involve multidimensional scaling of the correlation matrix.
Let's extract the weighted adjacency matrix, paying attention to the fact that `igraph` considers citations to travel from rows to columns.

```{r xtab}
xtab <- Matrix::t(as_adjacency_matrix(sj, attr = 'weight'))
```

We will scale the citations counts so no pair of nodes has more than 1000 arcs drawn between them on the final graphic.

```{r scaling}
small_xtab <- as.matrix(xtab)
diag(small_xtab) <- 0 # ignore self-citations
scalefactor <- 1000 / max(small_xtab)
small_xtab <- ceiling(small_xtab * scalefactor)
```

Then we can turn this back into an `igraph` object for visualisation.
A useful attribute is the PageRank, which we can use to make nodes with greater total influence appear larger on the plot.

```{r igraph_viz}
viz_ig <- graph_from_adjacency_matrix(t(small_xtab))
V(viz_ig)$PageRank <- page.rank(sj)$vector
```

```{r add_labels, echo = FALSE}
labels <- read.csv('data/cluster_names.csv', stringsAsFactors = FALSE)
V(viz_ig)$field <- V(sj)$field <- labels$field[match(V(viz_ig)$name, labels$community)]
```

Now let's perform multidimensional scaling to generate some coordinates.


```{r layout}
my_mds <- create_layout(viz_ig,
                        layout = 'igraph',
                        algorithm = 'mds',
                        dist = 1 - cor(as.matrix(xtab)))
```

Let's make some plots!

```{r superjournals, dev = 'cairo_pdf', fig.width = 10, fig.height = 10, fig.show = 'hide'}
ggraph(my_mds) +
  geom_edge_fan0(alpha = .01, colour = '#4F94CD') +
  geom_node_point(aes(size = PageRank), fill = '#4F94CD', pch = 21, colour = 'white') +
  coord_fixed() +
  scale_x_reverse() + # flip horizontally
  theme_graph() +
  theme(legend.position = 'none')
```

On its own, our visualisation does not imply much because we don't know which community is which.
We need labels on the communities for that.
I have gone through and manually assigned plausible labels to the communities generated by Infomap.

```{r ref.label = 'add_labels', echo = TRUE, eval = FALSE}
```

We can then have a labelled plot for reference, or even show a selection of "interesting" labels on the main graph, while omitting most of them to avoid clutter.

```{r labelled, dev = 'png', fig.width = 10, fig.height = 10, dev.args = list(type = 'cairo'), dpi = 600, warning = FALSE}
# Hide plot labels for very small fields
my_mds$flabel <- ifelse(rank(my_mds$PageRank) > 32, as.character(my_mds$field), NA)

ggraph(my_mds) +
  geom_edge_fan0(alpha = .007, colour = 'tomato2') + #4F94CD #2D7B95
  geom_node_point(aes(size = PageRank), fill = 'tomato2', pch = 21, colour = 'white') +
  geom_node_label(aes(label = flabel),
                  size = 2,
                  repel = TRUE,
                  family = 'Gill Sans MT Condensed',
                  fontface = 'bold',
                  colour = 'tomato2',
                  # Label options
                  segment.alpha = .5,
                  segment.size = 0.2,
                  label.r = unit(0.1, 'lines'),
                  label.size = NA, # no label border
                  label.padding = unit(0.1, 'lines'),
                  fill = rgb(1, 1, 1, .5)) +
  coord_fixed() +
  scale_x_reverse() +
  theme_graph() +
  theme(legend.position = 'none')
```

```{r labelled2, ref.label = 'labelled', fig.width = 10, fig.height = 10, include = FALSE, dev = 'cairo_pdf'}
```

## Within-field analysis

Now, let's take a particular field out of the network and examine its inner structure.
Varin et al. (2016) set a precedent to study statistics journals, so let's have a look at statistics.

```{r stats_subgraph}
statistics <- labels$community[labels$field == 'statistics']
stats_ig <- induced_subgraph(ig, which(membership(infomap) == statistics))
V(stats_ig)$PageRank <- page.rank(stats_ig)$vector
```

There are `r vcount(stats_ig)` journals in the statistics subgraph, in `r count_components(stats_ig, mode = 'strong')` strongly-connected component(s).

At the moment we only have Thomson Reuters' abbreviations for the journal names.
Some are more obvious than others! We have the full journal titles, but some are a bit too long, so we will use a table of custom short-ish titles.

```{r journaltitles}
titles <- read.csv('data/stats_titles.csv', stringsAsFactors = FALSE)
V(stats_ig)$title <- titles$Short[match(V(stats_ig)$name, titles$JCR)]
```

We can visualise the network as before.

```{r stats_viz_web, dev = 'png', fig.width = 10, fig.height = 10, dev.args = list(type = 'cairo'), dpi = 600}
stats_xtab <- Matrix::t(as_adjacency_matrix(stats_ig))
stats_layout <- create_layout(stats_ig,
                              layout = 'igraph',
                              algorithm = 'mds',
                              dist = 1 - cor(as.matrix(stats_xtab)))

ggraph(stats_layout) +
  geom_edge_fan0(alpha = .05, colour = 'tomato2') +
  geom_node_point(aes(size = PageRank), fill = 'tomato2', pch = 21, colour = 'white') +
  geom_node_text(aes(label = title), size = 3,
                 repel = TRUE,
                 family = 'Gill Sans MT Condensed',
                 fontface = 'bold',
                 colour = 'tomato2',
                 segment.alpha = .2) +
  coord_fixed() +
  theme_graph() +
  theme(legend.position = 'none')
```

```{r stats_viz, dev = 'cairo_pdf', fig.width = 10, fig.height = 10, include = FALSE, dev.args = list(bg = 'transparent')}
ggraph(stats_layout) +
  geom_edge_fan0(alpha = .05, colour = '#2D7B95') +
  geom_node_point(aes(size = PageRank), fill = '#2D7B95', pch = 21, colour = 'white') +
  geom_node_text(aes(label = title), size = 3,
                 repel = TRUE,
                 family = 'Gill Sans MT Condensed',
                 fontface = 'bold',
                 colour = '#2D7B95',
                 segment.alpha = .2) +
  coord_fixed() +
  theme_graph() +
  theme(legend.position = 'none',
        plot.background = element_rect(fill = 'transparent', colour = NA),
        panel.background = element_rect(fill = 'transparent', colour = NA))
```

And here is a ranking of these statistical journals, by Bradley–Terry score.

```{r stats_rank, message = FALSE}
(stats_ranks <- dplyr::arrange(
  data.frame(journal = V(stats_ig)$title,
             Scroogefactor = scrooge::Scroogefactor(stats_xtab),
             PageRank = V(stats_ig)$PageRank,
             BradleyTerry = scrooge::BTscores(stats_xtab),
             rank = rank(-scrooge::BTscores(stats_xtab), ties.method = 'first')),
  desc(BradleyTerry)
))
```

And now for something completely different.
What if we actually have a single non-statistics "super-journal" in the statistics network?

```{r stats_others}
mapping <- match(V(ig)$name, V(stats_ig)$name)
mapping[is.na(mapping)] <- max(mapping, na.rm = TRUE) + 1
stats_others <- contract.vertices(ig, mapping, list(name = 'first', 'ignore'))
V(stats_others)$name[vcount(stats_others)] <- "(All others)"
others_xtab <- Matrix::t(as_adjacency_matrix(stats_others))
diag(others_xtab) <- 0

dplyr::arrange(
  data.frame(journal = V(stats_others)$name,
             Scroogefactor = scrooge::Scroogefactor(others_xtab),
             PageRank = page.rank(stats_others)$vector,
             BradleyTerry = scrooge::BTscores(others_xtab)),
  desc(BradleyTerry)
)
```

We can try plotting it, too.

```{r plot_stats_others, dev = 'png', fig.width = 10, fig.height = 6, dev.args = list(type = 'cairo'), dpi = 600}
stats_others <- simplify(stats_others, remove.loops = TRUE, remove.multiple = FALSE) # remove self-citations
V(stats_others)$PageRank <- page.rank(stats_others)$vector

other_layout <- create_layout(stats_others,
                              layout = 'igraph',
                              algorithm = 'mds',
                              dist = 1 - cor(as.matrix(others_xtab)))

ggraph(other_layout) +
  geom_edge_fan0(alpha = .01, colour = '#4F94CD') +
  geom_node_point(aes(size = PageRank), fill = '#4F94CD', pch = 21, colour = 'white') +
  geom_node_text(aes(label = name), size = 3,
                 repel = TRUE,
                 family = 'Gill Sans MT Condensed',
                 fontface = 'bold',
                 colour = '#4F94CD',
                 segment.alpha = .2) +
  coord_fixed() +
  theme_graph() +
  theme(legend.position = 'none')
```

## Ranking

We can penalise the Bradley–Terry model by adding a "player zero" who cites/is cited by every player/journal/field at a constant rate (say 1/2). This will help reduce the chance of outliers (such as journals for which we have very little citation data, or fields containing very few journals) from shooting to the top or the bottom of a Bradley–Terry scores league table.

```{r ranking, dev = 'cairo_pdf', fig.show = 'hide', fig.width = 5, fig.height = 3, dev.args = list(bg = 'transparent')}
zero_cite <- .15 * sum(xtab) / 2 / nrow(xtab) ## 2ka = .15n
penalised_xtab <- rbind(zero_cite, cbind(zero_cite, xtab)) # add zeroth player

library(scrooge)
field_ranks <- data.frame(
  field = V(sj)$field,
  PageRank = PageRank(penalised_xtab)[-1],
  BradleyTerry = BTscores(penalised_xtab)[-1],
  Scroogefactor = Scroogefactor(penalised_xtab)[-1]
)
field_ranks$rank <- rank(-field_ranks$BradleyTerry, ties.method = 'first')
dplyr::arrange(field_ranks, desc(BradleyTerry))

lastplace <- max(field_ranks$rank)
interestingfields <- c('mathematics',
                       'informatics',
                       'biomedical sciences',
                       'medicine',
                       'chemistry & physics',
                       'zoology',
                       'textiles')

library(ggrepel)
ggplot(field_ranks) +
  aes(rank, 100*BradleyTerry, label = field) +
  geom_point(colour = 'tomato2', size = 1) +
  geom_text_repel(data = subset(field_ranks, rank < 5 | field %in% interestingfields),
                  nudge_y = .15,
                  nudge_x = -5,
                  segment.alpha = .25,
                  family = 'Gill Sans MT Condensed', 
                  fontface = 'bold',
                  colour = 'tomato2',
                  point.padding = unit(0.1, 'lines')
                  ) +
  scale_x_reverse(name = NULL,
                  labels = scales::ordinal,
                  breaks = c(1, seq(20, lastplace, by = 20)),
                  minor_breaks = c(1, seq(5, lastplace, by = 5)),
                  limits = c(lastplace, 1),
                  expand = c(0.02, 0)) +
  scale_y_continuous(NULL,
                     position = 'right') +
  theme_bw() +
  theme(text = element_text(family = 'Gill Sans MT Condensed', face = 'bold'),
        axis.text = element_text(colour = 'tomato2'),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.ticks = element_line(colour = 'tomato2'),
        axis.line = element_line(colour = 'tomato2'),
        plot.background = element_rect(fill = 'transparent', colour = NA),
        panel.background = element_rect(fill = 'transparent', colour = NA))
```

```{r stats_ranking, dev = 'cairo_pdf', fig.show = 'hide', fig.width = 5, fig.height = 3, dev.args = list(bg = 'transparent')}
interestingjournals <- c('R Journal',
                         'Stat Science',
                         'Biometrika',
                         #'J Stat Software',
                         'Annals',
                         'Stats in Medicine',
                         #'Biostatistics',
                         'Machine Learning',
                         'JRSS-B',
                         'Statistics',
                         #'J Applied Stats',
                         'Metrika',
                         'JCGS')
lastplace2 <- max(stats_ranks$rank)

ggplot(stats_ranks) +
  aes(rank, 100*BradleyTerry, label = journal) +
  geom_point(colour = '#2D7B95', size = 1) +
  geom_text_repel(data = subset(stats_ranks, rank <= 3 | journal %in% interestingjournals),
                  nudge_y = .15,
                  nudge_x = -5,
                  segment.alpha = .25,
                  family = 'Gill Sans MT Condensed',
                  fontface = 'bold',
                  colour = '#2D7B95',
                  point.padding = unit(0.2, 'lines')
                  ) +
  scale_x_reverse(name = NULL,
                  labels = scales::ordinal,
                  breaks = c(1, seq(20, lastplace2, by = 20)),
                  minor_breaks = c(1, seq(5, lastplace2, by = 5)),
                  limits = c(lastplace2, 1),
                  expand = c(0.02, 0)) +
  scale_y_continuous(NULL,
                     position = 'right') +
  theme_bw() +
  theme(text = element_text(family = 'Gill Sans MT Condensed', face = 'bold'),
        axis.text = element_text(colour = '#2D7B95'),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.ticks = element_line(colour = '#2D7B95'),
        axis.line = element_line(colour = '#2D7B95'),
        plot.background = element_rect(fill = 'transparent', colour = NA),
        panel.background = element_rect(fill = 'transparent', colour = NA))
```