# CODY: A graph-based framework for the analysis of COnversation DYnamics in
online social networks ††thanks: We thank the Klaus Tschira Foundation for
funding this research in the framework of the EPINetz project:
https://epinetz.de.

1st John Ziegler  Institute of Computer Science
Heidelberg University
Heidelberg, Germany
<EMAIL_ADDRESS>2nd Fabian Kneissl  Institute of
Computer Science
Heidelberg University
Heidelberg, Germany
<EMAIL_ADDRESS>3rd Michael Gertz  Institute of
Computer Science
Heidelberg University
Heidelberg, Germany
<EMAIL_ADDRESS>
###### Abstract

Conversations are an integral part of online social media, and gaining
insights into these conversations is of significant value for many commercial
as well as academic use cases. From a computational perspective, however,
analyzing conversation data is complex, and numerous aspects must be
considered. Next to the structure of conversations, the discussed content – as
well as their dynamics – have to be taken into account. Still, most existing
modelling and analysis approaches focus only on one of these aspects and, in
particular, lack the capability to investigate the temporal evolution of a
conversation. To address these shortcomings, in this work, we present CODY, a
content-aware, graph-based framework to study the dynamics of online
conversations along multiple dimensions. Its capabilities are extensively
demonstrated by conducting three experiments based on a large conversation
dataset from the German political Twittersphere. First, the posting activity
across the lifetime of conversations is examined. We find that posting
activity follows an exponential saturation pattern. Based on this activity
model, we develop a volume-based sampling method to study conversation
dynamics using temporal network snapshots. In a second experiment, we focus on
the evolution of a conversation’s structure and leverage a novel metric, the
temporal Wiener index, for that. Results indicate that as conversations
progress, a conversation’s structure tends to be less sprawling and more
centered around the original seed post. Furthermore, focusing on the dynamics
of content in conversations, the evolution of hashtag usage within
conversations is studied. Initially used hashtags do not necessarily keep
their dominant prevalence throughout the lifetime of a conversation. Instead,
various ”hashtag hijacking” scenarios are found. Also, during the conduction
of this work, we developed a software package to facilitate the analysis of
heterogeneous information networks. It is made publicly available to be used
by the research community.

###### Index Terms:

social media analytics, temporal heterogeneous information networks,
conversation analysis, Twitter data

Figure 1: The graph-based CODY model incorporates three central dimensions of
a conversation: structure, content and dynamics. Different node and edge types
are represented by different symbols and line styles, respectively. Temporal
snapshots enable the dynamics of the conversation to be shown.

##  I Introduction

Interpersonal communication is integral to many social media platforms and is
well-known through conversational features such as Twitter replies or YouTube
comments. Extracting valuable information from these conversations can be
beneficial for many commercial as well as academic use cases. For this,
elaborate modelling and analysis approaches are required. Still, existing
approaches show several deficiencies. In a recent work by Brambilla et al.
[1], the authors state that ”[m]ost studies on social networks have focused
only on user relationships _or_ [emphasis added] on the shared content, while
ignoring the valuable information hidden in the digital conversations, in
terms of structure of the discussion and relation between contents, which is
essential for understanding online communication behavior.” They highlight the
need for an analysis model to incorporate a conversation’s content _and_
structure. While this insight leads towards a more holistic conversation
analysis model, it still lacks an essential aspect. As conversations are
fundamentally characterized by their temporal evolution, a respective analysis
model should also consider a conversation’s dynamics. To the best of our
knowledge, existing conversation analysis approaches are still missing this
aspect [1, 2, 3]. We overcome these shortcomings by proposing a graph-based
conversation analysis model, named CODY, which incorporates the conversation’s
_structure_ , _content_ and _dynamics_. Its capabilities are demonstrated by
focusing on three main research questions:

  1. RQ1

How does the posting activity change over the lifetime of a conversation? Does
a generalizing model describe the observed activity pattern?

  2. RQ2

How do conversations structurally evolve? Is there a metric to quantitatively
describe the structural evolution?

  3. RQ3

How does the focus of discussed topics change over the lifetime of a
conversation? Do initially used hashtags keep their dominant prevalence?

To answer these research questions, an analysis model needs to be developed
that is versatile enough to cope with the three dimensions of a conversation –
structure, content and dynamics – as no current approach exists for that. The
proposed CODY model is based on the concept of temporal heterogeneous
information networks (THINs) (cf. [4, 5, 6]) and therefore comes with proper
flexibility to model a conversation’s evolution, along with its structure and
content. By leveraging a large conversation dataset from the German political
Twittersphere, the model’s applicability is shown, and proposed research
questions are investigated. The developed methodology and the results of our
experiments make up several contributions:

  * •

The graph-based CODY model constitutes a versatile framework to study the
structure, content and dynamics of conversations. Its flexibility in terms of
integrated node/edge types allows to model various (content) objects and their
interplay within a conversation.

  * •

RQ1: Based on empirical results, the posting activity of a conversation
follows an exponential saturation pattern. Contributions occur mainly at the
early stage of a conversation, and participation becomes less towards the end.

  * •

The developed activity model lays the theoretical foundation for a volume-
based sampling approach used to model a conversation’s evolution based on
temporal snapshots. The common issue of finding an appropriate sampling
technique is overcome with that.

  * •

RQ2: The newly proposed temporal Wiener index metric allows to quantitatively
measure a conversation’s structural evolution. Experimental results show that
as conversations progress, their structure tends to be less sprawling and more
centered around the original post.

  * •

RQ3: Focused on the content of conversations, investigations regarding the
temporal prevalence of hashtags show that initially used hashtags do not
necessarily keep their dominant role across a conversation’s lifetime.
Instead, various ”hashtag hijacking” scenarios can be observed.

  * •

Along with the conducted experiments, we developed a software package to
facilitate the analysis of heterogeneous information networks. Its source
code111GitHub: https://github.com/codingfabi/hetpy is publicly available to be
used by the research community.

Regarding the used terminology, in line with previous work [3], we focus on
conversations that start from an initial social media post and that are made
up of direct user interactions related to this post rather than
”conversations” that are centered around an event or other social media
entities, such as @mentions. In our understanding, conversations must be based
on direct interpersonal communication instead of latent social interactions.
In the context of the Twitter platform, replies are considered ”[…] the most
direct communication sign between two users” [3]. Therefore, Twitter reply
threads, as leveraged for the conducted experiments, make up conversations.

##  II Related Work

Most closely related to our work are studies about the modelling and analysis
of online conversations. As such, Brambilla et al. [1] conduct an in-depth
study of various online conversations. For their analysis, they propose a
graph-based framework and mainly focus on intention analysis as well as
network construction. In contrast to us, they do not consider the dynamics of
the conversation graph itself. Still, they consider temporal aspects of the
analyzed conversations, such as their duration and typical reply times.
Similar to our approach, they also view the network as heterogeneous by
including different node and edge types. Recently, Botzer and Weninger [7]
leverage entity graphs extracted from online conversation data, in their case
Reddit, to study whether online discourse is predictable, how the online
conversations are structured and whether theories on spreading activation can
be successfully applied in the context of online discourse analysis. Compared
to our work, they do not investigate the conversation graphs’ dynamics or
consider the heterogeneity of the conversations’ content. Further, in one of
the early works on conversation graph analysis, Cogan et al. [3] investigate
the characteristics of conversation networks emerging from an initial social
media seed post. This understanding of a conversation graph is in line with
our definition. In their work, they propose a method to gather conversation
data as completely as possible and conduct an experimental analysis based on
Twitter data. Their outlook highlights the importance of future work
considering a conversation graph’s temporal evolution. Focused on the aspect
of toxicity, Saveski et al. [2] also study the structure of Twitter
conversations. They conduct their analysis on different levels, from
individual users to the entire conversation group. Even though applied in a
static setting, they also leverage the Wiener index to investigate the
topology of analyzed reply trees. In two prediction tasks, they aim at
forecasting future toxic conversation characteristics based on previous
conversation data. In another application-focused work, Mathew et al. [8]
analyze Twitter conversation data related to e-commerce promotional events.
They model the data as a user-centric network by considering replies,
mentions, and retweets for the weighted edges and Twitter users as nodes.
Regarding the network’s temporal evolution, they determine the users’ flow
between different network components across time. Unfortunately, they do not
elaborate on their temporal graph model but only propose a sampling approach
using fixed time windows. Also, they do not investigate the dynamics of any
network metric and solely focus on users, not the content, as part of a
conversation.

Information diffusion is a topic that is also extensively studied in the
context of online social networks. The survey by Guille et at. [9] and the
more recent one by Yujie [10] give a good overview of this research field.
Commonly used information diffusion models also leverage topological
information of the underlying social network. The spreading process is
frequently modelled as a network, referred to as diffusion graph. In contrast
to our work, information diffusion mainly investigates how information spreads
across the entire network. We, however, are less interested in the information
diffusion process but instead focus on the graph-centric evolution of
individual online conversations, which includes not only the discussed
information but also the participating actors and the structure of a
conversation. Still, based on our model, temporal meta paths related to the
concept of THINs could be leveraged to study information diffusion processes
within individual conversations.

As part of our methodology, we leverage THINs to model the evolution of online
conversations. Even though, to the best of our knowledge, none of the existing
work applied THINs in the context of conversation analysis, they are
extensively used for other use cases. For example, Sun et al. [4] propose a
method to detect communities in THINs. They understand THINs as a sequence of
temporal network snapshots, which aligns with our approach. Similarly,
Cuzzocrea and Folino [11] deal with community detection in temporal
information networks. They specifically focus on the temporal tracking of the
detected communities and propose a detection method to analyze structural
changes regarding the community structure of the evolving network. Further,
again in the context of THINs, Fard et al. [6] tackle the task of meta path
prediction. They leverage temporal snapshots to model a network’s evolution.
Sajadmanesh et al. [12] also predict relationships in THINs. Additionally,
they predict the time when these relationships will occur. Further, multiple
other use cases are elaborated in the context of THINs. These include the
engineering of node embeddings [13] [14], sequential recommendation [15],
network motifs [5], and link inference [16] [17].

##  III CODY Model

