Example: Topic modelling Covid preprints
Source:vignettes/covid-preprint-topics.Rmd
covid-preprint-topics.RmdIntroduction
This is an example topic modelling analysis using the Structural Topic Model (STM) (Roberts, Stewart, and Tingley 2019) to explore a subset of bioRxiv and medRxiv preprints covering research related to Covid-19. The medRxiv website actually provides access to a dedicated collection of preprints on both medRxix and bioRxiv covering COVID-19 SARS-CoV-2. Here we will however work with the whole preprint collection and create a subset using keyword matching.
Note that this analysis is intended to illustrate the method of topic modeling. A complete analysis of the sample preprints would be more detailed and would likely feature a larger number of topics.
This document assumes some familiarity with Structural Topic
Modeling (STM). Please consult the stm package
vignette (Roberts, Stewart, and
Tingley 2019) for background.
Getting the document data
The bioRxiv and medRxiv servers provide API access to retrieve
preprint meta-data. One option to access the API is to use the medrxivr
package. The code snippet below would for example collect all
bioRxiv preprint meta-data for the year 2015.
library(medrxivr)
biorxiv_raw <- mx_api_content(server = "biorxiv",
from_date = "2015-12-01",
to_date = "2015-12-31")This package has been setup to support exploring different thematic
preprint subsets and assumes usage of all preprints from both servers
published from 2013 to 2023. This data can not be shared in this package
as the usage and API terms of
bioRxiv/medRxiv do not permit redistribution and rehosting of
the complete data; therefore this package provides scripts to
collect and clean the data in order to replicate this analysis (check
/data-raw in the package repository).
The documented examples assume that you have obtained the data yourself and created a local copy for replication of the analyses. Below is the code that would allow you to do this. NOTE that at the time of writing this returned a dataset with 365526 records. Retrieving this data would take several hours.
When replicating an analysis like this use the bioRxiv/medRxiv API responsibly. Consider collecting smaller datasets (the example shown below would for example only require preprints from 2020 to 2023). Alternatively, use the bulk snapshot as detailed here.
library(dplyr)
library(medrxivr)
# get publications from medRxiv and bioRxiv
pubs_biorxiv_raw <- medrxivr::mx_api_content(server = "biorxiv",
#from_date = "2019-01-01",
to_date = "2023-12-31")
pubs_medrxiv_raw <- medrxivr::mx_api_content(server = "medrxiv",
#from_date = "2019-01-01",
to_date = "2023-12-31")
pubs_biorxiv_raw <- pubs_biorxiv_raw %>%
mutate(server = "biorxiv")
pubs_medrxiv_raw <- pubs_medrxiv_raw %>%
mutate(server = "medrxiv")
preprints_raw <- dplyr::bind_rows(pubs_biorxiv_raw, pubs_medrxiv_raw)
save(preprints_raw, file = "./data-raw/preprints_raw.Rdata")In addition to the 15 meta-data variables returned via the API the
code above adds server as a variable to indicate the origin
of the preprint. This is an important variable for the examples shown
below.
Cleaning, filtering and annotating the data
As a first step the raw preprint data has to be cleaned.
Specifically, we want to retain only one unique record per preprint,
multiple versions for each unique doi may
exist. Below we retain only the latest version per DOI and also
ensure that there are no duplicate dois. This is crucial as
we need a unique identifier for each document that we want to include in
the topic modelling.
library(dplyr)
preprints_cleaned <- preprints_raw %>%
group_by(doi) %>%
filter(version == max(version)) %>%
ungroup() %>%
distinct(doi, .keep_all = TRUE)Next we annotate the data with additional variables, specifically we
want to reduce the publication date to the publication
year, and add an additional variable
is_published, which is inferred from the
published variable. The latter provides the DOI of
the journal where the preprint has been published after peer review.
We also limit the data to preprints from 2020 to 2023 (for an analysis of Covid topics preprints prior to 2020 are not relevant), and only retain variables that we may want to explore as document covariates, i.e. document properties that potentially influence the prevalence of topics.
preprints <- preprints_cleaned %>%
mutate(published = stringr::str_trim(published)) %>%
mutate(published = na_if(published, "NA")) %>%
mutate(is_published = as.numeric(!is.na(published))) %>%
mutate(is_published = case_when(is_published == 1 ~ "published",
is_published == 0 ~ "not published",
TRUE ~ "undefined")) %>%
mutate(year = lubridate::year(date)) %>%
filter(year >= 2020 & year <= 2023) %>%
select(doi, server, title, abstract, date, year, version, is_published)Finally, we define keywords and reduce our preprint set to preprints that contain one of those keywords in either the title or abstract, resulting in a subset of 29692 publications.
library(stringr)
keywords <- c("sars-cov", "covid")
search_pattern <- stringr::regex(paste(keywords, collapse = "|"),
ignore_case = TRUE)
covid_preprints <- preprints %>%
filter(stringr::str_detect(title, pattern = search_pattern) |
stringr::str_detect(abstract, pattern = search_pattern))Preparing and preprocessing the documents for text analysis
In order to analyze the preprints with the stm package we
need to create a representation of the documents and document meta-data
that stm can utilize. stm offers built-in
methods to support this (see specifically the
textProcessor() and prepDocuments() methods
(Roberts, Stewart, and
Tingley 2019)). Here we use instead the quanteda package
(Benoit et al. 2018)
which provides a broad range of methods for text pre-processing and
analysis, creating formats also supported by stm.
Create a corpus
First, we create a corpus object from the dataframe of preprints. The corpus is essentially a library of documents that will be used for the next steps. It specifies which variable should be used to uniquely identify documents and which variable holds the textual content (here the preprint abstracts) that should be processed.
Echoing the corpus will provide some basic information. All other variables in the original dataframe will be interpreted and included as document metadata (‘docvars’), which could later be included in the STM topic modelling process.
Tokenize and preprocess
For further analysis the corpus documents have to be
tokenized, i.e. further processing the texts have to be
broken into semantic units that are relevant for our analysis. The most
common approach is to interpret each word (typically designated
by whitespaces or punctuation) as a token. This is applied here as well.
quanteda offers several alternative approaches. Instead of
individual words, sequences of words (n-grams) could for
example be used.
The tokenization method also provides several options for preprocessing and filtering the tokens. Here for example, while tokenizing, we will simultaneously remove punctuation, numbers, special symbols and URLs. Furthermore, we split words containing hyphens, a word like social-ecological will thus be split into two individual tokens (social and ecological).
The text preprocessing choices could strongly influence the results of a text analysis and should be thoroughly explained, carefully evaluated and ideally be based on theory (see (Denny and Spirling 2018)).
Create a Document-feature matrix
The tokens object is then used to create a
document-feature matrix. For further statistical analysis this
reduces the tokens to a matrix of documents (rows) and unique
terms (columns) that counts the number of occurrences
for each term in each document. quanteda captures this as
features which supports more general options than
terms (see the quanteda documentation for
details).
Filter terms and documents
This is followed by other (optional) processing and filtering steps. A common option for example — to reduce the size of the data or assist in the interpretation — is the removal of so-called stopwords (e.g. “the”, “and”, “or” etc).
A step omitted here is reducing words (terms) to their word stem. The stemming algorithm (several are available) reduces words to its word stem. The terms “universal”, “university” and “universe” would for example be reduced to the same word stem of “univers”; this example indicates that this approach may require careful consideration.
Stemming has the advantage that it could potentially reduce the size of the matrix substantially.
pubs_dfm <- pubs_dfm %>%
quanteda::dfm_remove(pattern = quanteda::stopwords("english")) #%>%
#quanteda::dfm_wordstem()# echo the result
> pubs_dfm
Document-feature matrix of: 29,692 documents, 82,472 features (99.87% sparse) and 6 docvars.
features
docs nitric oxide synthesised three isoforms synthases viz nnos neurons enos
10.1101/038398 6 6 1 1 1 1 1 1 2 1
10.1101/058511 0 0 0 0 0 0 0 0 0 0
10.1101/292979 0 0 0 2 0 0 0 0 0 0
10.1101/402370 0 0 0 0 0 0 0 0 0 0
10.1101/420737 0 0 0 0 0 0 0 0 0 0
10.1101/596700 0 0 0 0 0 0 0 0 0 0
[ reached max_ndoc ... 29,686 more documents, reached max_nfeat ... 82,462 more features ]Further options may be considered to reduce noise and/or the size of the matrix. The following code removes for example terms (or features) that consist only of one character, terms that do not appear in at least two different documents, and furthermore would remove documents that do not contain at least 5 tokens. In our example this drops one document and reduces the number of retained features by more than half.
pubs_dfm <- pubs_dfm %>%
quanteda::dfm_remove(min_nchar = 2) %>%
quanteda::dfm_trim(min_docfreq = 2, docfreq_type = "count") %>%
quanteda::dfm_subset(quanteda::ntoken(.) > 4)# echo the result
> pubs_dfm
Document-feature matrix of: 29,691 documents, 37,093 features (99.72% sparse) and 6 docvars.
features
docs nitric oxide synthesised three isoforms synthases viz neurons enos endothelial
10.1101/038398 6 6 1 1 1 1 1 2 1 2
10.1101/058511 0 0 0 0 0 0 0 0 0 0
10.1101/292979 0 0 0 2 0 0 0 0 0 0
10.1101/402370 0 0 0 0 0 0 0 0 0 0
10.1101/420737 0 0 0 0 0 0 0 0 0 0
10.1101/596700 0 0 0 0 0 0 0 0 0 0
[ reached max_ndoc ... 29,685 more documents, reached max_nfeat ... 37,083 more features ]Topic modeling
Fitting the STM topic model
The key input to any topic modelling algorithm is the number of
topics (K) that the model should be fit to, which we here
set to 20 (see the separate document for a discussion on suitable
choices for the number of topics).
Before fitting the topic model we convert the document-feature matrix
into the native STM format. In order to fit a topic model with the
stm() function we need the set of documents,
the vocabulary of which these documents are composed and a
dataframe specifying the values of all document meta-data variables
(data) which can be used in the process as “covariates”
that might influence the prevalence of topics in a document.
In the example below, we ask stm to incorporate the
origin of the document (server) and the publication
year when fitting the topic model. The argument
prevalence = ~ server * s(year) expresses that we assume
that the prevalence of topics in a document is influenced by these two
variables, and that they also interact, i.e. we work with the hypothesis
that different temporal trends could be expected for documents published
on either of the two preprint servers1.
The consideration of covariates is optional. If omitted the model reduces to a Correlated Topic Model (Blei and Lafferty 2007; Roberts, Stewart, and Tingley 2019).
We also supply a seed, which allows to replicate the
results of the topic modeling.
Estimating the effect of document covariates
Once the model has converged we can estimate the effect of document
covariates on the topic prevalence. The estimateEffect()
function allows to run regressions based on the formula specified as the
first argument. It is here identical to the formula used when fitting
the topic model, and regressions are run for all 20 topics. The same
metadata as used previously needs to be supplied for this function in
addition to the topic model object.
covid_effect_K20 <- estimateEffect(1:20 ~ server * s(year),
stmobj = covid_model_K20,
metadata = covid_stm_docs$meta)This concludes fitting the model. The following sections step through a sample exploration of this topic model.
Analysing and interpreting the topic model
Basic topic model information
The topic model is defined by two matrices that capture probability
distributions of topics over documents (gamma matrix) and words
(or terms) over topics (beta matrix). We can start exploring
these with some of the built-in functions of stm.
The plot() function plots a chart showing topic
proportions for all topics in the model. A topic is identified by a
unique ID (1-20) and in the plot below the five words (or terms) that
have the highest probability of being associated with the given topic.
This gives an early indication of the distinct latent topics in the
analysed subset of preprints.
plot(covid_model_K20, n = 5)
The summary() function provides a more detailed view of
the topics and can help to begin interpreting and labeling the 20
topics. Specifically, the output shows four different sets of words
associated with a topic. ‘Highest Prob’ lists the words that
have the highest probability of being associated with a topic. A
comparison of different topics highlights that a term such as
covid has a high probability for several topics. The list of
‘FREX’ words summarizes words that are frequent and exclusive
in a topic, i.e. characterize a topic in comparison to other topics
(consult stm::labelTopics() for details as well as
Lift and Score word sets).
summary(covid_model_K20)
># A topic model with 20 topics, 29691 documents and a 37093 word dictionary.
># Topic 1 Top Words:
># Highest Prob: model, can, covid, transmission, epidemic, data, disease
># FREX: npis, mathematical, compartmental, scenarios, seir, reproduction, sir
># Lift: 1we, abms, ao_scplowbstractc_scplowas, apt, artefact, asilv, asymptotically
># Score: distancing, social, epidemic, reproduction, npis, model, r0
># Topic 2 Top Words:
># Highest Prob: cov, sars, drug, antiviral, activity, drugs, covid
># FREX: src, figdir, o_linksmallfig, c_fig, m_fig, o_fig, gif
># Lift: k777, gif, pmmov, wwtps, 13k, 17k, 18k
># Score: mpro, antiviral, drug, inhibitors, protease, drugs, compounds
># Topic 3 Top Words:
># Highest Prob: covid, risk, age, mortality, ci, associated, years
># FREX: hispanic, ethnicity, pregnant, racial, black, smoking, preterm
># Lift: 65s, asmr, assault, backgroundethnic, backgroundracial, backgroundsocio, brunt
># Score: ci, age, women, mortality, ethnicity, aor, hispanic
># Topic 4 Top Words:
># Highest Prob: covid, health, pandemic, mental, social, study, survey
># FREX: loneliness, emotional, attitude, depression, insecurity, mental, anxiety
># Lift: insecurity, accelerometers, amhara, angry, anovas, asd, asleep
># Score: mental, anxiety, depression, respondents, social, psychological, students
># Topic 5 Top Words:
># Highest Prob: sars, cov, infection, testing, transmission, cases, children
># FREX: ifr, seroprevalence, schools, school, contacts, household, attack
># Lift: inmates, 19y, 35y, 39y, 9a, abidjan, addscovid
># Score: school, seroprevalence, children, schools, household, transmission, testing
># Topic 6 Top Words:
># Highest Prob: cov, sars, virus, viral, coronavirus, respiratory, infection
># FREX: bats, cats, deer, animals, covs, wildlife, hcov
># Lift: cats, 5x106, aav6, aegyptiacus, aethiops, affinis, agm
># Score: cov, sars, mice, rna, viruses, coronaviruses, animals
># Topic 7 Top Words:
># Highest Prob: patients, covid, hospital, disease, clinical, severe, admission
># FREX: acei, admission, arbs, admitted, icu, aki, aceis
># Lift: 1.1x109, 2020r1g1a1a01006229, 2l, 4.0x109, 40y, ahmad, ahrq
># Score: patients, admission, icu, hospital, admitted, ci, hospitalized
># Topic 8 Top Words:
># Highest Prob: covid, cases, countries, number, deaths, data, pandemic
># FREX: cfr, italy, cities, country, countries, fatality, china
># Lift: 1000m, 1th, 55th, abysmally, abyss, adhanom, adminstat
># Score: countries, cases, lockdown, deaths, country, cfr, daily
># Topic 9 Top Words:
># Highest Prob: protein, binding, spike, sars, cov, ace2, rbd
># FREX: conformational, cryo, conformation, glycans, conformations, nanobodies, residues
># Lift: 13c, 6lzg, 6m0j, 6vw1, 6vxx, aabpu, abdab
># Score: binding, rbd, protein, spike, ace2, proteins, epitopes
># Topic 10 Top Words:
># Highest Prob: data, can, learning, covid, using, model, based
># FREX: aerosol, n95, aerosols, respirators, decontamination, airborne, machine
># Lift: elastomeric, forehead, papr, radiomics, exhaled, singing, 0.3m
># Score: learning, masks, aerosol, machine, respirators, n95, mask
># Topic 11 Top Words:
># Highest Prob: vaccine, vaccination, covid, middle, vaccines, dot, dose
># FREX: hesitancy, dot, vaccinate, middle, hesitant, rollout, ve
># Lift: #949850, acceptant, adjrr, aesis, amparo, analysesthe, andersen
># Score: vaccination, vaccine, dot, booster, dose, vaccinated, middle
># Topic 12 Top Words:
># Highest Prob: studies, care, covid, health, research, data, pandemic
># FREX: reviews, telemedicine, preprints, scoping, articles, blacksquare, publications
># Lift: preprints, 1.2m, aas, abbreviating, accustomed, activists, advisor
># Score: review, care, services, articles, reviews, pubmed, service
># Topic 13 Top Words:
># Highest Prob: sars, cov, genome, mutations, sequencing, viral, variants
># FREX: phylogenetic, gisaid, clades, wgs, genomes, genomic, haplotype
># Lift: clades, snvs, 1.1.7s, 11083g, 14408c, 17del, 20a
># Score: mutations, genome, wastewater, genomes, sequences, sequencing, genomic
># Topic 14 Top Words:
># Highest Prob: cells, cell, sars, cov, expression, infection, ace2
># FREX: autophagy, mirnas, mirna, at2, ifns, ciliated, scrna
># Lift: 25hc, angiotensinogen, antagonizes, apcs, arf6, asgr1, at2s
># Score: cells, expression, ace2, cell, genes, epithelial, tmprss2
># Topic 15 Top Words:
># Highest Prob: antibody, sars, cov, antibodies, igg, responses, vaccine
># FREX: iga, bau, igg, humoral, immunogenicity, as03, reactogenicity
># Lift: 1x1011, 28d, 30ug, ad26cov2, addas03, adhu5, atellica
># Score: igg, antibody, antibodies, neutralizing, rbd, vaccine, spike
># Topic 16 Top Words:
># Highest Prob: sars, cov, pcr, samples, rt, test, testing
># FREX: ag, rdt, rdts, lod, rt, panbio, kits
># Lift: poct, cobas, panbio, #yomecorono, 1.6x104, 10min, 15min
># Score: rt, pcr, assay, samples, saliva, detection, assays
># Topic 17 Top Words:
># Highest Prob: variants, omicron, variant, delta, ba, cov, sars
># FREX: omicron, ba, xbb, subvariants, delta, bq, voc
># Lift: 1.5s, 129s2, 1f11, 2.86s, 3b8, 417n, 75d30121c11061
># Score: omicron, ba, variants, variant, delta, mutations, voc
># Topic 18 Top Words:
># Highest Prob: covid, patients, disease, severe, immune, inflammatory, associated
># FREX: autoantibodies, ipf, balf, neutrophils, il, fibrosis, autoantibody
># Lift: 18f, 24hr, a2ar, aab, actinobacteria, activin, adiponectin
># Score: inflammatory, il, patients, cytokine, inflammation, cytokines, endothelial
># Topic 19 Top Words:
># Highest Prob: patients, treatment, covid, group, days, day, trial
># FREX: placebo, randomized, hcq, soc, azithromycin, arm, tocilizumab
># Lift: 200mg, 400mg, 500mg, 600mg, 800mg, aureobasidium, ayush
># Score: placebo, hcq, trial, patients, randomized, tocilizumab, hydroxychloroquine
># Topic 20 Top Words:
># Highest Prob: covid, symptoms, long, workers, infection, participants, study
># FREX: hcws, taste, smell, hcw, fatigue, workers, headache
># Lift: chemesthetic, dirty, eyewear, firefighters, ohs, principality, psychophysical
># Score: symptoms, hcws, workers, participants, symptom, hcw, fatigueTopic-document and term-topic distributions
As mentioned, the topic model is defined by the gamma
(distribution of topics over words) and beta (distribution of
terms over topics) matrices. With the help of the tidytext
package we can extract those into dataframes for a more detailed
analysis.
Each row in the following dataframe lists the probability
(gamma) of a given topic occurring in a given
document2.
># Rows: 593,820
># Columns: 3
># $ document
[3m
[38;5;246m<int>
[39m
[23m 1
[38;5;246m,
[39m2
[38;5;246m,
[39m3
[38;5;246m,
[39m4
[38;5;246m,
[39m5
[38;5;246m,
[39m6
[38;5;246m,
[39m7
[38;5;246m,
[39m8
[38;5;246m,
[39m9
[38;5;246m,
[39m10
[38;5;246m,
[39m11
[38;5;246m,
[39m12
[38;5;246m,
[39m13
[38;5;246m,
[39m14
[38;5;246m,
[39m15
[38;5;246m,
[39m16
[38;5;246m,
[39m17
[38;5;246m,
[39m18…
># $ topic
[3m
[38;5;246m<int>
[39m
[23m 1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1
[38;5;246m,
[39m1…
># $ gamma
[3m
[38;5;246m<dbl>
[39m
[23m 0.001283143
[38;5;246m,
[39m0.003547179
[38;5;246m,
[39m0.003745510
[38;5;246m,
[39m0.005495076
[38;5;246m,
[39m0.0630417…Similarly, the beta matrix (extracted into a dataframe)
lists for each row the probability (beta) of a given
term occurring in a given topic.
># Rows: 741,860
># Columns: 3
># $ topic
[3m
[38;5;246m<int>
[39m
[23m 1
[38;5;246m,
[39m2
[38;5;246m,
[39m3
[38;5;246m,
[39m4
[38;5;246m,
[39m5
[38;5;246m,
[39m6
[38;5;246m,
[39m7
[38;5;246m,
[39m8
[38;5;246m,
[39m9
[38;5;246m,
[39m10
[38;5;246m,
[39m11
[38;5;246m,
[39m12
[38;5;246m,
[39m13
[38;5;246m,
[39m14
[38;5;246m,
[39m15
[38;5;246m,
[39m16
[38;5;246m,
[39m17
[38;5;246m,
[39m18
[38;5;246m,
[39m1…
># $ term
[3m
[38;5;246m<chr>
[39m
[23m "#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"
[38;5;246m,
[39m"#1"…
># $ beta
[3m
[38;5;246m<dbl>
[39m
[23m 1.874650e-112
[38;5;246m,
[39m2.519509e-119
[38;5;246m,
[39m1.248970e-118
[38;5;246m,
[39m6.665355e-160
[38;5;246m,
[39m3.21…Starting with the beta matrix we can create word clouds to explore useful semantic labels for each topic.
Understanding and labeling topics
As mentioned previously stm can use the word
probabilities to also compute FREX (frequent and exclusive)
words per topic. Those can be retrieved with the
stm::labelTopics() function, which returns ordered lists of
the different word sets characterizing a topic. We can combine the
information about high probability and FREX words to
create word clouds for each topic, which might help to assign a summary
label for each topic.
library(tidyr)
library(tibble)
# get the top FREX words
frex_top20 <- as.data.frame(labelTopics(covid_model_K20, n = 20)$frex) %>%
rownames_to_column(var = "topic") %>%
pivot_longer(starts_with("V"), values_to = "term") %>%
mutate(is_frex = 1) %>%
select(-name)
topic_words <- tidytext::tidy(covid_model_K20, matrix = "beta") %>%
#filter(!(term %in% c("sars", "cov", "covid"))) %>%
mutate(topic = as.character(topic)) %>%
group_by(topic) %>%
arrange(-beta) %>%
slice_head(n = 50) %>%
mutate(beta_norm = (beta - min(beta)) / (max(beta) - min(beta))) %>%
ungroup() %>%
left_join(frex_top20, by = c("topic", "term")) %>%
mutate(is_frex = ifelse(is.na(is_frex), "0", "1")) %>%
filter(!(term %in% c("sars", "cov", "covid"))) %>%
mutate(topic = paste("Topic", topic))
ggplot(topic_words, aes(label = term, size = beta_norm, color = is_frex)) +
ggwordcloud::geom_text_wordcloud_area(shape = "square",
family = "Arial",
rm_outside = TRUE) +
scale_radius(range = c(4, 15)) +
scale_color_manual(values = c("0" = "black", "1" = "#D55E00")) +
facet_wrap(~topic, ncol = 4)
Word clouds showing the 50 most probable terms per topic. Words are scaled by normalized probability per topic, terms that are also among the top 20 FREX terms are highlighted in orange. For readability, the terms ‘sars’, ‘cov’ and ‘covid’ have been removed as they occur with high probability in most topics.
From the review of this combination of FREX and high probability terms distinct topics are emerging, such as: “epidemic models” (Topic 1), “vaccines” (Topic 11), “testing” (Topic 16), “virus variants” (Topic 17), “treatments” (Topic 2), “mortality risks” (Topic 3), “mental health” (Topic 4), “country-wise case reports” (Topic 8), “virus molecular structure” (Topic 9).
Deriving semantic labels will also benefit from a review of sample
documents that are representative for certain topics, i.e. where the
topic model assigns a high topic proportion of a given topic to the
document. Furthermore, the consistency of assigned semantic labels
should be checked prior to further analysis (see for example the oolong
package for a systematic approach).
Covariate effects
A key feature of STM is the incorporation of document covariates into
the topic model. In our example we considered the publication year and
the preprint server as covariates that might influence the prevalence of
a topic in a document. We also applied a regression for these covariates
to the fit model, which were extracted with
stm::estimateEffect().
As a first exploration we can print the regression tables for all or selected topics.
summary(covid_effect_K20)
>#
># Call:
># estimateEffect(formula = 1:20 ~ server * s(year), stmobj = covid_model_K20,
># metadata = covid_stm_docs$meta)
>#
>#
># Topic 1:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0234449 0.0036920 6.350 2.18e-10 ***
># servermedrxiv 0.1239685 0.0043664 28.392 < 2e-16 ***
># s(year)1 -0.0020910 0.0145920 -0.143 0.886
># s(year)2 -0.0004224 0.0152862 -0.028 0.978
># s(year)3 0.0083459 0.0067173 1.242 0.214
># servermedrxiv:s(year)1 -0.0694235 0.0172904 -4.015 5.96e-05 ***
># servermedrxiv:s(year)2 -0.0968399 0.0189587 -5.108 3.28e-07 ***
># servermedrxiv:s(year)3 -0.0705738 0.0085752 -8.230 < 2e-16 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 2:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.098647 0.002875 34.317 < 2e-16 ***
># servermedrxiv -0.085906 0.003158 -27.202 < 2e-16 ***
># s(year)1 -0.020088 0.011420 -1.759 0.0786 .
># s(year)2 -0.011175 0.013462 -0.830 0.4065
># s(year)3 -0.028594 0.004526 -6.318 2.69e-10 ***
># servermedrxiv:s(year)1 0.021779 0.012678 1.718 0.0858 .
># servermedrxiv:s(year)2 0.008106 0.014518 0.558 0.5766
># servermedrxiv:s(year)3 0.030343 0.005266 5.762 8.37e-09 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 3:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.011077 0.002809 3.943 8.08e-05 ***
># servermedrxiv 0.053945 0.003304 16.327 < 2e-16 ***
># s(year)1 -0.002580 0.011404 -0.226 0.821
># s(year)2 -0.002023 0.011888 -0.170 0.865
># s(year)3 0.001222 0.004902 0.249 0.803
># servermedrxiv:s(year)1 0.015001 0.013884 1.080 0.280
># servermedrxiv:s(year)2 0.014291 0.015040 0.950 0.342
># servermedrxiv:s(year)3 0.029818 0.006562 4.544 5.54e-06 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 4:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.008847 0.003281 2.696 0.00702 **
># servermedrxiv 0.057761 0.003844 15.028 < 2e-16 ***
># s(year)1 0.004904 0.013367 0.367 0.71372
># s(year)2 -0.005263 0.013979 -0.376 0.70656
># s(year)3 0.004530 0.005814 0.779 0.43584
># servermedrxiv:s(year)1 0.020017 0.016348 1.224 0.22078
># servermedrxiv:s(year)2 -0.003802 0.017986 -0.211 0.83258
># servermedrxiv:s(year)3 0.020303 0.007644 2.656 0.00791 **
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 5:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0049084 0.0023124 2.123 0.0338 *
># servermedrxiv 0.0543009 0.0027565 19.699 <2e-16 ***
># s(year)1 -0.0007422 0.0094201 -0.079 0.9372
># s(year)2 0.0001310 0.0098815 0.013 0.9894
># s(year)3 0.0011181 0.0041599 0.269 0.7881
># servermedrxiv:s(year)1 0.0260117 0.0114839 2.265 0.0235 *
># servermedrxiv:s(year)2 -0.0017266 0.0120470 -0.143 0.8860
># servermedrxiv:s(year)3 -0.0085558 0.0053387 -1.603 0.1090
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 6:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.148622 0.003823 38.880 < 2e-16 ***
># servermedrxiv -0.115802 0.004189 -27.644 < 2e-16 ***
># s(year)1 -0.021038 0.013367 -1.574 0.1155
># s(year)2 -0.049342 0.012242 -4.031 5.58e-05 ***
># s(year)3 -0.043739 0.006045 -7.235 4.76e-13 ***
># servermedrxiv:s(year)1 0.005868 0.014890 0.394 0.6935
># servermedrxiv:s(year)2 0.035529 0.013942 2.548 0.0108 *
># servermedrxiv:s(year)3 0.031022 0.006966 4.453 8.50e-06 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 7:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.010661 0.002895 3.682 0.000231 ***
># servermedrxiv 0.086061 0.003444 24.989 < 2e-16 ***
># s(year)1 -0.003943 0.011603 -0.340 0.734027
># s(year)2 -0.007907 0.012165 -0.650 0.515703
># s(year)3 -0.003280 0.005141 -0.638 0.523432
># servermedrxiv:s(year)1 -0.049427 0.014367 -3.440 0.000582 ***
># servermedrxiv:s(year)2 -0.019463 0.015633 -1.245 0.213120
># servermedrxiv:s(year)3 -0.038588 0.006495 -5.941 2.86e-09 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 8:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.027259 0.003009 9.058 < 2e-16 ***
># servermedrxiv 0.108699 0.003536 30.739 < 2e-16 ***
># s(year)1 -0.012945 0.012138 -1.066 0.28621
># s(year)2 -0.018204 0.012476 -1.459 0.14455
># s(year)3 -0.017085 0.005252 -3.253 0.00114 **
># servermedrxiv:s(year)1 -0.088375 0.014650 -6.032 1.64e-09 ***
># servermedrxiv:s(year)2 -0.046710 0.015507 -3.012 0.00260 **
># servermedrxiv:s(year)3 -0.060176 0.006582 -9.142 < 2e-16 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 9:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.170536 0.004559 37.408 < 2e-16 ***
># servermedrxiv -0.163997 0.004759 -34.462 < 2e-16 ***
># s(year)1 0.017930 0.018942 0.947 0.3439
># s(year)2 -0.034783 0.020036 -1.736 0.0826 .
># s(year)3 -0.028006 0.006534 -4.286 1.82e-05 ***
># servermedrxiv:s(year)1 -0.010857 0.020026 -0.542 0.5877
># servermedrxiv:s(year)2 0.033191 0.021831 1.520 0.1284
># servermedrxiv:s(year)3 0.027910 0.007115 3.923 8.77e-05 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 10:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0447172 0.0033446 13.370 < 2e-16 ***
># servermedrxiv 0.0118019 0.0036545 3.229 0.00124 **
># s(year)1 0.0043187 0.0147319 0.293 0.76941
># s(year)2 0.0002224 0.0133051 0.017 0.98666
># s(year)3 0.0262228 0.0058475 4.484 7.34e-06 ***
># servermedrxiv:s(year)1 -0.0178913 0.0164581 -1.087 0.27701
># servermedrxiv:s(year)2 -0.0303929 0.0155716 -1.952 0.05097 .
># servermedrxiv:s(year)3 -0.0409950 0.0070742 -5.795 6.90e-09 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 11:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.005054 0.002338 2.162 0.03066 *
># servermedrxiv 0.008234 0.002686 3.066 0.00217 **
># s(year)1 0.010111 0.009858 1.026 0.30505
># s(year)2 0.010014 0.010120 0.990 0.32240
># s(year)3 0.006443 0.004282 1.505 0.13238
># servermedrxiv:s(year)1 0.056938 0.011637 4.893 9.99e-07 ***
># servermedrxiv:s(year)2 0.073684 0.012418 5.934 2.99e-09 ***
># servermedrxiv:s(year)3 0.046987 0.005958 7.886 3.23e-15 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 12:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.021980 0.002843 7.731 1.10e-14 ***
># servermedrxiv 0.041958 0.003407 12.314 < 2e-16 ***
># s(year)1 -0.008566 0.011366 -0.754 0.451
># s(year)2 -0.005784 0.012314 -0.470 0.639
># s(year)3 -0.002578 0.004877 -0.529 0.597
># servermedrxiv:s(year)1 -0.006652 0.014012 -0.475 0.635
># servermedrxiv:s(year)2 0.009497 0.015954 0.595 0.552
># servermedrxiv:s(year)3 0.034803 0.006470 5.379 7.55e-08 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 13:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.104435 0.003148 33.180 < 2e-16 ***
># servermedrxiv -0.086419 0.003440 -25.123 < 2e-16 ***
># s(year)1 -0.031774 0.015613 -2.035 0.0419 *
># s(year)2 -0.027972 0.015196 -1.841 0.0657 .
># s(year)3 -0.023382 0.005570 -4.198 2.70e-05 ***
># servermedrxiv:s(year)1 0.077963 0.018016 4.327 1.51e-05 ***
># servermedrxiv:s(year)2 0.028315 0.017009 1.665 0.0960 .
># servermedrxiv:s(year)3 0.041166 0.006570 6.266 3.76e-10 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 14:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.124625 0.003458 36.041 <2e-16 ***
># servermedrxiv -0.112753 0.003650 -30.887 <2e-16 ***
># s(year)1 0.004254 0.015133 0.281 0.779
># s(year)2 -0.007602 0.015631 -0.486 0.627
># s(year)3 0.007547 0.006365 1.186 0.236
># servermedrxiv:s(year)1 -0.003166 0.016058 -0.197 0.844
># servermedrxiv:s(year)2 0.005410 0.016769 0.323 0.747
># servermedrxiv:s(year)3 -0.005132 0.007174 -0.715 0.474
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 15:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.053518 0.003476 15.398 < 2e-16 ***
># servermedrxiv -0.027449 0.003902 -7.035 2.03e-12 ***
># s(year)1 0.044388 0.014863 2.986 0.00283 **
># s(year)2 0.025077 0.015488 1.619 0.10544
># s(year)3 0.017662 0.006493 2.720 0.00652 **
># servermedrxiv:s(year)1 0.003930 0.017384 0.226 0.82115
># servermedrxiv:s(year)2 0.033328 0.017772 1.875 0.06076 .
># servermedrxiv:s(year)3 -0.003307 0.007628 -0.433 0.66466
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 16:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0531676 0.0034466 15.426 < 2e-16 ***
># servermedrxiv 0.0176342 0.0040500 4.354 1.34e-05 ***
># s(year)1 -0.0493934 0.0129310 -3.820 0.000134 ***
># s(year)2 -0.0296889 0.0131932 -2.250 0.024436 *
># s(year)3 -0.0278944 0.0057571 -4.845 1.27e-06 ***
># servermedrxiv:s(year)1 0.0607324 0.0150374 4.039 5.39e-05 ***
># servermedrxiv:s(year)2 0.0014310 0.0160447 0.089 0.928934
># servermedrxiv:s(year)3 -0.0006634 0.0072656 -0.091 0.927245
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 17:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.012514 0.002371 5.278 1.31e-07 ***
># servermedrxiv -0.010843 0.002658 -4.079 4.54e-05 ***
># s(year)1 0.062388 0.011969 5.212 1.88e-07 ***
># s(year)2 0.167403 0.015884 10.539 < 2e-16 ***
># s(year)3 0.074739 0.005334 14.012 < 2e-16 ***
># servermedrxiv:s(year)1 -0.043502 0.013431 -3.239 0.0012 **
># servermedrxiv:s(year)2 -0.080333 0.017880 -4.493 7.05e-06 ***
># servermedrxiv:s(year)3 -0.037898 0.006093 -6.220 5.03e-10 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 18:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.055678 0.003065 18.163 < 2e-16 ***
># servermedrxiv -0.019896 0.003490 -5.701 1.20e-08 ***
># s(year)1 0.005633 0.011949 0.471 0.63732
># s(year)2 0.004103 0.012755 0.322 0.74767
># s(year)3 0.029434 0.005806 5.070 4.01e-07 ***
># servermedrxiv:s(year)1 -0.004577 0.014626 -0.313 0.75432
># servermedrxiv:s(year)2 -0.004592 0.015460 -0.297 0.76646
># servermedrxiv:s(year)3 -0.020692 0.006932 -2.985 0.00284 **
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 19:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0135345 0.0025715 5.263 1.43e-07 ***
># servermedrxiv 0.0273562 0.0029744 9.197 < 2e-16 ***
># s(year)1 0.0005924 0.0100014 0.059 0.9528
># s(year)2 -0.0063075 0.0104244 -0.605 0.5451
># s(year)3 -0.0030460 0.0043654 -0.698 0.4853
># servermedrxiv:s(year)1 -0.0031975 0.0122727 -0.261 0.7945
># servermedrxiv:s(year)2 0.0282285 0.0130551 2.162 0.0306 *
># servermedrxiv:s(year)3 0.0027456 0.0057006 0.482 0.6301
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
>#
># Topic 20:
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 0.0067704 0.0019407 3.489 0.000486 ***
># servermedrxiv 0.0313936 0.0022126 14.189 < 2e-16 ***
># s(year)1 -0.0015268 0.0076128 -0.201 0.841049
># s(year)2 -0.0002849 0.0081858 -0.035 0.972241
># s(year)3 0.0004597 0.0034188 0.134 0.893036
># servermedrxiv:s(year)1 0.0087563 0.0095568 0.916 0.359554
># servermedrxiv:s(year)2 0.0130479 0.0106322 1.227 0.219756
># servermedrxiv:s(year)3 0.0213364 0.0046178 4.620 3.84e-06 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The following sections illustrate the exploration of covariate effects with several examples.
Publication year
stm offers several built-in methods to explore the
covariate effects visually. The visualization of the effect of preprint
publication year on expected topic proportions confirm
some trends that appear to match intuitively with different phases of
the Covid-19 pandemic. Modelling the spread of infections (Topic 1) had
a high relevance initially, but declined later, the same applies to
(PCR-)testing (Topic 16), while vaccines (Topic 11) received limited
coverage in earlier preprints, but became more important in later years
as they were developed and distributed.
plot(covid_effect_K20,
covariate = "year",
method = "continuous",
model = covid_model_K20,
topics = c(1, 16, 11),
xaxt = "n",
main = 'Effect of publication year on prevalence of Topic 1 ("epidemic \nmodels"), Topic 16 ("testing") and Topic 11 ("vaccines")',
labeltype = "prob",
xlab = "Publication year")
axis(1, at = c("2020","2021","2022","2023"), labels = c(2020, 2021, 2022, 2023))
Alternatively, the stminsights package could be used to
extract the same regression information from the stm
effects object and create customized charts.
library(stminsights)
year_effect <- get_effects(estimates = covid_effect_K20,
variable = "year",
type = "continuous")
year_effect %>%
mutate(topic = as.character(topic)) %>%
mutate(topic = paste("Topic", topic)) %>%
ggplot(aes(x = value, y = proportion)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2, linetype = 0) +
xlab("Publication year") +
ylab("Topic prevalence") +
facet_wrap(~topic, ncol = 4) +
theme_minimal()
Preprint server
The model also incorporated the preprint server (bioRxiv and
medRxiv) as a covariate. We hypothesized that the prevalence of
certain topics will be influenced by this covariate, i.e. that we are
more likely to see certain topics on either bioRxiv and
medRxiv. STM interprets this as a treatment (see (Roberts
et al. 2014) for background and examples) and we can extract
the effect of this treatment from the regression object returned by
estimateEffect().
The stm::plot() offers different methods to explore
those effects. The figure below lists the effect of ‘treatment
medrxiv’ for all topics in the model. The treatment can have a
positive or negative effect. Here we can for example see that Topic 4
(“mental health”) can be expected with higher prevalence in
medRxiv preprints, whereas Topic 9 (“virus molecular
structure”) will have a much lower prevalence in medRxiv
preprints, i.e. should be expected with higher prevalence in
bioRxiv preprints.
plot(covid_effect_K20,
covariate = "server",
#topics = c(9, 10, 16, 1),
model = covid_model_K20,
method = "difference",
cov.value1 = "medrxiv", cov.value2 = "biorxiv",
xlab = "higher biorxiv prevalence ... higher medrxiv prevalence",
xlim = c(-0.19, 0.1),
#labeltype = "prob",
main = "Effect of preprint server ('treatment medrxiv')")
We can also compare the topical prevalence effect of covariate values for a single topic. The figure below confirms the previous observation for Topic 9 (“virus molecular structure”) which has a prevalence close to zero in medRxiv preprints but a prevalence of approximately 0.17 for bioRxiv preprints, i.e. it seems to appear almost exclusively in the latter.
plot(covid_effect_K20,
covariate = "server",
topics = c(9),
model = covid_model_K20,
method = "pointestimate",
xlab = "Topical prevalence",
xlim = c(-0.04, 0.2),
#labeltype = "prob",
main = "Effect of preprint server covariate for Topic 9")
Combination of preprint server and publication year
Finally, we can also explore the combined effect of covariates; we
assumed an interacting effect of the two covariates server
and year. The plot below visualizes this for Topic 17
(“virus variants”) and illustrates that the topical prevalence increased
throughout the pandemic in all preprints but that the topic apparently
received more coverage in bioRxiv preprints.
Again, stminsights can be used to retrieve this
information and create customized visualizations for the combined
covariate effect.
biorxiv_effect <- get_effects(covid_effect_K20,
variable = "year", type = "continuous",
moderator = "server", modval = "biorxiv")
medrxiv_effect <- get_effects(covid_effect_K20,
variable = "year", type = "continuous",
moderator = "server", modval = "medrxiv")
server_effects <- bind_rows(biorxiv_effect, medrxiv_effect)
server_effects %>%
mutate(topic = as.character(topic)) %>%
mutate(topic = paste("Topic", topic)) %>%
ggplot(aes(x = value, y = proportion, color = moderator,
group = moderator, fill = moderator)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper,
fill = moderator), alpha = 0.2, linetype = 0) +
xlab("Publication year") +
ylab("Topic prevalence") +
facet_wrap(~topic, ncol = 4) +
theme_minimal() +
theme(legend.position = "bottom") 
Exploring the topic structure
Topic correlations
Finally, considering that STM is extending (Roberts, Stewart, and
Tingley 2019) a correlated topic model (Blei and
Lafferty 2007) we may also explore if topics are frequently
cooccuring. A a topic correlation matrix is retrieved with
stm::topicCorr() which can then be plotted. The resulting
figure below indicates at least two clusters of topics.
The correlation matrix can be used for further analysis (e.g. networks and clusters) or alternative visualizations. Both are outside the scope of this example and the course module. For illustration, the following figure creates an alternative network visualization, with nodes colored according to a community analysis of the topical network, which here reveals four distinct topical clusters.