This vignette performs and documents the analysis for the visualisation Measuring influence in scientific fields, available online. Unfortunately the raw data is not freely available, so only the methodology may be published for now.
Looking for the D3.js visualisation code? See http://selbydavid.com/influence/influence.js
We will use the following packages.
library(BradleyTerryScalable) # 0.1.0
library(dplyr)
library(ggplot2); theme_set(theme_bw())
library(ggrepel)
library(igraph)
library(qvcalc)
library(tidyr)
Citation data has been provided by Thomson Reuters/Clarivate Analytics. It includes data from Journal Citation Reports (JCR) year 2015.
JCR_raw <- read.table('data/2015_Cited.txt', stringsAsFactors = FALSE,
col.names = c('Cited', 'Citing', 'AllYears', 2014:2005, 'Earlier'))
We have 14 columns but are only interested in the identities of the citing and cited journals and the number of citations issued in 2015 to publications from the preceding 10 years (2005–2014). Any pairs of journals with no citations in this period may be omitted.
JCR <-
JCR_raw %>%
transmute(Citing, Cited, weight = AllYears - Earlier) %>%
filter(weight > 0)
We have arranged a data frame into a format suitable for conversion to an igraph
object: three columns, corresponding to “source”, “destination” and edge weight
.
G <- graph_from_data_frame(JCR)
Journals that give out zero citations or receive none themselves are not really interesting and can cause problems by making the network not to be strongly connected. We will only keep the main strongly-connected component of the graph.
components <- components(G, 'strong')
G <- delete_vertices(G, components$membership != which.max(components$csize))
There are two different techniques we might use to discover communities in the citation network: the directed Louvain algorithm and the map equation.
Respective authors have published software to run each algorithm, but they do not use a common interface; each requires input data to be in a slightly different format. Thus our next task is to format our citation graph into files appropriate for each programme. This is a bit annoying — hopefully in the future we can unify this into a single interface.
The following function converts an igraph
object into pajek
format and writes it to a file. The Louvain algorithm doesn’t require the vertices and arcs to be declared, but the Infomap algorithm does, so we have a slightly different format for each procedure.
write_pajek <- function(g, file = stdout(), louvain = FALSE) {
nvertices <- igraph::vcount(g)
output <- paste('*Vertices', nvertices)
output <- c(output, paste0(1:nvertices, ' "', V(g)$name, '"'))
output <- c(output, "*Arcs")
edges <- igraph::ends(g, 1:igraph::ecount(g))
# Convert vertex names to indices
edges <- apply(edges, 2, function(x) match(c(x), V(g)$name))
edges <- apply(edges, 1, paste, collapse = ' ')
if (louvain) # Only requires `src dest weight`
output <- NULL
output <- c(output, paste(edges, E(g)$weight))
if (is.character(file)) {
file <- file(file)
on.exit(close(file))
}
writeLines(output, file)
}
Using this function, we write the citations data to disk.
write_pajek(G, 'data/JCR_infomap.net')
write_pajek(G, 'data/JCR_louvain.txt', louvain = TRUE)
Now we can run the directed Louvain algorithm. This requires an extra step: first we have to convert from plaintext to the Louvain algorithm’s dedicated binary format. This input is then passed to the Louvain algorithm proper.
# Convert to binary
system('bin/convert -i data/JCR_louvain.txt -o data/JCR_louvain.bin -w data/JCR_louvain.weights')
# Run algorithm
system('./bin/community.exe data/JCR_louvain.bin -w data/JCR_louvain.weights -l -1',
intern = TRUE, ignore.stderr = TRUE) -> lv_output
The output ‘tree’ in lv_output
should now be processed into a suitable data frame. As our input file was indexed from one but C++ indexes from zero, the community.exe
programme has introduced a disjoint node, 0
which lives in its own category. We need to omit this. (Although the convert.exe
has an option -r
that supposedly renumbers nodes from zero, it doesn’t seem to work.)
louvain <- lapply(lv_output, strsplit, ' ') %>%
unlist() %>% as.numeric() %>%
matrix(byrow = TRUE, ncol = 2) %>%
as.data.frame() %>%
setNames(c('node', 'community')) %>%
filter(node != 0) %>% # remove nonexistent 0 node, added as an artifact by C++.
mutate(level = cumsum(node == 1)) # every time node counter resets, go to next level
A useful check: does our output have the same number of journals as the input did?
stopifnot(vcount(G) == nrow(filter(louvain, level == 1)))
There are 122 communities in level 1, 14 communities in level 2, 10 in level 3 and 10 in level 4. The latter corresponds to the maximum modularity solution. The following code walks down the hierarchy and finds the level 4 community to which each original node has been assigned.
louvain_communities <- louvain %>%
inner_join(filter(louvain, level == 2) %>%
select(-level, community2 = community),
by = c('community' = 'node')) %>%
inner_join(filter(louvain, level == 3) %>%
select(-level, community3 = community),
by = c('community2' = 'node')) %>%
inner_join(filter(louvain, level == 4) %>%
select(-level, community4 = community),
by = c('community3' = 'node')) %>%
filter(level == 1) %>% select(-level)
Now we need to match the communities up with the original journal names. A table of journal abbreviations and full names has been downloaded from the Journal Citation Reports web site.
JCR_abbreviations <- read.csv('data/2015_JCR_abbreviations.csv',
skip = 1,
stringsAsFactors = FALSE)[, 2:3]
The raw journal names are in a mix of all caps and title case, so we will use some tools to get everything looking a bit more presentable. The function capitaliseJournalTitle()
converts a character string to Title Case.
capitaliseJournalTitle <- function(y) {
y <- tolower(y)
unname(
sapply(y, function(x) {
x <- gsub('-', ': ', x)
s <- strsplit(x, " ")[[1]]
paste(toupper(substring(s, 1, 1)), substring(s, 2),
sep = "", collapse = " ")
})
)
}
Now, we join the directed Louvain output with the abbreviated journal titles, which are joined in turn with the full journal titles (converted to Title Case). The resulting data can be written to disk.
louvain_named <- louvain_communities %>%
mutate(journal = V(G)$name[node]) %>%
left_join(JCR_abbreviations, by = c('journal' = 'JCR.Abbreviated.Title')) %>%
mutate(Full.Journal.Title = capitaliseJournalTitle(Full.Journal.Title)) %>%
distinct(Full.Journal.Title, .keep_all = TRUE) # as left_join makes duplicates
louvain_named %>%
select(community1 = community,
community2:community4,
journal_abbr = journal,
journal_full = Full.Journal.Title) %>%
write.csv('output/louvain_communities.csv', row.names = FALSE)
Now, we can run the Infomap algorithm and save the output as follows. Adding -2
optimises for a two-level partition (i.e. community, node) without a deeper hierarchy. The -d
argument indicates the network has directed edges. The -k
argument takes self-loops into account (as they are ignored by default).
We can set the --num-trials
with -N
. This means we run the outer algorithm N
times and pick the run with the best solution. We can set the --seed
with -s
(for reproducibility).
The programme doesn’t let you specify the output filename; only the folder. To avoid having a convoluted folder structure, we will run the algorithm, saving the file as the same name each time to output/
, immediately renaming it using file.rename()
to something more sensible.
infomap_output_files <- paste0('output/infomap_',
c('undirected', 'directed', 'directed_self'),
'.tree')
if (any(!file.exists(infomap_output_files))) { # Only run if not run already
# Undirected
system('bin/Infomap data/JCR_infomap.net output -2u -N100 -s2018')
file.rename('output/JCR_infomap.tree', infomap_output_files[1])
# Directed, ignore self-loops
system('bin/Infomap data/JCR_infomap.net output -2d -N100 -s2018')
file.rename('output/JCR_infomap.tree', infomap_output_files[2])
# Directed, consider self-loops
system('bin/Infomap data/JCR_infomap.net output -2dk -N100 -s2018')
file.rename('output/JCR_infomap.tree', infomap_output_files[3])
}
stopifnot(!file.exists('output/JCR_infomap.tree')) # check nothing left behind
info_undirected <- read.table(infomap_output_files[1],
col.names = c('path', 'flow', 'name', 'node'))
info_directed <- read.table(infomap_output_files[2],
col.names = c('path', 'flow', 'name', 'node'))
info_directed_self <- read.table(infomap_output_files[3],
col.names = c('path', 'flow', 'name', 'node'))
Let’s see what the output looks like.
head(info_undirected)
path | flow | name | node |
---|---|---|---|
1:1 | 0.0197050 | PLOS ONE | 73 |
1:2 | 0.0076432 | P NATL ACAD SCI USA | 283 |
1:3 | 0.0060845 | SCI REP-UK | 46 |
1:4 | 0.0060106 | NATURE | 36 |
1:5 | 0.0053352 | SCIENCE | 307 |
1:6 | 0.0039210 | J BIOL CHEM | 559 |
The flow
column gives us a journal-level PageRank ranking of the network for free. Notice the top few journals are Nature, PNAS, Science, PLoS One et al.
The path
is hierarchical, representing the community and location within that community. This can be separated out into columns.
We have three different configurations: undirected, directed and directed with self-loops. The first one is included mainly for completeness, but to compare properly with Louvain we will assume we want one of the directed variants.
Since the hierarchical structure is stored as a string in the path
column, our first task is to separate this into multiple columns.
info_undirected <- info_undirected %>%
separate(path,
paste0('level', 1:2),
sep = ':',
convert = TRUE)
info_directed <- info_directed %>%
separate(path,
paste0('level', 1:2),
sep = ':',
convert = TRUE)
info_directed_self <- info_directed_self %>%
separate(path,
paste0('level', 1:2),
sep = ':',
convert = TRUE)
head(info_undirected)
level1 | level2 | flow | name | node |
---|---|---|---|---|
1 | 1 | 0.0197050 | PLOS ONE | 73 |
1 | 2 | 0.0076432 | P NATL ACAD SCI USA | 283 |
1 | 3 | 0.0060845 | SCI REP-UK | 46 |
1 | 4 | 0.0060106 | NATURE | 36 |
1 | 5 | 0.0053352 | SCIENCE | 307 |
1 | 6 | 0.0039210 | J BIOL CHEM | 559 |
Undirected Infomap returns 70 communities, directed Infomap returns 69 communities and directed Infomap with self-loops returns 450 communities. It looks like we should go for the directed Infomap algorithm, ignoring self-loops.
As with the Louvain algorithm, we need to bind the community membership structure with journal titles.
infomap_communities <- info_directed %>%
left_join(JCR_abbreviations, by = c(name = 'JCR.Abbreviated.Title')) %>%
mutate(Full.Journal.Title = capitaliseJournalTitle(Full.Journal.Title)) %>%
distinct(Full.Journal.Title, .keep_all = TRUE)
## Warning: Column `name`/`JCR.Abbreviated.Title` joining factor and character
## vector, coercing into character vector
A quick unit test to check that all the journals appear in the abbreviation translation table:
stopifnot(all(info_directed$name %in% infomap_communities$name))
Then save the data to disk.
infomap_communities %>%
select(community = level1,
journal_abbr = name,
journal_full = Full.Journal.Title) %>%
write.csv('output/infomap_communities.csv',
row.names = FALSE)
To rank communities against one another it is useful to combine communities into ‘super-journals’ and study the citations between them.
Let’s (optionally) save the Infomap citations to disk:
infomap_superjournals <- JCR %>%
inner_join(infomap_communities %>%
select(from = level1, name),
by = c(Citing = 'name')) %>%
inner_join(infomap_communities %>%
select(to = level1, name),
by = c(Cited = 'name')) %>%
group_by(from, to) %>%
summarise(weight = sum(weight)) %>%
ungroup()
# write.csv(infomap_superjournals,
# 'output/citations_infomap.csv',
# row.names = FALSE)
With a quick unit test that we have the right number of communities:
ncomms <- length(unique(infomap_communities$level1))
stopifnot(length(unique(infomap_superjournals$from)) == ncomms)
stopifnot(length(unique(infomap_superjournals$to)) == ncomms)
And the Louvain citations:
louvain_superjournals <- JCR %>%
inner_join(louvain_named %>% select(from = community, journal), by = c(Citing = 'journal')) %>%
inner_join(louvain_named %>% select(to = community, journal), by = c(Cited = 'journal')) %>%
group_by(from, to) %>%
summarise(weight = sum(weight)) %>%
ungroup()
# write.csv(louvain_superjournals,
# 'output/citations_louvain.csv',
# row.names = FALSE)
Now we have aggregated the communities into ‘super-journals’, let’s rank them. The remainder of this document is concerned with Infomap communities only.
We use the quasi-Stigler model described in Varin et al. (2016).
Stigler_superjournals <-
infomap_superjournals %>%
xtabs(weight ~ to + from, data = .) %>%
scrooge::BradleyTerry()
Stigler_superjournal_scores <-
Stigler_superjournals %>%
coef() %>%
c(0, .) %>%
setNames(1:ncomms) %>%
(function(x) x - mean(x))
One way to calculate approximate confidence intervals is with quasivariances, for which we shall use the qvcalc
package.
covariance_matrix <- matrix(0, ncomms, ncomms)
covariance_matrix[-1, -1] <- vcov(Stigler_superjournals)
## Warning in summary.glm(object, ...): observations with zero weight not used
## for calculating dispersion
interfield_quasivariances <- qvcalc(covariance_matrix,
estimates = Stigler_superjournal_scores)
First we need to filter the Journal Citation Reports down to the strongly-connected core.
JCR_core <- JCR %>%
filter(Citing %in% V(G)$name,
Cited %in% V(G)$name)
Then we bind the individual journals in the JCR data to their respective communities.
JCR_Infomap <- JCR_core %>%
left_join( infomap_communities %>%
select(Citing = name, fromCommunity = level1) ) %>%
left_join( infomap_communities %>%
select(Cited = name, toCommunity = level1) ) %>%
select(Cited, Citing, weight, fromCommunity, toCommunity)
## Joining, by = "Citing"
## Joining, by = "Cited"
Then we can fit a Bradley–Terry model to each community. This uses Ella Kaye’s BradleyTerryScalable package.
The following function attempts to fit a quasi-Stigler model to each of the 69 communities.
tryBradleyTerry <- function(d, player0 = TRUE) {
union(pull(d, Cited), pull(d, Citing)) %>%
unique %>% length -> nUniqueJournals
if ( nUniqueJournals < 2 )
return(NA)
tryCatch(
btfit(btdata(d), a = 1),
error = function(e) e
)
}
Once each model has been fitted (or an error returned, in the case of singleton communities), we can then extract the model coefficients (else return NA
).
tryBTscores <- function(model) {
if (!inherits(model, 'btfit'))
return (NA)
tryCatch({
if ( is.list(coef(model)) )
unname(coef(model)) # removes spurious prefixes like '1.PHILOS PUBLIC AFF'
# which can occur if network does not appear to be fully connected
else
coef(model)
},
error = function(e) e)
}
Algorithms such as PageRank have a damping factor that ensures ergodicity when networks may not be well-connected. This is a concern for us, too. How can we make sure that the quasi-Stigler model fitted to each community is stable?
By introducing a ‘player 0’ or ‘journal 0’, who wins and loses against every other player a constant number of times. This effectively adds shrinkage to the model; any journals with a very small number of citations given or received are likely to see their scores shrink towards the mean (i.e. 0).
For a feeling of familiarity with PageRank, but not for any particularly theoretically-justified reason, we will use the parameter 0.15, so that each player zero cites and is cited \[0.15 \times \frac{\sum_{i=1}^n \sum_{j=1}^n c_{ij}} {2n}\] times, where \(c_{ij}\) is the number of citations from journal \(i\) to journal \(j\), and \(n\) is the total number of journals in that community (excluding the 0th journal, of course).
player0_weights <- JCR_Infomap %>%
filter(fromCommunity == toCommunity) %>% # Within-field citations only
group_by(community = fromCommunity) %>%
select(-fromCommunity, -toCommunity) %>%
summarise(weight = .15 * sum(weight) / length(union(Cited, Citing)) / 2)
inPlayer0Weights <- seq_len(max(JCR_Infomap$fromCommunity)) %in% player0_weights$community
if (!all(inPlayer0Weights))
warning('Community ', which(!inPlayer0Weights), ' not in `player0_weights`')
## Warning: Community 69 not in `player0_weights`
player0_grid <- expand.grid(
Cited = 'Player 0',
Citing = union(JCR_Infomap$Cited, JCR_Infomap$Citing)
) %>%
left_join(select(infomap_communities, community = level1, name),
by = c(Citing = 'name')) %>%
rbind(select(., Cited = Citing, Citing = Cited, community)) %>%
left_join(player0_weights, by = 'community')
## Warning: Column `Citing`/`name` joining factor and character vector,
## coercing into character vector
Especially interesting is the introduction of a super-journal called ‘All others’, representing the rest of the network. The magnitude of the shift in journal ranking this induces can be used as some measure of how inter-disciplinary (or insular) a journal is. We won’t bother with a ‘player 0’ in this case, since the ‘All others’ journal will effectively play a similar role.
It is safe to ignore the number of self-citations that ‘All others’ has, since they will cancel out in the model anyway.
other_weights <- JCR_Infomap %>%
filter(fromCommunity != toCommunity) %>%
group_by(Cited) %>%
summarise(Citing = 'Other fields',
weight = sum(weight),
toCommunity = unique(toCommunity),
fromCommunity = toCommunity) %>%
bind_rows(
JCR_Infomap %>%
filter(fromCommunity != toCommunity) %>%
group_by(Citing) %>%
summarise(Cited = 'Other fields',
weight = sum(weight),
fromCommunity = unique(fromCommunity),
toCommunity = fromCommunity)
)
Now we fit the models.
Stigler_models <- JCR_Infomap %>%
bind_rows(other_weights) %>%
# Intra-community citations only:
filter(fromCommunity == toCommunity) %>%
group_by(community = fromCommunity) %>%
# The two community cols are now redundant:
select(-fromCommunity, -toCommunity) %>%
# Add `player 0` to each community:
bind_rows(player0_grid) %>%
# Fit quasi-Stigler models with/without player 0
do(Stigler = select(., -community) %>%
filter(Cited != 'Player 0',
Citing != 'Player 0',
Cited != 'Other fields',
Citing != 'Other fields') %>%
tryBradleyTerry,
Stigler0 = select(., -community) %>%
filter(Cited != 'Other fields',
Citing != 'Other fields') %>%
tryBradleyTerry,
StiglerXF = select(., -community) %>% # XF = 'external fields'
filter(Cited != 'Player 0',
Citing != 'Player 0') %>%
tryBradleyTerry) %>%
mutate(BTscores = list(tryBTscores(Stigler)),
BTscores0 = list(tryBTscores(Stigler0)),
BTscoresXF = list(tryBTscores(StiglerXF)))
## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector
See Firth and de Menezes (2004). ‘Quasi-variances’. Biometrika, 91, 1, pp. 65–80. The qvcalc R package is not well optimised for large datasets. We can get around this by writing our own quasi-variance loss functions in C++ and minimising using BFGS.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double logloss2 (NumericVector q, NumericMatrix B) {
int N = B.nrow();
double v, e, loss = 0.0;
for (int i = 0; i < N; ++i) {
for (int j = 0; j < i; ++j) {
v = B(i, i) + B(j, j) - 2 * B(i, j); // Variance of simple contrast
e = (q[i] + q[j]) / v; // Relative error
loss += log(e) * log(e); // Squared log loss
}
}
return loss;
}
// [[Rcpp::export]]
NumericVector Dlogloss2 (NumericVector q, NumericMatrix B) {
int N = B.nrow();
double v, e;
NumericVector grad(N);
for (int i = 0; i < N; ++i) {
for (int j = 0; j < i; ++j) {
v = B(i, i) + B(j, j) - 2 * B(i, j);
e = (q[i] + q[j]) / v;
grad[i] += 2 * log(e) / (q[i] + q[j]);
}
}
return grad;
}
With the functions defined, let’s optimise:
tryQV <- function(model) {
if (!inherits(model, 'btfit'))
return(NA)
vCov <- as.matrix(vcov(model))
q0 <- diag(vCov)
qvars <- optim(q0, logloss2, method = 'BFGS', gr = Dlogloss2, B = vCov)$par
return (sqrt(qvars)) # But sometimes quasi-variances are negative!
}
tryQV2 <- function(model) {
# Return standard errors of differences with 'Other fields'
# as a reference point. Because quasi-variances unstable.
# Var(mu1 - mu2) = Var(mu1) + Var(mu2) - 2 * Cov(mu1, mu2)
# where Var('Other fields') = Cov('Other fields', .) = 0
if (!inherits(model, 'btfit'))
return(NA)
vCov <- vcov(model, ref = 'Other fields')
SEs <- sqrt(Matrix::diag(vCov))
return(SEs)
}
Stigler_models <- Stigler_models %>%
mutate(quasiSE0 = list(tryQV(Stigler0)),
quasiSEother = list(tryQV2(StiglerXF)))
Quasi-variances don’t work so well at approximating the comparison intervals when the ‘Other fields’ superjournal is included. Therefore for the ‘wider influence’ case, we use tryQV2
, which calculates the variances of comparisons relative to the ‘Other fields’ superjournal. This means our ‘wider influence’ comparison intervals are only valid for comparisons with the ‘Other fields’ journal; not between regular journals. At least we don’t end up with negative quasi-variances!
Now we need to do a bit of housekeeping. Where communities contained a single journal, no model has been fitted and instead NA
has been returned. Also, some communities are not strongly connected, in which case the dangling nodes have been omitted without coefficients. This doesn’t really apply to the damped case, because all the nodes should have been artificially connected by player 0.
# Export scores
Stigler_scores <-
Stigler_models %>%
rowwise() %>%
mutate(BTscores = list(unlist(BTscores))) %>%
pull(BTscores) %>% unlist() %>%
data_frame(journal = names(.),
influence = .)
Stigler_scores0 <-
Stigler_models %>%
pull(BTscores0) %>% unlist() %>%
data_frame(journal = names(.),
influence = .)
# Quasi variances
Stigler_qvars0 <-
Stigler_models %>%
pull(quasiSE0) %>% unlist() %>%
data_frame(journal = names(.),
quasiSE = .)
Parsing the ‘all other’ journals requires some extra work, because it is a journal which has the same name in every community.
Stigler_scores_other <-
Stigler_models %>%
pull(BTscoresXF) %>% unlist() %>%
data_frame(journal = names(.),
influence2 = .)
Stigler_qvars_other <-
Stigler_models %>%
pull(quasiSEother) %>% unlist() %>%
data_frame(journal = Stigler_scores_other$journal,
quasiSE2 = .)
# For binding to the end of the results
external_influence <- Stigler_scores_other %>%
filter(journal == 'Other fields') %>%
cbind(Stigler_qvars_other %>%
filter(journal == 'Other fields') %>%
select(-journal)) %>%
mutate(memberOf = row_number())
What effect does player 0 have on scores?
Stigler_scores0 %>%
left_join(Stigler_scores, by = 'journal', suffix = c('.Player0', '.mle')) %>%
ggplot(aes(influence.mle, influence.Player0)) +
geom_point(alpha = .5) +
labs(x = 'Ordinary Stigler-model export scores',
y = 'Damped with \'Player 0\'',
caption = 'Damping set to .15*N/d/2')
## Warning: Removed 117 rows containing missing values (geom_point).
And the ‘All others’ journals?
Stigler_scores_other %>%
left_join(Stigler_scores, by = 'journal') %>%
ggplot(aes(influence, influence2)) +
geom_point(alpha = .5) +
labs(x = 'Within-field Stigler-model export scores',
y = 'Including inter-field citations')
## Warning: Removed 115 rows containing missing values (geom_point).
Though our data are arguably hierarchical, we will produce a ‘tidy’ structure, where each row describes the super-journal to which that journal belongs, and where every super-journal belongs to a pseudo-super-journal 0.
This file will be saved to csv so it can be read into d3.js
, or other visualisation software.
interfield <- data_frame(
journal = names(Stigler_superjournal_scores),
influence = Stigler_superjournal_scores,
quasiSE = interfield_quasivariances$qvframe$quasiSE,
memberOf = 0L
)
intrafield <- infomap_communities %>%
select(journal = name,
memberOf = level1,
full_name = Full.Journal.Title) %>% # = NA for superjournals
left_join(Stigler_scores0, by = 'journal') %>%
left_join(Stigler_qvars0, by = 'journal') %>%
left_join(Stigler_scores_other, by = 'journal') %>%
left_join(Stigler_qvars_other, by = 'journal') %>%
bind_rows(external_influence)
if (nmissing <- sum(is.na(intrafield$influence))) {
warning(nmissing, ' journal/s have no influence score\n',
sum(is.na(intrafield$quasiSE)), ' journal/s have no quasi-variance\n')
}
## Warning: 70 journal/s have no influence score
## 70 journal/s have no quasi-variance
output <- bind_rows(interfield, intrafield) %>%
group_by(memberOf) %>%
mutate(rank = dense_rank(desc(influence)),
rank2 = dense_rank(desc(influence2))) %>%
ungroup()
write.csv(output, 'docs/data/rankings.csv', row.names = FALSE)