As already described, a conversation’s structure, content and dynamics must be
considered for a conversation analysis model to be leverageable in different
analysis scenarios. Following the same pattern, this section introduces the
CODY model in three steps. First, Section III-A explains how the CODY model
considers the structure of conversations. Secondly, the temporal evolution of
a conversation is integrated into the model in Section III-B. Finally, Section
III-C complements the model by additionally considering the content of
conversations. Apart from the conversation model, metrics to investigate a
conversation’s structural evolution are needed. Therefore, Section III-D
introduces a novel metric, denoted as _temporal_ Wiener index, to
quantitatively study how a conversation’s structure evolves over time.

###  III-A Structural conversation model

Considering the structure of a conversation is an essential requirement for a
holistic analysis model. On an abstract level, an online conversation consists
of posts (e.g., tweets or comments) and links between these posts that
represent semantic relationships (e.g., replying or commenting). The concept
of graphs222Both terms, _graph_ and _network_ , are used interchangeably in
this work. The same is true for the terms ”object”, ”vertex”, and ”node”, as
well as ”edge”, ”link”, and ”relationship”. is particularly suitable for
modelling objects, in our case posts, and their relationships. We resort to
these graphs as the primary modelling approach for the CODY framework. Posts
contained in the conversation are represented as nodes and are linked via
semantic relationships represented as edges. For a conversation graph with
nodes $V$ and edges $E$ to take the different semantic types into account, the
network definition comes with a node type mapping $\phi:V\rightarrow A$ and a
link type mapping $\psi:E\rightarrow R$ with $A$ and $R$ denoting the
network’s node and edge types respectively. This network typing aligns with
the concept of HINs (cf. Sun and Han [18]).

######  Definition III.1 (Structural conversation network)

The structural network $G_{struc}=(V,E)$ of a conversation consists of a set
of posts modelled as nodes $V$ ($\phi(v)=\text{post}$ with $v\in V$) and a set
of links $E\subset V\times V$ among these posts. Each link $e=(v_{i},v_{j})$
is a tuple of the two adjacent posts $v_{i}$ and $v_{j}$.

A single conversation is treated as an individual network based on Definition
III.1. Generally, interpersonal communication modelled by the structural
conversation network implies directionality and a post is usually meant as a
response to another, e.g., reply $\rightarrow$ seed post. Therefore, the
structural conversation network is a _directed_ network. Also, the outgoing
degree of a node is one except for the root node, which is not a response to
any other post. Therefore, if one would not consider the edges’
directionality, two posts would be connected by exactly one path, and the
structural conversation network could be described as a _tree_ (cf. Cogan et
al. [3]).

###  III-B Temporal conversation model

Conversations are inherently dynamic and evolve over time. It is crucial to
consider these dynamics for modelling a conversation. Typically the posts
contained in a conversation come with a timestamp, such as the publishing
date. Formally, the assignment of timestamps is described by a node time
mapping $\tau:V\rightarrow T$ with $T$ representing the domain of time. Given
that the posts’ timestamps impose a natural order on a conversation’s
elements, the respective structural conversation network can be split into a
temporal sequence of snapshots.

######  Definition III.2 (Temporal conversation model)

Each post in a structural conversation network comes with a timestamp:
$\tau(v)=t$. Given a sequence of $n$ timestamps $(t_{i})_{i=0}^{n}$ with
$t_{i}<t_{i+1}$, the structural network $G_{struc}$ can be split into $n$
temporal snapshots: $G_{temp}=(G_{i})_{i=1}^{n}$. For each snapshot
$G_{i}=(V_{i},E_{i})\in G_{temp}$ the contained nodes and edges are defined as
follows: $V_{i}=\\{v\in V\;|\;t_{0}\leq\tau(v)=t<t_{i}\\}$ and
$E_{i}=\\{e=(v_{i},v_{j})\in E\;|\;v_{i},v_{j}\in V_{i}\\}$.

Notably, the time windows used for the snapshots do not have to be of equal
size. The condition $|t_{i+1}-t_{i}|=|t_{i+3}-t_{i+2}|$ does not necessarily
hold. The network $G_{struc}$ might also be sampled according to the volume of
conversation posts as an alternative to equally sized time windows. In such a
case, the condition $|V_{i}^{\prime}|=|\\{v\in
V\;|\;t_{i-1}\leq\tau(v)=t<t_{i}\\}|=|V_{i+1}^{\prime}|=|\\{v\in
V\;|\;t_{i}\leq\tau(v)=t<t_{i+1}\\}|$ holds, meaning that the amount of newly
added posts is constant across the network snapshots. This _volume-based
sampling_ approach is especially suitable in cases where conversation posts
are not equally distributed over time. Therefore, a time-based sampling, as
suggested by Cogan et al. [3] and used by Mathew et al. [8], would lead to a
highly skewed distribution of newly added posts across the network snapshots.

###  III-C Content-aware conversation model

Next to structure and dynamics, the actual content of a conversation must be
considered by a conversation model as well. Various types of semantically
meaningful content can be part of a conversation, such as written text, shared
URLs, posted images or involved users. These content objects are extracted
from the analyzed posts to be incorporated into the conversation network.
Figure 2 illustrates this extraction process. In the content-aware
conversation network, the extracted content objects are linked to their
belonging post node.

Figure 2: Modelled content objects are extracted from the analyzed posts.
Extracted content, the two hashtags and the image, is linked to the respective
post node in the content-aware conversation network.

######  Definition III.3 (Content-aware conversation network)

The content-aware conversation model is based on the network $G_{struc}$,
enriched by additional content objects modelled as nodes. Additional node
types $a\in A$ indicate the different content types. The actual content
objects are extracted from the respective conversation posts:
$\epsilon:V_{P}\rightarrow 2^{V_{C}}$ with $V_{P}=\\{v\in
V\;|\;\phi(v)=\text{post}\\}$ and $V_{C}=\\{v\in V\;|\;\phi(v)\in
A\backslash\\{\text{post}\\}\\}$. To indicate the extraction process, in the
content-aware conversation network $G_{cont}$ the content objects are linked
to the original post node(s): $v_{p}\rightarrow v_{c}$ with $v_{c}\in V_{C}$
and $v_{c}\in\epsilon(v_{p})$. These links might as well be typed. Respective
types are modelled as additional edge types $r\in R$.

The same content object might belong to different posts in the content-aware
conversation network (e.g., the same URL might be shared in multiple posts).
Therefore, two formerly disconnected post nodes might become indirectly linked
via shared content nodes (e.g., URL, hashtag, username).

The content enrichment of the structural conversation network also applies to
the individual snapshots of the temporal conversation network $G_{temp}$. Each
network snapshot then contains the content objects belonging to the contained
posts. Together, the temporal and content-aware conversation model
incorporates a conversation’s structure, dynamics and content. Formally, the
CODY model can be seen as a THIN.

###  III-D Temporal Wiener index

To quantitatively describe the structure of conversation networks, tailored
metrics are needed. In a recent work, the Wiener index, known from the field
of chemistry [19], has been leveraged to measure the structure of reply trees
[2]. According to Saveski et al. [2], the Wiener index describes two
characteristic topologies of a tree-like graph in its extreme cases. On the
one hand, a low Wiener index indicates a conversation in which all non-seed
posts respond to the original seed post, and these non-seed posts do not show
any interactions among each other. On the other hand, a high Wiener index
indicates a conversation tree with a single dominant interaction branch. To
the best of our knowledge, no work extends the Wiener index to fit temporally
evolving conversation trees. Based on the definition of the static Wiener
index employed by Saveski et al. [2], we define the _temporal_ Wiener index as
formalized in the following:

######  Definition III.4 (Temporal Wiener index)

The temporal Wiener index $w(G_{temp})$ of a temporal sequence
$(G_{i})_{i=1}^{n}$ of tree-like structured conversation networks is
determined by the average distance between post nodes in each network
snapshot:

$w(G_{temp})=\left(\frac{1}{|V_{i}|\cdot(|V_{i}|-1)}\sum_{m\in
V_{i}}\sum_{n\in V_{i}}\delta_{i}(m,n)\right)_{i=1}^{n}$ (1)

with $\delta_{i}(m,n)$ denoting the length of the shortest path between the
nodes $m$ and $n$ in the temporal snapshot $G_{i}$.

Given a temporal sequence of conversation networks, the temporal Wiener index
thus itself is given as a temporal sequence of static Wiener indices.
Therefore, the time series of Wiener indices allows for studying the temporal
changes of the Wiener index across the temporal evolution of a conversation.

##  IV Experiments

To demonstrate the capabilities of the CODY model, this section discusses the
results of various conducted experiments. For all experiments, a Twitter
conversation dataset that we collected is leveraged (cf. Section IV-A). Also,
our developed software package (cf. Section I) is used to analyze the
experiments’ results with regard to HIN analysis.

###  IV-A Twitter conversation dataset

The dataset used for the subsequent experiments is based on a collection of
German political Twitter accounts published as EPINetz Twitter dataset [20].
We collected conversations that started from a seed post of one of these
accounts by leveraging the official Twitter search API v2333Twitter Developer
Platform: https://developer.twitter.com/en/docs/twitter-
api/tweets/search/introduction, accessed 10/07/23. According to Maireder and
Ausserhofer, the political domain is particularly suitable for examining the
structure of conversations ([21], as cited in [22]): ”Retweets and
mention/reply functions have been found to influence the structure of
conversational discourse if it is especially related to politics.” In total,
the dataset consists of more than 2M tweets as part of 5k conversations, with
each conversation including at least 50 posts. Figure 3 shows the distribution
of conversation sizes. The tweets were published between December 2020 and
January 2023, and the median conversation duration in the dataset is about
seven days.

Figure 3: Distribution of conversation sizes as contained in the leveraged
dataset of the German political Twittersphere.

According to Sun and Han, the Twitter information network is well suited to be
modelled as HIN [18]. In line with the CODY model, we view each conversation
as a temporal and content-aware conversation network with the contained tweets
as posts and hashtags as content objects. Figure 4 illustrates the network
schema of the conversation model used for the subsequent experiments. In the
analyzed networks, links between tweets indicate replies which are seen as
”[…] the most direct communication sign between two users” [3] and therefore,
make up the structure of a conversation.

Regarding data quality, there is no guarantee that the complete lifetime of a
conversation is captured. This would require a continued re-crawl of the
conversation data, which is unfeasible due to API usage restrictions.
Nevertheless, the reported high average conversation duration suggests that
most of the conversation’s evolution is captured. This assumption is also
underlined by the conducted experiment, which examines the posting activity of
the collected conversations (cf. Section IV-B). Most activities seem to occur
within the first 10% of a conversation’s lifetime. Therefore, sufficient posts
should be contained in the dataset even for shorter conversations, for which
data of less than a day is captured. In their analysis, Cogan et al. state
that conversation activity usually stops after 6 hours [3]. Further,
additional data quality issues related to online conversations have been
reported in the past, such as replies to private profiles or the deletion of
posts [3]. However, less than 5% of conversations seem to be affected by that
[3]. Missing posts lead to a split of the conversation graph into multiple
components. Therefore, we take the weakly connected, giant component for
creating the reply trees to handle the cases in which already deleted posts or
those that link to private profiles are part of a conversation. Also, in a
prior step to generate these reply trees, tweets are grouped into
conversations based on the conversation_id field contained in the raw API
payload. Among the tweets of a conversation, reply links are determined by the
additional referenced_tweets field and the contained entries of type
replied_to.

Even though the following experiments have been conducted based on the
described dataset, we argue that the CODY framework is applicable to a wide
range of use cases and conversation data of platforms different from Twitter
could be used as well. The primary assumption of the model, which is to have
interrelated and timestamped conversation posts, is fulfilled by a large
variety of online conversation data sources.

Figure 4: Network schema used for the conducted experiments. A conversation’s
reply tree consist of tweet nodes connected by replies. Further, the ”usage”
relationship indicates that a hashtag is used within a tweet.

###  IV-B Posting activity

Analyzing the posting activity across the lifetime of a conversation is
crucial for a more general understanding of conversation dynamics (cf. RQ1).
It lays the foundation for deciding which temporal sampling approach must be
leveraged to model the conversation dynamics correctly (cf. Section III-B).
Figure 5 shows the volume completion rate (relative number of posts) of the
conversations in our dataset plotted against the conversations’ lifetime (time
completion rate). The dark and lighter value bands indicate mean values and
their standard deviations. Values are calculated using a $10^{-5}$ resolution
of the completion rate. The large fluctuations show that temporal posting
activity varies a lot between conversations. Still, a general activity pattern
across the lifetime of a conversation can be observed. If $\lambda_{volume}$
denotes the volume completion rate of a conversation, which can also be
described as the fraction of the total amount of conversation posts, and
$\lambda_{time}$ denotes the time completion rate (lifetime) of a
conversation, a conversation’s posting activity can be modelled as

$\lambda_{volume}=1-e^{-\alpha\lambda_{time}}$ (2)

with $\alpha$ as the saturation rate. As shown by the fitted function in
Figure 5 and indicated by a reduced chi-square value, $\chi^{2}_{\nu}$, of
$0.4$, this model, despite its simplicity, describes the process in a precise
manner. For the curve fitting, the lmfit Python package is used [23]. Most of
the posting activity happens during the early stage of a conversation. After
an exponential increase in conversation posts, the posting activity slows down
and converges towards its final saturation level. This saturating process is
also characterized by the saturation time constant $\Gamma$, defined as
$\Gamma=\frac{1}{\alpha}$. For our dataset, the saturation time is about
$0.03$, meaning that on average, after 3% of the conversation lifetime, about
$1-\frac{1}{e}\approx 0.63\;\widehat{=}\;63\%$ of the posts are published.

Figure 5: The evolution of a conversation regarding the number of posts is
characterized by an exponential saturation function.

The non-linear activity pattern must be considered when selecting the sampling
method for the temporal modelling of conversations (cf. Section III-B). Given
the exponential increase in posts during the early stage of a conversation, a
sampling approach using fixed time windows would lead to a highly skewed
distribution of newly added posts across the network snapshots of a
conversation. In the beginning, the network snapshots would grow
significantly, whereas in later stages, only minor changes would be captured.
Quantitatively the change rate of sampling snapshots can be described by the
derivative of the posting activity:

$\lambda_{volume}^{\prime}=\frac{d}{d\lambda_{time}}\lambda_{volume}=\alpha
e^{-\alpha\lambda_{time}}.$ (3)

Therefore, the ratio of snapshot change rates during two points in time
$\lambda_{time_{1}}$ and $\lambda_{time_{2}}$ of a conversation’s lifetime is
given as:

$\frac{\lambda_{volume}^{\prime}(\lambda_{time_{1}})}{\lambda_{volume}^{\prime}(\lambda_{time_{2}})}=e^{\alpha(\lambda_{time_{2}}-\lambda_{time_{1}})}.$
(4)

For our conversation activity model ($\alpha=32.5$), about $131$ times more
posts would be added to a snapshot right at the beginning of a conversation
($\lambda_{time_{1}}=0.05$) compared to a later stage
($\lambda_{time_{2}}=0.2$) by using fixed time windows. Therefore, to
correctly capture changes across the entire lifetime of a conversation, the
exponential saturation regarding the posting activity needs to be considered.
This can easily be done using fixed volume sampling instead of fixed time
windows. The conversation size difference between two successive snapshots is
kept constant for this. In our case, we determine a sampling rate of five new
posts per snapshot. This sampling approach is also used for the subsequent
experiments. Given that the relative volume of posted tweets can be seen as
_normalized_ conversation time, we use a conversation’s completion rate as the
unit for the temporal dimension whenever appropriate.

###  IV-C Temporal Wiener index

Figure 6: Exemplary temporal Wiener index. After a steep increase at the
beginning of the conversation, the value rapidly decreases. Structural changes
in the conversation also reflect this evolution. Over time most posts directly
link to the original seed post, and offshoots are rarely observed.

The CODY framework views conversations as networks. This approach allows us to
analyze the conversations’ topological structure. Leveraging the _temporal_
Wiener index, this can also be done by considering the conversation’s dynamics
(cf. RQ2). For the following analysis, we calculate the temporal Wiener
indices for the conversations contained in our dataset. Only their reply trees
are considered without additional content objects (i.e., hashtags) to comply
with the definition of the Wiener index metric. For the software
implementation, the igraph Python package is used [24]. To handle the edge
case of a single-node reply tree, we define the Wiener index as zero at this
conversation state.

Figure 6 shows a temporal Wiener index example. As one can observe, the
structural evolution of the network is also reflected by the temporal change
of the Wiener index. The initial rise of the Wiener index value is caused by
the first reply, which directly links to the original seed post. Given that
only one shortest path of length one is contained in the network at this
stage, namely from the first reply post to the original post, a maximum Wiener
index value of $0.5$ is derived for this network snapshot. From there on, most
posts are directly linked to the seed post, reflected by the rapid decrease of
the temporal Wiener index. Overall, off-shooting branches are rarely observed
in the conversation network, which explains its low final value of about
$0.03$.

To not just focus on a single conversation example but to analyze the
structural evolution of conversations in a more general way, the average
temporal Wiener index across the lifetime of a conversation is shown in Figure
7 in the form of a box plot. First, it should be noted that temporal Wiener
index values vary a lot between the different conversations, as shown by the
large quartile ranges, the high variability, and the outliers. Still, a
general evolution pattern can be observed. Over time, the Wiener index
decreases and converges towards a final median value of about $0.01$.
Structurally, this means that over time, most replies directly link to the
original seed post, and sprawling conversation branches do not seem to play a
significant role in the conversation network. Nevertheless, to return to the
high data variability, in rare cases, offshoots might be observed in the
conversation structure, reflected by the higher outlier values as shown in the
box plot. Interestingly, Cogan et al. report that more than 60% of their
analyzed Twitter conversation graphs follow the structure of paths [3]. Apart
from the described outliers, this structure is less prevalent in our case as
it would lead to significantly higher final Wiener index values. Many reasons
are conceivable for this structural difference. One might be that, in our
case, most conversation participants wanted to interact directly with the
political actor who posted the original seed post. Thus, more star-like
conversation structures are observed.

Figure 7: Box plot showing the averaged temporal evolution of Wiener index
values across the lifetime of the analyzed conversations. Decreasing values
indicate the tendency of conversations to become less sprawling.

###  IV-D Hashtag Hijacking

In a holistic conversation analysis framework, the discussed content and its
dynamics should also be considered next to the structure of the conversation.
On the Twitter platform, used hashtags are an indicator of the discussed
topics and the semantics of posted content. In that sense, hashtags can be
leveraged to analyze how the topical focus of a conversation changes over
time. More specifically, one can investigate whether initially used hashtags
keep their dominant prevalence over the lifetime of a conversation (cf. RQ3).
For this, we leverage the proposed Twitter conversation dataset and, given
that hashtags are represented as nodes in the content-aware conversation
network, temporally track their degree centrality as an indicator of the
importance they play in the conversation. In the second step, the temporal
degree evolution of the hashtag(s) used in the initial seed post is compared
to the evolution of the other hashtags used later in the conversation. Suppose
a non-initial hashtag reaches a higher degree centrality value than the
initially used hashtag(s) at the end of a conversation; we denote the scenario
as _hashtag hijacking_. In case of an overtake that is followed by a ”retake”
of an initial hashtag so that it reclaims the highest prevalence with regards
to its degree centrality in the conversation network at the conversation’s
ending, the scenario is called _failed hijacking_. It is not a hijacking
scenario if no non-initial hashtags are used in a conversation, or no overtake
happens. If multiple hashtags take over the initial hashtag(s), the scenario
is only counted as a single hijacking. Semantically, a hashtag hijacking
scenario can be seen as an indicator of a _topical shift_ within a
conversation. Nevertheless, it has to be noted that at this point, we do not
employ any natural language processing techniques to analyze the semantic
similarity of the prevalent hashtags in a conversation and whether these are
used synonymously. Takeovers might as well be performed by semantically
similar hashtags, and therefore, only a slight semantic shift would be caused.
An extension to overcome this shortcoming would be a valuable future
contribution.

To apply the methodology, we focus on conversations where the initial seed
post contains at least one hashtag. Also, conversations should have at least
one hashtag that is used a sufficient number of times. In our case, we set
this threshold to five usages in a conversation. About 1.6k of the 5k
conversations in the analyzed Twitter corpus fulfil the outlined requirements.
In total, about 1k hijacking scenarios ($\sim 67\%$) are detected in these
conversations. Further, 111 failed hijacking scenarios ($\sim 7\%$) and 427
scenarios without hijacking ($\sim 26\%$) are found. Compared to the total
number of conversations in the dataset, successful hashtag hijackings comprise
a considerable portion of about $22\%$.

Figure 8: Example of a hashtag hijacking scenario. The initial hashtag
_Klimaschutz_ (Eng. climate protection) is overtaken by the hashtag
_Tempolimit_ (Eng. speed limit) in terms of degree centrality.

As an example of a hijacking scenario, Figure 8 shows the temporal degree
centrality of three hashtags used in the related Twitter
conversation444Twitter conversation:
https://twitter.com/Wissing/status/1596101821994606592, accessed 10/07/23. The
original tweet of the conversation criticizes a climate protest by the ”Letzte
Generation” (Eng. ”Last Generation”) movement, in which activists glued
themselves to the tarmac of an airport in order to interrupt flight
operations. The _Klimaschutz_ (Eng. climate protection) hashtag is used in
this tweet. Interestingly, the initial thematic positioning of the tweet is
overtaken. This takeover is reflected by the hashtag hijacking of the
_Tempolimit_ (Eng. speed limit) hashtag, which represents the discussion about
a speed limit on German autobahns. In this case, the hijacking can be
understood as a counter-positioning of many conversation participants.
Alternatively, their intention might have been to relate this conversation to
another climate-related conversation.

Apart from a hashtag hijacking scenario occurring within a conversation, some
use cases might also benefit from an analysis of _when_ the takeovers happen.
Figure 9 shows a heatmap of the occurrence probability of the takeovers
related to the (failed) hijacking scenarios across the lifetime of a
conversation. Dark patches indicate a high takeover occurrence probability,
whereas bright ones indicate the opposite. All heatmap rows are normalized
based on the total number of occurrences of the respective takeover. As
timestamps of the occurrences, we take, on the one hand, the first point at
which a takeover is happening by the finally dominating hashtag. For both
scenarios, this is when a non-initial hashtag takes over the initial
hashtag(s). On the other hand, for the failed hijacking scenario, we also keep
track of the last time a ”retake” of one of the initial hashtags occurs. At
this point, one of the initial hashtags is reclaiming its position as the most
prevalent hashtag in the conversation. As shown by Figure 9, many initial
takeovers happen early in the conversations. A slight shift towards the end of
the conversation can be observed for the failed hijacking scenario.
Nevertheless, the ”competition” for the most prevalent position is not
finished at this point. Towards the end of the conversations, most retakes
occur.

Figure 9: Initial hashtag overtakes in hijacking scenarios mainly happen at
the beginning of a conversation. Unlike in successful scenarios, in failed
ones and for retakes, the overtake tends to occur slightly later in the
conversation.

##  V Conclusion and Future Work

Existing work that focuses on analyzing conversations in online social
networks only partially considers the three dimensions of a conversation:
structure, content and dynamics. To address this shortcoming, the CODY model
is presented in this work. Based on the concept of temporal heterogeneous
information networks, it allows to model a conversation’s topological
structure, incorporate various content objects, and consider its temporal
evolution. To demonstrate the model’s capabilities, results of several
experiments are presented. These experiments cover the examination of posting
activity across the lifetime of a conversation, the study of our novel
temporal Wiener index metric, and the analysis of various ”hashtag hijacking”
scenarios. Further work will employ more advanced natural language processing
methods to improve the semantic understanding of analyzed conversations. Also,
framing the hashtag hijacking phenomenon as a prediction task might be an
interesting future extension. In this regard, one might also examine the
relationship between a conversation’s structural and content-wise evolution.

## References

  * [1] M. Brambilla, A. Javadian Sabet, K. Kharmale, and A. E. Sulistiawati, “Graph-Based Conversation Analysis in Social Media,” _Big Data and Cognitive Computing_ , vol. 6, no. 4, p. 113, 2022.
  * [2] M. Saveski, B. Roy, and D. Roy, “The structure of toxic conversations on Twitter,” in _Proc. of the Web Conference 2021_ , 2021, pp. 1086–1097.
  * [3] P. Cogan, M. Andrews, M. Bradonjic, W. S. Kennedy, A. Sala, and G. Tucci, “Reconstruction and analysis of Twitter conversation graphs,” in _Proc. of the First ACM International Workshop on Hot Topics on Interdisciplinary Social Networks Research_ , 2012, pp. 25–31.
  * [4] Y. Sun, J. Tang, J. Han, M. Gupta, and B. Zhao, “Community evolution detection in dynamic heterogeneous information networks,” in _Proc. of the Eighth Workshop on Mining and Learning with Graphs_ , 2010, pp. 137–146.
  * [5] Y. Li, Z. Lou, Y. Shi, and J. Han, “Temporal motifs in heterogeneous information networks,” in _MLG Workshop@KDD_ , 2018.
  * [6] A. Milani Fard, E. Bagheri, and K. Wang, “Relationship Prediction in Dynamic Heterogeneous Information Networks,” in _Advances in Information Retrieval_. Cham: Springer, 2019, pp. 19–34.
  * [7] N. Botzer and T. Weninger, “Entity graphs for exploring online discourse,” _Knowledge and Information Systems_ , pp. 1–19, 2023.
  * [8] B. Mathew, T. Chakraborty, N. Ganguly, and S. Datta, “Mining Twitter conversations around E-commerce promotional events,” in _Proc. of the 19th ACM Conference on Computer Supported Cooperative Work and Social Computing Companion_ , 2016, pp. 345–348.
  * [9] A. Guille, H. Hacid, C. Favre, and D. A. Zighed, “Information diffusion in online social networks: A survey,” _ACM Sigmod Record_ , vol. 42, no. 2, pp. 17–28, 2013.
  * [10] Y. Yujie, “A survey on information diffusion in online social networks,” in _Proc. of the 2020 European Symposium on Software Engineering_ , 2020, pp. 181–186.
  * [11] A. Cuzzocrea and F. Folino, “Community evolution detection in time-evolving information networks,” in _Proc. of the Joint EDBT/ICDT 2013 Workshops_ , 2013, pp. 93–96.
  * [12] S. Sajadmanesh, S. Bazargani, J. Zhang, and H. R. Rabiee, “Continuous-time relationship prediction in dynamic heterogeneous information networks,” _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , vol. 13, no. 4, pp. 1–31, 2019.
  * [13] X. Wang, Y. Lu, C. Shi, R. Wang, P. Cui, and S. Mou, “Dynamic heterogeneous information network embedding with meta-path based proximity,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 34, no. 3, pp. 1117–1132, 2020.
  * [14] R. Bian, Y. S. Koh, G. Dobbie, and A. Divoli, “Network embedding and change modeling in dynamic heterogeneous networks,” in _Proc. of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , 2019, pp. 861–864.
  * [15] T. Xie, Y. Xu, L. Chen, Y. Liu, and Z. Zheng, “Sequential recommendation on dynamic heterogeneous information network,” in _2021 IEEE 37th International Conference on Data Engineering (ICDE)_. IEEE, 2021, pp. 2105–2110.
  * [16] Y. Jia, Y. Wang, X. Jin, Z. Zhao, and X. Cheng, “Link inference in dynamic heterogeneous information network: A knapsack-based approach,” _IEEE Transactions on Computational Social Systems_ , vol. 4, no. 3, pp. 80–92, 2017.
  * [17] C. Aggarwal, Y. Xie, and P. S. Yu, “On dynamic link inference in heterogeneous networks,” in _Proc. of the 2012 SIAM International Conference on Data Mining_. SIAM, 2012, pp. 415–426.
  * [18] Y. Sun and J. Han, “Mining heterogeneous information networks: a structural analysis approach,” _ACM SIGKDD Explorations Newsletter_ , vol. 14, no. 2, pp. 20–28, 2013.
  * [19] H. Wiener, “Structural determination of paraffin boiling points,” _Journal of the American chemical society_ , vol. 69, no. 1, pp. 17–20, 1947.
  * [20] T. König, W. J. Schünemann, A. Brand, J. Freyberg, and M. Gertz, “The EPINetz Twitter Politicians Dataset 2021. A New Resource for the Study of the German Twittersphere and Its Application for the 2021 Federal Elections,” _Polit. Vierteljahresschrift_ , vol. 63, no. 3, pp. 529–547, 2022.
  * [21] A. Maireder and J. Ausserhofer, “Political discourses on Twitter: Networking topics, objects, and people,” _Twitter and society_ , vol. 89, pp. 305–318, 2014.
  * [22] C. Koylu, “Modeling and visualizing semantic and spatio-temporal evolution of topics in interpersonal communication on Twitter,” _International Journal of Geographical Information Science_ , vol. 33, no. 4, pp. 805–832, 2019.
  * [23] M. Newville, T. Stensitzki, D. B. Allen, M. Rawlik, A. Ingargiola, and A. Nelson, “LMFIT: Non-linear least-square minimization and curve-fitting for Python,” _Astrophysics Source Code Library_ , pp. ascl–1606, 2016.
  * [24] G. Csardi, T. Nepusz _et al._ , “The igraph software package for complex network research,” _InterJournal, complex systems_ , vol. 1695, no. 5, pp. 1–9, 2006.


*format=pnasformat

# A COVINDEX based on a GAM beta regression model with an application to the
COVID-19 pandemic in Italy

Luca Scrucca
Dipartimento di Economia
Università degli Studi di Perugia
Via A. Pascoli 20, 06123 Perugia, Italy
<EMAIL_ADDRESS>
https://orcid.org/0000-0003-3826-0484

###### Abstract

Detecting changes in COVID-19 disease transmission over time is a key
indicator of epidemic growth. Near real-time monitoring of the pandemic growth
is crucial for policy makers and public health officials who need to make
informed decisions about whether to enforce lockdowns or allow certain
activities. The effective reproduction number $R_{t}$ is the standard index
used in many countries for this goal. However, it is known that due to the
delays between infection and case registration, its use for decision making is
somewhat limited. In this paper a near real-time COVINDEX is proposed for
monitoring the evolution of the pandemic. The index is computed from
predictions obtained from a GAM beta regression for modelling the test
positive rate as a function of time. The proposal is illustrated using data on
COVID-19 pandemic in Italy and compared with $R_{t}$. A simple chart is also
proposed for monitoring local and national outbreaks by policy makers and
public health officials.

Abstract

Keywords pandemic surveillance; GAM beta regression; COVINDEX; public-health
decision-making

## 1\. Introduction

The World Health Organization (WHO) declared coronavirus disease (COVID-19) a
pandemic on 11 March 2020. Since then, most countries around the world have
addressed this threat by implementing various strategies to fight the
pandemic. From simple preventive measures, such as case identification and
contact tracing, quarantine and isolation, to more severe strategies based on
general lockdowns of all non-essential economical and social activities. Since
public health decision-making requires the balancing of numerous, and often
conflicting, factors, a timely and data-informed decision making process
appears crucial.

Several studies have been recently devoted to the analysis of COVID-19 data.
Referring to the Italian situation, Sebastiani et al., (2020) evaluated the
impact of government measures on the evolution of pandemic. Girardi et al.,
(2020) used robust dose-response curves to predict the contagion dynamics of
COVID-19, while Alaimo Di Loro et al., (2021) proposed an extended Generalized
Linear Model based on the Richards’ curve to model and predict incidence
indicators. A Poisson autoregressive model was discussed by Agosto et al.,
(2021) to monitor the time evolution of the COVID-19 contagion curve, while
Bartolucci and Farcomeni, (2021) introduced a spatio-temporal model based on
discrete latent variables for the analysis of weekly positive rates. Finally,
Farcomeni et al., (2021) investigated an ensemble approach for short-term
prediction of occupancy of intensive care units due to COVID-19 outbreak.

The basic reproduction number, $R_{0}$, is an indicator of the epidemic’s
virulence. It is defined as the average number of infections caused by an
infected person when the whole population is susceptible, and for SARS-CoV-2
is between 2 and 3 (Li et al.,, 2020; Hilton and Keeling,, 2020). As the
pandemic evolves, the effective reproduction number $R_{t}$ is a more useful
measure. This is the average number of infections that an infected person will
cause. An $R_{t}$ above 1.0 indicates that the outbreak is growing, and below
1.0 means that it is shrinking. As a simple understood measure, $R_{t}$ is
regularly published and discussed by the media, and it has been used in many
countries, including Italy, to decide whether to tighten or loosen control
measures. However, $R_{t}$ suffers from several drawbacks when used to monitor
the transmission of the disease over time, the main one being the delay with
which it signals the evolution of the pandemic (Gostic et al.,, 2020; Adam,,
2020). Therefore, with a delay on the estimate of $R_{t}$ between ten days to
two weeks, the use of $R_{t}$ as a near real-time decision-making tool appears
rather pointless. For further discussion on the risks caused by the misuse of
the reproduction number in the COVID-19 surveillance see Maruotti et al.,
(2021).

This paper introduces a COVID-19 index, called COVINDEX, which tries to assess
whether the epidemic is growing, shrinking, or holding steady. The proposed
index is estimated by modelling the test positive rate (TPR) with a GAM beta
regression model. TPR is an easily computed statistic, defined as the fraction
of all COVID-19 tests performed on a given day that are actually positive.
This metric can be used to understand the spread of the virus, but it also
offers a measure of how adequately a country is testing. TPR can be high if
the number of positive tests is too high, but also if the number of total
tests is too low. Most developed countries faced limited testing capacity
during the initial phase of the pandemic, which resulted in high TPR values
due to testing conducted primarily on symptomatic individuals. In the
following months the ability to administer tests using PCR (polymerase chain
reaction) or molecular swabs largely increased, leading to a situation that
allows both symptomatic and asymptomatic individuals to be tested. Although
TPR can’t be used for estimating incidence of the virus in the general
population, a fundamental epidemic parameter that would require a carefully
designed sampling plan, it can be used for monitoring the evolution of
infection and transmission in the community. Higher positive rates suggest the
need for further restrictions, such as wearing masks and physical distancing,
to slow down the spread of the disease. As a rule of thumb, World Health
Organization recommended 5% as the threshold for the percent positive rate to
declare the COVID-19 transmission under control.

The main advantages of the proposed COVINDEX is the use of data routinely
collected and its timely estimation which provides a near real-time tool to
assess the effectiveness of interventions and to inform policy. Furthermore,
since it is based on a statistical model, the associated uncertainty can be
estimated.

The paper is organized as follows. Section 2 introduces the GAM beta
regression model, its estimation and uncertainty assessment. Section 3
describes the proposed COVINDEX and its usage for monitoring the pandemic
evolution. Section 4 includes a detailed analysis of COVID-19 pandemic in
Italy, including the estimation of COVINDEX, from early March 2020 to the end
of March 2021. Section 5 contains a comparison between the proposed COVINDEX
and the effective reproduction number, showing the advantages of COVINDEX as
near real-time monitoring tool. The final section provides some concluding
remarks.

## 2\. Statistical Model for the Test Positive Rate

### 2.1. GAM beta regression

Let $y_{t}$ be the test positive rate (TPR), defined as the ratio of the
number of new positive cases $P_{t}$ to the number of tests $T_{t}$ at time
$t$, with $t$ assuming integer values between $1$ and $n$, respectively, for
the first and last day of the analyzed period. As a proportion TPR is
naturally limited in the range $[0,1]$. Several approaches and models can be
used for response variables that are expressed as proportions (Douma and
Weedon,, 2019), and perhaps the most popular statistical model is the beta
regression model (Ferrari and Cribari-Neto,, 2004; Zeileis and Cribari-Neto,,
2010).

Assume that TPR can be modelled by a beta distribution written as

$y_{t}\sim\text{Beta}(\mu_{t},\phi),$

with mean and variance of the beta distribution given, respectively, by

$\displaystyle\operatorname{\mathbb{E}}[y_{t}]$ $\displaystyle=\mu_{t},$ and
$\displaystyle\operatorname{\mathbb{V}}[y_{t}]$
$\displaystyle=\frac{\mu_{t}(1-\mu_{t})}{1+\phi}.$

Strictly speaking, the beta distribution can only model data in the open set
$(0,1)$. If extreme values 0 and 1 can actually be observed, the inflated
zero- and/or one beta distribution of Ospina and Ferrari, (2010) could be
used. Since in practice TPR rarely assumes a value of 0 and almost never 1, if
needed, the simple approach proposed by Smithson and Verkuilen, (2006) can be
adopted by applying the data transformation $(y_{t}(n-1)+0.5)/n$. The latter
is the approach followed in this paper.

The mean $\mu_{t}$ can be expressed as a function of the linear predictor
$\eta_{t}=\boldsymbol{\beta}{}^{\\!\top}\boldsymbol{x}_{t}$, where
$\boldsymbol{\beta}$ is a $(p+1)$-dimensional vector of unknown regression
coefficients (including the intercept), and $\boldsymbol{x}_{t}$ is the vector
of observed values on $p$ predictors plus a one for the intercept. Usually,
the logistic function is used in beta regression, so we can write

$\mu_{t}=\operatorname{logistic}(\eta_{t})=\frac{\exp(\eta_{t})}{1+\exp(\eta_{t})}=\frac{1}{1+\exp(-\eta_{t})}.$

The inverse of the logistic function is the logit function, the so-called
_link_ function in GLM terminology (McCullagh and Nelder,, 1989), given by:

$\operatorname{logit}(\mu_{t})=\log\left(\frac{\mu_{t}}{1-\mu_{t}}\right)=\eta_{t}.$

Generalized Additive Models (GAMs; Hastie and Tibshirani,, 1990) allows to
model the dependence of the response variable in a flexible way using smooth
functions of the predictors by defining the linear predictor as

$\eta_{t}=\beta_{0}+\sum_{j=1}^{p}f_{j}(x_{tj}),$

where $f_{j}(x_{tj})=\sum_{k=1}^{K_{j}}\beta_{jk}B_{jk}(x_{tj})$ is the
smoothing term for the $j$th predictor with $\\{B_{jk}()\\}_{k=1}^{K_{j}}$ a
set of known basis functions associated to unknown parameters $\beta_{jk}$.

Several smoothers can be defined by adopting different basis functions, such
as penalized regression splines, cubic regression splines, etc. For an
overview of the several smoothing functions available using splines bases see
Wood, (2017, Chapter 5). Among the various possibilities, thin plate
regression splines (TPRS; Wood,, 2003) represents a convenient form because
TPRS (i) do not require to specify the “knots”, (ii) use a low rank
approximation of the full basis expansion, and (iii) are isotropic smoothers,
so they are unaffected by any rotation or reflection of the covariates.

In our application the only feature included as smoothing term in the linear
predictor is time, so $x_{1}$ is an integer counting the days since the first
day of the analysis. To some extent, the coding of such feature has no
practical consequence, and other equivalent forms could have been used as
well. In addition, to account for the reduced tracing activity during weekends
(Saturday and Sunday) and holidays, a dummy variable $x_{2}$ is included
taking value 1 for data referring to weekends or holidays, and 0 otherwise.
The rationale behind the inclusion of such term is that the number of swabs
processed is noticeably limited during weekends and holidays, so a significant
increase in the test positivity rate is often observed due to the limited
testing capacity and the higher probability of testing only symptomatic cases.

Thus, in our case the linear predictor of GAM simplifies to

$\eta_{t}=\beta_{0}+\sum_{k=1}^{K}\beta_{1k}B_{1k}(x_{t1})+\beta_{2}x_{t2},$

where $\\{B_{1k}\\}_{k=1}^{K}$ represents the basis of thin plate regression
splines. Note, however, that other smooth functions would have given nearly
equivalent results.

### 2.2. Estimation

Estimation of the GAM model introduced in previous section can be pursued by
REstricted Maximum Likelihood (REML), which amounts to maximize the penalized
log-likelihood

$\ell_{P}(\boldsymbol{\beta})=\ell(\boldsymbol{\beta})-\frac{1}{2}\lambda\boldsymbol{\beta}{}^{\\!\top}\boldsymbol{S}\boldsymbol{\beta},$
(1)

where $\ell(\boldsymbol{\beta})=\sum_{t=1}^{n}\ell(y_{t}|\boldsymbol{\beta})$
is the log-likelihood for the observed values $y_{t}$ of the response
variable. The last term in the right-hand side represents the smoothing
penalty, with $\lambda$ a smoothing parameter and $\boldsymbol{S}$ a known
penalty matrix.

As reported in Wood, (2011) and Wood et al., (2016), REML is equivalent to
marginal likelihood estimation of $\boldsymbol{\beta}$ when the model contains
Gaussian random effects, and it also leads to more stable estimates of
$\lambda$ with much reduced risk of under-smoothing compared to GCV.
Furthermore, as discussed in the next section, the REML estimates of
regression coefficients have an asymptotically MAP Bayesian interpretation
that is very useful for obtaining simulated credible intervals for
predictions. For a recent review on inference and computation in GAMs see
Wood, (2020).

The selection of the smoothing parameter can be obtained, among many other
proposals, by minimizing the conditional Akaike’s information criterion (AIC).
This version of AIC for GAMs uses the log-likelihood evaluated at the
penalized MLE, and with the effective degrees of freedom computed as discussed
in Wood et al., (2016).

However, because the number of administered swabs is not constant over time,
we must take into account this fact when modelling the test positive rate.
There are several reasons for this empirical evidence. First of all, during
the weekends (particularly on Sundays) and holidays the number of swabs drops
drastically. Furthermore, during periods of strong expansion of the pandemic,
the monitoring system is unable to carry out effective surveillance and only
symptomatic patients are likely to be tested. Accounting for the different
number of swabs in the model for the positive rate can be achieved by adopting
a weighted penalized log-likelihood criterion. This amounts to replace the
log-likelihood $\ell(\boldsymbol{\beta})$ in (1) with the weighted version

$\ell_{W}(\boldsymbol{\beta})=\sum_{t=1}^{n}w_{t}\ell(y_{t}|\boldsymbol{\beta}),$

where $w_{t}$ are prior weights specifying the contribution of each data point
to the log-likelihood. In particular, indicating with $\bar{T}$ the average
number of administered swabs over the period, weights can be defined as
$w_{t}=T_{t}/\bar{T}$ so that positive rates $y_{t}$ computed from number of
swabs larger than the average have proportionally larger weights, and vice
versa for those rates based on number of swabs smaller than the average.
Furthermore, with the adopted definition for the weights the contribution of
each datum is specified without changing the overall magnitude of the log-
likelihood.

Once the model is fitted, the predicted TPR can be computed as

$\widehat{\mu}_{t}=\operatorname{logistic}\left(\widehat{\beta}_{0}+\sum_{k=1}^{K}\widehat{\beta}_{1k}B_{1k}(x_{t1})+\widehat{\beta}_{2}x_{t2}\right).$
(2)

On certain occasions, for instance when computing the COVINDEX discussed in
Section 3, we may want to compute predictions for the TPR with the
weekends/holidays effect ruled out. This is easily accomplished by setting
$x_{t2}=0$ for all $t$.

### 2.3. Uncertainty and inference

The penalized likelihood approach described above has also a Bayesian
interpretation by assuming an improper multivariate normal prior on
$\boldsymbol{\beta}$. In this case, the REML estimates of $\boldsymbol{\beta}$
coefficients are asymptotically the maximum a posteriori (MAP) of the Bayesian
posterior distribution, with the latter given by

$\boldsymbol{\beta}|(\boldsymbol{y},\lambda)\sim
N(\widehat{\boldsymbol{\beta}},(\widehat{\boldsymbol{I}}+\lambda\boldsymbol{S})^{-1}),$
(3)

where $\widehat{\boldsymbol{I}}$ is the observed information matrix (Hessian
of the negative log-likelihood) at $\widehat{\boldsymbol{\beta}}$ (Wood,,
2017, Section 6.10). This result is useful for computing approximate credible
intervals for any function of $\boldsymbol{\beta}$ by simulating from the
posterior (Gelman and Hill,, 2006, Section 7.2). Wood, (2017, p. 294) reported
good frequentist coverage properties for such Bayesian credible intervals,
with empirical coverage close to the nominal level when averaged across the
domain of the function.

In practice, coefficients $\boldsymbol{\beta}^{*}$ are simulated from (3), and
then plugged in equation (2) to get the simulated means $\mu_{t}^{*}$. The
process is replicated a large number of times, say $10\,000$ or more, and the
percentiles of the simulated distributions at different values of $x_{t}$ can
be used to compute the limits of approximate credible intervals for the mean.
To compute approximate credible intervals for the single prediction we
simulate response values as
$y^{*}_{t}\sim\text{Beta}(\mu^{*}_{t},\widehat{\phi})$, where $\mu_{t}^{*}$ is
the simulated mean as described above, $\widehat{\phi}$ is the model estimate
of the precision parameter, and then we compute the percentiles of the
simulated distribution of predicted values for the response. The empirical
coverage of the prediction intervals will be assessed in Section 4.2.

## 3\. COVINDEX as a Monitoring and Decision-Making Tool

The COVINDEX proposed in this paper is an attempt to compute a synthetic index
summarizing the evolution of the COVID-19 pandemic, which can be useful to
policy makers and public health officials for monitoring local and national
outbreaks. In our proposal this is simply computed as

$\text{COVINDEX}_{t}=\frac{\widehat{\mu}_{t}}{\widehat{\mu}_{t-7}},$ (4)

the ratio of the predicted positive rate at time $t$ to the prediction 7 days
earlier. The value of 7 is chosen because it is approximately the expected
incubation time for COVID-19 (Nazar and Elfadil,, 2021), and because it
corresponds to the observed weekly fluctuation in testing. A COVINDEX value
larger than 1.0 means that the pandemic is growing, while a value smaller than
1.0 indicates that new infections are slowing down.

The COVINDEX estimate is clearly affected by uncertainty and to account for it
the approach outlined in Section 2.3 for TPR can be used here as well. In
particular, for each simulated series of values $\mu_{t}^{*}$, simulated
COVINDEX series can be obtained by applying equation (4) to get the simulated
values $\text{COVINDEX}^{*}_{t}=\mu^{*}_{t}/\mu^{*}_{t-7}$. Approximate
credible intervals can then be computed from the percentiles of the simulated
distribution.

We argue that decisions made by policy makers should be based both on the
COVINDEX, which provides an outlook on the likely behaviour of the pandemic in
the near future, and on the level of the estimated TPR, which represents its
current status. Following this idea, a TPR-COVINDEX risk quadrant chart can be
drawn (see Figure 1). This chart illustrates four potential scenarios which
represent a useful tool for a decision maker. The quadrants are defined by the
dashed lines drawn at selected threshold values. For COVINDEX the natural
reference value is 1.0, with values below it indicating a shrinking outbreak,
and values higher than 1.0 indicating epidemic situations that are
increasingly worrying and out of control. Note that, since the index is a
ratio, the $y$-axis is expressed in logarithmic scale. For the positive rate,
the threshold value can be set according to the World Health Organization,
which published a set of criteria to inform whether the epidemic is under
control. In particular, one criterion states that “[…] less than 5% of samples
positive for COVID-19, at least for the last 2 weeks, assuming that
surveillance for suspected cases is comprehensive” (World Health
Organization,, 2019).

According to the above mentioned threshold values, the upper-right quadrant
represents the worst-case scenario, with high values of both TPR and COVINDEX.
On the contrary, the best-case scenario is the lower-left quadrant which has
both low TPR and COVINDEX less than 1.0 indicating a decreasing circulation of
the virus. The remaining quadrants are intermediate cases. Typical situations
will move in a clockwise direction, moving from the worst-case, represented by
the red quadrant on top-right, to the orange quadrant at bottom-right, and
eventually reaching the yellow quadrant indicating an outbreak under control.
However, in some cases the pandemic could regain strength by getting COVINDEX
values greater than 1.0, thus moving towards the top-left orange quadrant or
directly towards the worst-case situation described by the red quadrant. A
description of the Italian situation since March 2020 is discussed in Section
4.

Figure 1: TPR-COVINDEX risk quadrant chart.

## 4\. Application to Italian COVID-19 Pandemic

### 4.1. Data

The Italian Department of Protezione Civile provides daily information on the
COVID-19 pandemic, both at the national and the regional level, in a public
GitHub repository (Presidenza del Consiglio dei Ministri – Dipartimento della
Protezione Civile,, 2020). Among the data contained in this repository, the
cumulative number of naso-pharyngeal or molecular swabs and the corresponding
positive tests are provided. Starting with January 15th, 2021, antigen tests
are also officially recorded, while previously only some regions included them
in the recorded statistics since autumn 2020. The reliability of such
information is at best questionable and not available uniformly for the year
2020. For these reasons, in our analyses we considered the information from
daily molecular swabs (not persons tested) to compute the test positive rate
(TPR), a commonly used screening and diagnostic tool for COVID-19 (World
Health Organization,, 2020).

The plot on Figure 2 shows the observed TPR over time with points proportional
to the administered swabs.

Figure 2: Plot of test positive rate from beginning of COVID-19 pandemic in
Italy to the end of observational period with size of points proportional to
the number of molecular swabs administered.

### 4.2. GAM beta regression model estimate

Table 1 reports the summary output of the estimated GAM beta regression model
for the test positive rate in Italy from March 1st 2020 to June 30th 2021. The
parametric terms include the intercept and a dummy variable for the days
following the weekends (Saturday and Sunday) and holidays. The smooth term
captures the evolution of underlying trend in the observed test positive rate.
The amount of smoothing applied to the time predictor is selected by
minimizing the AIC, as shown in Figure 3. The graphs of the autocorrelation
and partial autocorrelation functions for the deviance residuals in Figure 4
show no significant remaining correlation at different lags.

Table 1: GAM beta regression model summary. Num. of obs. = 487 | Dispersion par. = 1467.3
---|---
Log-likelihood = 1828.7 | Deviance expl. = 0.9891
REML = 1726.3 | AIC = -3583.1
Parametric coefficients:
| Estimate | Std. error | z-value |  $p$-value
(Intercept) | -3.1904 | 0.01177 | -271.07 | < 0.001
Weekend | 0.1681 | 0.01023 | 16.44 | < 0.001
Smooth terms:
| edf | Ref. df | ChiSq-value |  $p$-value
$s(t)$ | 35.11 | 39.67 | 16820 | < 0.001
Figure 3: Trace plot of AIC as a function of $k$, the number of basis
functions used by the thin plate regression smoother, with the corresponding
effective degrees of freedom (EDF). The EDF expresses the complexity of the
smoother, with larger values indicating more wiggly smoothers, and it cannot
be larger than $k$. The vertical dashed line is drawn at the minimum AIC used
for the selection of the final model. Figure 4: Autocorrelation and partial
autocorrelation functions for the deviance residuals of the estimated GAM beta
regression in Table 1

Figure 5 reports the estimated curve for the test positive rate with 95%
credible intervals for the mean and the single value obtained by simulating
from the posterior distribution as described in Section 2.3. The highest
positive rates are achieved in March 2020 during the first wave of pandemic,
and on November 2020, corresponding to the second wave. A resurgence of spread
during the end of 2020 is followed by a quick decrease in earlier 2021.
Subsequently, the situation remained stable for about a month, but in the
second part of February another sharpe increase occurred due to the appearance
of COVID-19 variants in the Italian territory (in particular the Alpha or
english variant). Starting from the beginning of April, a marked decline in
the TPR can be observed, likely favored by the increased full vaccination
coverage of the Italian population.

Figure 5: GAM estimate of test positive rate (blue line) in Italy during 2020
and first half of 2021, with 95% simulated credible intervals for the mean
(dark grey area) and for the single value (light grey area).

### 4.3. Empirical coverage of credible intervals for predictions

To investigate the accuracy of the simulated prediction intervals, we
considered all the dates from February 1st 2021 to May 31st 2021 and the time
horizon for the predictions from 1 day to 14 days. For each date we fitted the
GAM beta regression model using all the data available up to that day, and
then we used the estimated model to compute the simulated credible intervals
for the following 14 days. This process was replicated for the all the days in
the specified range and the empirical coverage calculated. Figure 6 reports on
the left panel a graph showing the inclusion or exclusion of the observed
values of TPR in the simulated prediction intervals, while on the right a plot
of the empirical coverage for the time horizon from 1 up to 14 days. Overall
the coverage is close to the nominal level, with all the values above 90% for
the forecasts of the first week, and between 85% and 90% for the second week.
It is interest to note that most of the coverage errors occur in periods of
abrupt changes, for instance at the sharp rise of TPR in the last week of
February or at the beginning of TPR decline in mid-March. As expected, the
empirical coverage decreases as the prediction horizon increases.

Figure 6: Prediction intervals coverage for the GAM beta regression model of
TPR at 95% nominal level. The graph on the left shows the inclusion (blank
cells) or exclusion (red cells) of the observed value of TPR for the
prediction intervals in subsequent dates corresponding to the time horizon
from 1 up to 14 days. The empirical coverage percentages are shown on the
right graph, with the vertical dashed line representing the nominal level.

### 4.4. COVINDEX estimate

Based on the estimated model and uncertainty for the test positive rate, the
COVINDEX is computed following equation (4). Figure 7 shows the estimated
COVINDEX with 95% credible intervals. Notice that the $y$-axis is expressed on
logarithmic scale, the natural scale to visualize ratios (Wilke,, 2019, , Sec.
3.2). The index fluctuates widely throughout the year 2020, following the
periods of expansion and contraction of the spread of the pandemic. After the
first wave in spring 2020 we observe a quick decreasing trend, followed by a
slowly increase during the summer, corresponding to a relaxation of the
containment measures, with values significantly larger than 1.0 during August.
This represents the first signal of a resurgence of the pandemic. Sharpe and
large increases are also observed during October in conjunction with the
second wave that strongly affected Italy. Two additional peaks are detected at
the end of 2020 and on February 2021, corresponding to gradual relaxation of
containment measures before Christmas and mid January, with the latter that
occurred during the period of political instability associated with the change
of government.

Figure 7: Evolution of COVINDEX (on logarithmic scale) for Italy during 2020
and first half of 2021, with approximate 95% simulated credible intervals.

From the adopted definition of equation (4), COVINDEX is computed by taking
the ratio of the estimated TPR with respect to a 7-days-before estimate. In
Section 3 we provide the rationale for this choice. However, it may be
interesting to investigate how the index changes assuming different lags.
Figure 8 shows the COVINDEX estimates obtained when different lag values are
used. The general behaviour of the curves is similar across different lags,
but the amplitude of the oscillation increases as the lag increases. This
appears reasonable since, essentially, COVINDEX compares the estimated TPR at
time $t$ with the value at time $t-\text{lag}$. Thus, for smaller values of
the lag the index fluctuates less and is more stable than at higher lag
values. However, lag values that are too small cannot highlight the dynamics
of TPR because too close values tend to be quite similar. In this sense, the
selected 7-day lag appears to be a sensible choice.

Figure 8: A comparison of COVINDEX (on logarithmic scale) computed at
different lags for Italy on 2020 and first half of 2021.

### 4.5. TPR-COVINDEX risk quadrant chart analysis

Figure 9 shows the TPR-COVINDEX risk quadrant chart for Italy, with points
connected following the temporal path, a graph also known as connected
scatterplot (Haroz et al.,, 2015). The curve concisely represents the
evolution of both indices during the pandemic. Starting with the critical
situation in March 2020, the situation improved in the following months,
moving from the red quadrant to the orange bottom-right quadrant and then the
yellow quadrant during summer 2020. By the end of summer 2020 we observe a
worsening of the situation that lead to the red quadrant in November. In the
following months there has been a constant oscillation between the red and
right-orange quadrants, indicating a serious pandemic situation.

Figure 9: TPR-COVINDEX risk chart as connected scatterplot for Italy. The
first day of each month is highlighted to provide a temporal reference.

A scatterplot of TPR vs COVINDEX is also useful for surveillance of the
pandemic in different Italian regions. Figure 10 summarizes the status of the
pandemic for the Italian regions at selected time points. A high-risk
situation is observed at the beginning of November 2020, where all regions
belong to the red quadrant. The following month saw an improvement with most
regions moving towards the bottom-right orange quadrant. A more complex and
varied situation is observed between February and March 2021, with some
regions moving from the red to the orange quadrant, and vice versa for other
regions.

Figure 10: TPR-COVINDEX plot for Italian regions at different time points.

## 5\. A comparison of COVINDEX with the effective reproduction number

The main index used in Italy for pandemic surveillance is $R_{t}$, the
effective reproduction number. The procedure employed for estimating $R_{t}$
is described by Guzzetta and Merler, (2020) and it is based on the Bayesian
methodology of Cori et al., (2013). Details can be found at
https://www.epicentro.iss.it/coronavirus/sars-cov-2-sorveglianza-dati. An
archive containing both the data and the R script used by the Italian National
Institute of Health (ISS) for computing $R_{t}$ is available at
https://www.epicentro.iss.it/coronavirus/open-data/calcolo_rt_italia.zip.

In this Section we provide a comparison of the proposed COVINDEX with the
values of $R_{t}$ estimated following the procedure outlined above for Italy
from March 2020 to June 2021. Furthermore, since the effective reproduction
number does not provide a timely snapshot of the evolution of the pandemic, we
also provide two examples showing the failure of $R_{t}$ to highlight the
likely evolution of the pandemic and we compare its behaviour with the
proposed COVINDEX.

The top graph reported in Figure 11 shows the estimated curves for COVINDEX
and $R_{t}$. Overall, a similar trend can be observed for the two curves,
particularly since October 2020. $R_{t}$ appears to be more wiggly than
COVINDEX, especially during the summer 2020. Likely, this is related to the
large uncertainty in that period due to the relative small number of positive
cases (around few hundreds) observed in that period. One of main drawbacks of
using $R_{t}$ for real-time monitoring is shown in the final part of the
graph. In fact, if at the end of the June the COVINDEX curve seems to suggest
a resumption of the pandemic, the $R_{t}$ index continues to show a decreasing
trend This behaviour can also be seen in other time periods, as discussed
below.

Figure 11: Comparison of COVINDEX and $R_{t}$ for Italy. Bottom panels show
the comparison at early December 2020 (left) and at the end of February 2021
(right).

As mentioned in Section 4.4, the second wave of COVID-19 epidemic hit Italy
between the second half of October and the beginning of November 2020,
followed by a rapid decrease during the remainder of the month. However, from
the beginning of December 2020 it was evident that this decline had stopped
and that the situation was starting to get worse. This is clearly indicated by
the upward slope of the COVINDEX computed on December 5th and shown in the
bottom-left graph of Figure 11. On the contrary, the $R_{t}$ index calculated
on the same day, with estimates considered valid up to 14 days before,
produces a curve which erroneously suggests a decline in the spread of the
pandemic. However, if the $R_{t}$ curve is estimated a week later, we begin to
see an increase in the spread of the pandemic (see the dotted red curve in
bottom-left graph of Figure 11). The main problem is that such alert is
reported too late.

A similar situation is also faced at the beginning of March 2021. After a
period of almost constant positive rate during February 2021, with both
COVINDEX and $R_{t}$ oscillating around 1.0, by the end of the month there was
a clear increase of the test positive rate. This was immediately signalled by
COVINDEX computed on February 28th 2021 (see bottom-right graph in Figure 11),
but $R_{t}$ computed on the same day was still signalling a steady state and
only after a week an increasing value of $R_{t}$ would have signalled the
resurgence of the pandemic.

The comparison between COVINDEX and $R_{t}$ can also be conducted at the
regional level. Here we present a comparison for two Italian regions,
Lombardia and Umbria. These are two very different regions, both in size and
geographical position, but also in terms of pandemic history. If Lombardia was
the most affected region of Italy during the 1st wave of the COVID-19
pandemic, Umbria suffered only marginal effects in this phase. On the
contrary, the so-called 3rd wave that occurred in winter/spring 2021 hit
Umbria earlier than in the rest of Italian regions, including Lombardia.

Likewise the national level, there is a substantial similarity between the
trend of COVINDEX and $R_{t}$ for the two regions, with the former which
appears to have a smoother behaviour (see Figure 12). Both indices correctly
identified the peak of the pandemic in October 2020 and at the end of 2020.
But if for Umbria the beginning of 2021 is marked by the arrival of the third
wave caused by the circulation of SARS-CoV-2 variants, namely Alpha (or
English) and Gamma (or Brazilian), in Lombardia the presence of these variants
only occurred from mid-February. Subsequently, starting from spring 2021, a
decline in the epidemic can be observed in both regions.

However, there are also differences that are worth pointing out. For Lombardia
there are two values of $R_{t}$, in mid-February and early June, which appear
suspicious as they are placed outside the apparent trend, underestimating in
the first case and overestimating the trend in the second case. For Umbria,
the $R_{t}$ seems to increase starting from the last week of May, while the
COVINDEX still suggests a decreasing trend. This behaviour of $R_{t}$ is also
suspect as the TPR of the region remains substantially stable or slightly
decreasing in this period, with all the TPR values less than 1% in the last 10
days of June.

Figure 12: Comparison of COVINDEX and $R_{t}$ for the Italian regions
Lombardia (top) and Umbria (bottom), from September 2020 to June 2021.

## 6\. Final comments

In this paper we have proposed an index, named COVINDEX, that can be used for
near real-time monitoring of COVID-19 pandemic. The index is computed as the
ratio of the estimated test positive rate on a given day with respect to the
value estimated for a week before. Estimation of test positive rates is
obtained by statistical modelling the daily empirical positive rates
calculated from the observed data. To this end, a GAM beta regression model
with weights proportional to the administered tests is fitted. By exploiting
the relationship of penalized likelihood for GAMs with MAP Bayesian
estimation, credible intervals for COVINDEX can be obtained via simulation to
express the associated uncertainty.

We applied the proposed methodology to the Italian COVID-19 outbreak and we
compared the trend of COVINDEX to the effective reproduction number $R_{t}$.
The analyses carried out confirm that $R_{t}$ is a delayed index of epidemic
trend, and for this reason may provide a biased picture of the current
pandemic status. On the contrary, COVINDEX seems to provide a more up-to-date
information which can be used as a decision-making tool. This aspect is of
crucial importance for all policy makers and public health officials. We defer
to future research the evaluation of the implications deriving from the
adoption of the proposed index.

Although the main focus of the analysis in this paper was the national level,
similar considerations can be made for territorial administrative divisions,
such as regions and provinces. In these cases, however, it should be noted
that further assumptions are necessary, in particular the independence of the
epidemic trend between neighbouring territories. However, an improved approach
should account for the spatio-temporal dependency structure. For instance,
Mingione et al., (2021) fitted a spatio-temporal CAR model with spatial
dependence expressed by specifyicing an adjacency matrix derived from a
network model of links and transport exchanges among Italian regions (Della
Rossa et al.,, 2020). The study of these aspects is deferred to future
research.

All the analyses have been performed in R version 4.1.0 (R Core Team,, 2021),
using the package mgcv (Wood,, 2021) and functions written by the author. Code
to reproduce the analyses is available in a GitHub repository at
https://github.com/luca-scr/COVINDEX.

## References

  * Adam, (2020) Adam, D. (2020).  A guide to R – the pandemic’s misunderstood metric.  Nature, 583(7816):346–348.
  * Agosto et al., (2021) Agosto, A., Campmas, A., Giudici, P., and Renda, A. (2021).  Monitoring COVID-19 contagion growth.  Statistics in Medicine.
  * Alaimo Di Loro et al., (2021) Alaimo Di Loro, P., Divino, F., Farcomeni, A., Jona Lasinio, G., Lovison, G., Maruotti, A., and Mingione, M. (2021).  Nowcasting COVID-19 incidence indicators during the Italian first outbreak.  Statistics in Medicine, 40(16):3843–3864.
  * Bartolucci and Farcomeni, (2021) Bartolucci, F. and Farcomeni, A. (2021).  A spatio-temporal model based on discrete latent variables for the analysis of COVID-19 incidence.  Spatial Statistics.
  * Cori et al., (2013) Cori, A., Ferguson, N. M., Fraser, C., and Cauchemez, S. (2013).  A new framework and software to estimate time-varying reproduction numbers during epidemics.  American Journal of Epidemiology, 178(9):1505–1512.
  * Della Rossa et al., (2020) Della Rossa, F., Salzano, D., Di Meglio, A., De Lellis, F., Coraggio, M., Calabrese, C., Guarino, A., Cardona-Rivera, R., De Lellis, P., Liuzza, D., et al. (2020).  A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic.  Nature Communications, 11(1):1–9.
  * Douma and Weedon, (2019) Douma, J. C. and Weedon, J. T. (2019).  Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression.  Methods in Ecology and Evolution, 10(9):1412–1430.
  * Farcomeni et al., (2021) Farcomeni, A., Maruotti, A., Divino, F., Jona-Lasinio, G., and Lovison, G. (2021).  An ensemble approach to short-term forecast of COVID-19 intensive care occupancy in Italian regions.  Biometrical Journal, 63(3):503–513.
  * Ferrari and Cribari-Neto, (2004) Ferrari, S. and Cribari-Neto, F. (2004).  Beta regression for modelling rates and proportions.  Journal of Applied Statistics, 31(7):799–815.
  * Gelman and Hill, (2006) Gelman, A. and Hill, J. (2006).  Data analysis using regression and multilevel/hierarchical models.  Cambridge University Press.
  * Girardi et al., (2020) Girardi, P., Greco, L., Mameli, V., Musio, M., Racugno, W., Ruli, E., and Ventura, L. (2020).  Robust inference for non-linear regression models from the Tsallis score: Application to coronavirus disease 2019 contagion in Italy.  Stat, 9(1):e309.
  * Gostic et al., (2020) Gostic, K. M., McGough, L., Baskerville, E. B., Abbott, S., Joshi, K., Tedijanto, C., Kahn, R., Niehus, R., Hay, J. A., De Salazar, P. M., et al. (2020).  Practical considerations for measuring the effective reproductive number, $R_{t}$.  PLoS Computational Biology, 16(12):e1008409.
  * Guzzetta and Merler, (2020) Guzzetta, G. and Merler, S. (2020).  Stime della trasmissibilità di SARS-CoV-2 in Italia.  EpiCentro - Istituto Superiore di Sanità: https://www.epicentro.iss.it/coronavirus/open-data/rt.pdf.
  * Haroz et al., (2015) Haroz, S., Kosara, R., and Franconeri, S. L. (2015).  The connected scatterplot for presenting paired time series.  IEEE Transactions on Visualization and Computer Graphics, 22(9):2174–2186.
  * Hastie and Tibshirani, (1990) Hastie, T. J. and Tibshirani, R. J. (1990).  Generalized Additive Models, volume 43.  Chapman & Hall/CRC.
  * Hilton and Keeling, (2020) Hilton, J. and Keeling, M. J. (2020).  Estimation of country-level basic reproductive ratios for novel Coronavirus (SARS-CoV-2/COVID-19) using synthetic contact matrices.  PLoS Computational Biology, 16(7):e1008031.
  * Li et al., (2020) Li, Q., Guan, X., Wu, P., Wang, X., Zhou, L., Tong, Y., Ren, R., Leung, K. S., Lau, E. H., Wong, J. Y., et al. (2020).  Early transmission dynamics in Wuhan, China, of novel coronavirus–infected pneumonia.  New England Journal of Medicine, 382(13):1199–1207.
  * Maruotti et al., (2021) Maruotti, A., Ciccozzi, M., and Divino, F. (2021).  On the misuse of the reproduction number in the COVID-19 surveillance system in Italy.  Journal of Medical Virology, 93(5):2569–2570.
  * McCullagh and Nelder, (1989) McCullagh, P. and Nelder, J. A. (1989).  Generalized Linear Models.  Chapman and Hall, CRC, London, 2nd edition.
  * Mingione et al., (2021) Mingione, M., Di Loro, P. A., Farcomeni, A., Divino, F., Lovison, G., Maruotti, A., and Lasinio, G. J. (2021).  Spatio-temporal modelling of COVID-19 incident cases using Richards’ curve: An application to the Italian regions.  Spatial Statistics, page 100544.
  * Nazar and Elfadil, (2021) Nazar, Z. and Elfadil, A. (2021).  The estimations of the COVID-19 incubation period: A scoping reviews of the literature.  Journal of Infection and Public Health, 14(5):638–646.
  * Ospina and Ferrari, (2010) Ospina, R. and Ferrari, S. L. (2010).  Inflated beta distributions.  Statistical Papers, 51(1):111–126.
  * Presidenza del Consiglio dei Ministri – Dipartimento della Protezione Civile, (2020) Presidenza del Consiglio dei Ministri – Dipartimento della Protezione Civile (2020).  Dati COVID-19 Italia.  GitHub repository: https://github.com/pcm-dpc/COVID-19.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Sebastiani et al., (2020) Sebastiani, G., Massa, M., and Riboli, E. (2020).  COVID-19 epidemic in Italy: evolution, projections and impact of government measures.  European Journal of Epidemiology, 35(4):341–345.
  * Smithson and Verkuilen, (2006) Smithson, M. and Verkuilen, J. (2006).  A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables.  Psychological Methods, 11(1):54–71.
  * Wilke, (2019) Wilke, C. O. (2019).  Fundamentals of data visualization: a primer on making informative and compelling figures.  O’Reilly Media.
  * Wood, (2021) Wood, S. (2021).  mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation.  R package version 1.8-36.
  * Wood, (2003) Wood, S. N. (2003).  Thin plate regression splines.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(1):95–114.
  * Wood, (2011) Wood, S. N. (2011).  Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(1):3–36.
  * Wood, (2017) Wood, S. N. (2017).  Generalized Additive Models: an introduction with R.  Chapman & Hall/CRC, 2nd edition.
  * Wood, (2020) Wood, S. N. (2020).  Inference and computation with generalized additive models and their extensions.  Test, 29(2):307–339.
  * Wood et al., (2016) Wood, S. N., Pya, N., and Säfken, B. (2016).  Smoothing parameter and model selection for general smooth models.  Journal of the American Statistical Association, 111(516):1548–1563.
  * World Health Organization, (2019) World Health Organization (2019).  Public health criteria to adjust public health and social measures in the context of COVID-19.  Annex to Considerations in adjusting public health and social measures in the context of COVID-19, 12 May 2020: https://apps.who.int/iris/bitstream/handle/10665/332073/WHO-2019-nCoV-Adjusting_PH_measures-Criteria-2020.1-eng.pdf.
  * World Health Organization, (2020) World Health Organization (2020).  Considerations for implementing and adjusting public health and social measures in the context of COVID-19.  Iterim guidance, 4 November 2020: https://www.who.int/publications/i/item/considerations-in-adjusting-public-health-and-social-measures-in-the-context-of-covid-19-interim-guidance.
  * Zeileis and Cribari-Neto, (2010) Zeileis, A. and Cribari-Neto, F. (2010).  Beta regression in R.  Journal of Statistical Software, 34(2):1–24.