# Loglinear modeling with mixed numerical and categorical predictor variables
through an Extended Stereotype Model

Mark de Rooij
Methodology and Statistics department
Leiden University
The Netherlands
<EMAIL_ADDRESS>

###### Abstract

Loglinear analysis is most useful when we have two or more categorical
response variables. Loglinear analysis, however, requires categorical
predictor variables, such that the data can be represented in a contingency
table. Researchers often have a mix of categorical and numerical predictors.
We present a new statistical methodology for the analysis of multiple
categorical response variables with a mix of numeric and categorical predictor
variables. Therefore, the stereotype model, a reduced rank regression model
for multinomial outcome variables, is extended with a design matrix for the
profile scores and one for the dependencies among the responses. An MM
algorithm is presented for estimation of the model parameters. Three examples
are presented. The first shows that our method is equivalent to loglinear
analysis when we only have categorical variables. With the second example, we
show the differences between marginal logit models and our extended stereotype
model, which is a conditional model. The third example is more extensive, and
shows how to analyze a data set, how to select a model, and how to interpret
the final model.

_Keywords_ Multinomial $\cdot$ Reduced Rank Model $\cdot$ MM algorithm $\cdot$
Generalized SVD

## Declarations

Declarations of interest: none

##  1 Introduction

Categorical data are collected in many scientific disciplines and there is a
long history of statistical approaches to analyze such data (see Agresti,,
2013, chapter 17 for a historical tour). Loglinear models (Bishop et al.,,
1975; Wickens,, 1989; Agresti,, 2013) have been developed for the analysis of
categorical data which are represented as contingency tables, that are, multi-
way cross classifications of the variables with the number of occurrences of
any combination in the cells. Loglinear models are a special case of
generalized linear models (McCullagh and Nelder,, 1989) where the count
variable acts a the response and the variables as predictors in a Poisson
regression model.

When one of the variables is considered a response and the others predictors,
a loglinear analysis corresponds to a logit linear regression model. This
correspondence is often used to interpret the parameters of the loglinear
model where the focus is on the predictor-response relationship (see Agresti,,
2013, Section 9.5 for details). Logit linear regression is a special case of
logistic regression, where the latter allows for numerical predictor
variables. With numerical variables, the data can not be represented in a
contingency table anymore.

As noted by Agresti (2002, page 332) , loglinear models are most natural when
two or more variables are considered response variables. Having a single
response, it is more natural to consider logit linear or logistic regression
models directly. When some of the predictor variables are numerical, however,
loglinear models cannot be applied.

In the present paper, we generalize the loglinear model to the case where we
have numeric predictor variables or a mix of numeric and categorical predictor
variables. We focus on binary response variables, although the framework we
present can be extended to response variables with more than two categories.
Our endeavour starts with a multinomial logistic model for the response
profile, that is, any combination of the categories of the response variables.
The matrix with regression weights of such a multinomial model can be
constrained by using a design matrix for the profiles. Furthermore, a rank
constraint can be imposed on this matrix to decrease the number of parameters,
similar to reduced rank regression (Izenman,, 1975; Anderson,, 1984; Yee and
Hastie,, 2003; Reinsel et al.,, 2022). Our modelling framework has standard
logistic regression and loglinear models as special cases.

There have been other attempts to multivariate classification using mixtures
of dichotomous and continuous predictors. The location model (Olkin and Tate,,
1961; Afifi and Elashoff,, 1969; Krzanowski,, 1975, 1976, 1993) is such an
approach. In the location model, the joint distribution of the variables is
decomposed into a marginal distribution of the dichotomous variables, and a
conditional (multivariate) normal distribution of the continuous variables.
Krzanowski, (1976) shows in detail how to use the location model for
classification. The location model is a _generative_ model whereas the model
we propose in this paper is a _discriminative_ one (cf. Ng and Jordan,, 2001).
Each of these two paradigms has its advantages and disadvantages (Efron,,
1975; Rubinstein et al.,, 1997), depending on the circumstances. One advantage
of our model is that it needs no assumption about the distribution of the
continuous variables, whereas the location model assumes these variables to
follow a multivariate normal distribution.

After laying out our notation in the next section, we introduce the _extended
stereotype model_ and its properties in Section 3. Section 4, describes an
algorithm for maximum likelihood estimation of the model parameters and delves
into model selection. In Section 5, the extended stereotype model is applied
to three empirical data sets. The first analysis focuses on the relationship
between our extended stereotype model and loglinear analysis. The second
analysis contrast our approach with marginal models. The third example
discusses model selection and interpretation for a data set having five
response variables and eight predictors. We conclude this paper, in Section 6,
with a brief discussion.

##  2 Notation

For ease of reference, we will use the following notation:

  * •

Observations or participants are denoted by $i=1,\ldots,N$.

  * •

Predictor variables are indicated by $p=1,\ldots,P$.

  * •

Response variables are indicated by $r=1,\ldots,R$.

  * •

$\bm{X}$: Predictor matrix of size $N\times P$, with rows $\bm{x}_{i}$,
columns $\bm{x}_{p}$, and elements $x_{ip}$..

  * •

$\bm{Y}$: Response variable matrix of size $N\times R$, each
$y_{ir}=\\{0,1\\}$.

  * •

$\mathcal{G}$ response profile variable such that
$\mathcal{G}\in\\{1,\ldots,K\\}$, with $K=2^{R}$.

  * •

$\bm{G}$: response profile indicator matrix of size $N\times K$, with
$\max(K)=2^{R}$.

  * •

$\bm{Z}$: Design matrix of size $K\times Q$ for the profile scores. $Q$ is
dependent on the research questions. $Q\leq K-1$.

  * •

$\bm{W}$: Design matrix of size $K\times T$ for the residual dependencies. $T$
is dependent on the research questions. $T\leq K-1$.

  * •

$S^{*}$: maximum rank or dimensionality, $S^{*}=min(P,K-1)$.

  * •

$S$: chosen rank or dimensionality.

  * •

$\bm{B}_{x}$: with elements $b^{(x)}_{ps}$ a $P\times S$ matrix.

  * •

$\bm{B}_{z}$: with elements $b^{(z)}_{qs}$ a $Q\times S$ matrix.

  * •

$\bm{b}_{w}$: with elements $b^{(w)}_{t}$ a $T$ vector.

  * •

$\bm{A}$ matrix with implied regression coefficients, that is,
$\bm{A}=\bm{B}_{x}\bm{B}^{\prime}_{z}$.

  * •

$\bm{U}$: Object score matrix, equal to $\bm{U}=\bm{X}\bm{B}_{x}$.

  * •

$\bm{V}$: Response profile score matrix, equal to $\bm{V}=\bm{Z}\bm{B}_{z}$.

  * •

$\bm{m}$: Offset parameters $m_{k}$ for the $K$ profiles,
$\bm{m}=\bm{W}\bm{b}_{w}$.

##  3 The Extended Stereotype Model

###  3.1 Model

Let us first define the categorical variable $\mathcal{G}$ that summarizes the
response profile of a participant. With three dichotomous response variables,
possible categories of $\mathcal{G}$ are 000 and 101. The maximum number of
different profiles, or categories, is $2^{R}$. Sometimes, not all possible
profiles are observed and the number of categories of the variable
$\mathcal{G}$ is substantially smaller. Let us denote the number of observed
profiles by $K$. This multinomial response variable will be used in a
regression model of the form

$\pi_{k}(\bm{x}_{i})=P(\mathcal{G}=k|\bm{x}_{i})=\frac{\exp(\theta_{ik})}{\sum_{k^{\prime}=1}^{K}\exp(\theta_{ik^{\prime}})},\
\mathrm{for}\ k=1,\ldots,K.$

In the standard multinomial regression model $\theta_{ik}$ is parameterized as
$\theta_{ik}=m_{k}+\bm{x}^{\prime}_{i}\bm{a}_{k}$, where $\bm{a}_{k}$ is a
vector with regression weights for response profile $k$. Usually, one of the
elements $m_{k}$ and corresponding vector $\bm{a}_{k}$ are fixed to zero for
identification. The vector $\bm{m}$ with elements $m_{k}$ models the
(residual) dependencies among the $R$ response variables. The vectors
$\bm{a}_{k}$ can be collected in the matrix $P\times K$ matrix $\bm{A}$.

In our case, $K$ is an exponential function of the number of response
variables and grows quickly, that is, with 8 response variables $K$ is already
256 and with 11 response variables $K$ equals 2048. When $K$ is large, many
parameters need to be estimated for this model (i.e., $(P+1)\times(K-1)$),
that is, there are $K-1$ vectors $\bm{a}_{k}$ of length $P$.

The stereotype model (Anderson,, 1984), that is a reduced rank regression
model for a multinomial response variables, parameterized the $\theta_{ik}$ as

$\theta_{ik}=m_{k}+\bm{x}^{\prime}_{i}\bm{B}_{x}\bm{v}_{k},$

with $\bm{v}_{k}$ profile scores or loadings and $\bm{B}_{x}$ a matrix with
weights that find an optimal subspace of the predictors $\bm{X}$ that predict
the responses. Sometimes, the $\bm{u}_{i}=\bm{B}^{\prime}_{x}\bm{x}_{i}$ are
called object scores which are linear combinations of the predictor variables.
Let $S$ be the length of the vectors $\bm{u}_{i}$ and $\bm{v}_{k}$. When $S$
equals $S^{*}=\min(P,K-1)$ the model is equivalent to the standard multinomial
regression model. When $S<S^{*}$ the number of parameters is reduced. When the
$\bm{v}_{k}$ are concatenated into the matrix $\bm{V}$, we can write
$\bm{A}=\bm{B}_{x}\bm{V}^{\prime}$, which shows the reduced rank property.

For multivariate binary data, each of the $K$ profiles is characterised by the
underlying response variables. Let us define for response profile $k$ the
design variable $z_{kr}=\frac{1}{2}$ if profile $k$ has $Y_{r}=1$ and
$z_{kr}=-\frac{1}{2}$ if profile $k$ has $Y_{r}=0$. With these design
variables, we are able to impose a structure on the profile scores
($\bm{v}_{k}$) or the dependencies ($m_{k}$).

A _main effect structure_ can be imposed on the profiles scores, such that

$v_{ks}=\sum_{r}z_{kr}b_{rs}^{(z)},$ (1)

or a _main effect and two-way association structure_ where

$v_{ks}=\sum_{r}z_{kr}b_{rs}^{(z)}+\sum_{r,r^{\prime}>r}z_{kr}z_{kr^{\prime}}b_{rr^{\prime}s}^{(z)}=\sum_{r}z_{kr}\left(b_{rs}^{(z)}+\sum_{r^{\prime}>r}z_{kr^{\prime}}b_{rr^{\prime}s}^{(z)}\right).$
(2)

There are also other possibilities, such as all main effects and a few two-way
associations or extending to higher order associations. We therefore define a
$K\times Q$ design matrix $\bm{Z}$ that includes as columns the relevant
effects, such that we can write $\bm{V}=\bm{Z}\bm{B}_{z}$, where the weights
$\bm{B}_{z}$ are parameters to be estimated.

In a similar fashion, we can construct a design matrix $\bm{W}$ that imposes a
structure on the dependencies, that is

$\bm{m}=\bm{W}\bm{b}_{w}.$

Again, the columns of $\bm{W}$ are defined by the variables $z_{kr}$ and
possibly products of these for pairwise or higher order associations. To
obtain a hierarchical model, the columns of $\bm{Z}$ should be a subset of the
columns of $\bm{W}$. The matrix $\bm{W}$ is of size $K\times T$.

Summarizing, the extended stereotype model (ESM) can be written as

$\theta_{ik}=\bm{w}^{\prime}_{k}\bm{b}_{w}+\bm{x}^{\prime}_{i}\bm{B}_{x}\bm{B}^{\prime}_{z}\bm{z}_{k},$

where $\bm{w}_{k}$, $\bm{x}_{i}$, and $\bm{z}_{k}$ are known and can be
designed to answer certain research questions. The parameters of the model are
$\bm{b}_{w}$, $\bm{B}_{x}$ and $\bm{B}_{z}$.

###  3.2 Modelling the log odds

Let us consider the following two profiles $k=\bullet\bullet 1\bullet\bullet$
and $l=\bullet\bullet 0\bullet\bullet$, that only differ in the $r$-th
position, that is, the dots indicate that the corresponding responses in the
profile can take any value as long as they are the same in the two profiles.

In the Extended Stereotype Model the following expression pertains to this log
odds

$\lambda_{r}(\bm{x}_{i})=\log\left(\frac{\pi_{k}(\bm{x}_{i})}{\pi_{l}(\bm{x}_{i})}\right)=\log(\pi_{k}(\bm{x}_{i}))-\log(\pi_{l}(\bm{x}_{i}))=m_{k}-m_{l}+\bm{x}^{\prime}_{i}\bm{B}_{x}(\bm{v}_{k}-\bm{v}_{l}).$

The log odds can be further worked out, but the expression depend on the
structural form of the profile scores.

  * •

With a _main effect structure_ on the profile scores, we have for the last
term

$(\bm{v}_{k}-\bm{v}_{l})=\bm{B}^{\prime}_{z}(\bm{z}_{k}-\bm{z}_{l})=(z_{kr}-z_{lr})\sum_{s}b^{(z)}_{rs},$

because only the $z$-values of the $r$-th response variable differ between the
two profiles, the others are equal (see equation 1). Furthermore,
$z_{kr}-z_{lr}=1$. Therefore, when $x_{p}$ increases with a unit,
$\lambda_{r}(\bm{x}_{i})$ goes up by $\sum_{s}b_{ps}^{(x)}b_{rs}^{(z)}$.

  * •

With a _main and association effect structure_ on the profile scores, we have
for the last term (see equation 2)

$(\bm{v}_{k}-\bm{v}_{l})=(z_{kr}-z_{lr})\sum_{s}\left(b_{rs}^{(z)}+\sum_{r^{\prime}>r}z_{kr^{\prime}}b_{rr^{\prime}s}^{(z)}\right),$

which does not only depend on the main effect parameters anymore. Therefore,
the change in log-odds has with a unit increase in a predictor has no simple
expression assuming this structure and consequently, model interpretation
should focus on the highest order effects.

###  3.3 Modelling the log odds ratio

Let us consider the following four profiles $k=\bullet\bullet 11\bullet$,
$l=\bullet\bullet 01\bullet$, $n=\bullet\bullet 10\bullet$, and
$o=\bullet\bullet 00\bullet$, that only differ in the $r$ and $r^{\prime}$-th
position, that is, the dots indicate that the corresponding responses in the
profile can take any value as long as they are the same in the four profiles.

In the Extended Stereotype Model the following expression pertains to this log
odds ratio

$\displaystyle\theta_{r,r^{\prime}}(\bm{x}_{i})$ $\displaystyle=$
$\displaystyle\log\left(\frac{\pi_{k}(\bm{x}_{i})\pi_{o}(\bm{x}_{i})}{\pi_{l}(\bm{x}_{i})\pi_{n}(\bm{x}_{i})}\right)=\log(\pi_{k}(\bm{x}_{i}))+\log(\pi_{o}(\bm{x}_{i}))-\log(\pi_{l}(\bm{x}_{i}))-\log(\pi_{n}(\bm{x}_{i}))$
$\displaystyle=$ $\displaystyle
m_{k}+m_{o}-m_{l}-m_{n}+\bm{x}^{\prime}_{i}\bm{B}_{x}(\bm{v}_{k}+\bm{v}_{o}-\bm{v}_{l}-\bm{v}_{n}).$

Again, the expression depends on the imposed structure on the profile scores.

  * •

With a _main effect structure_ on the profile scores, all
$\theta_{r,r^{\prime}}(\bm{x}_{i})=m_{k}+m_{o}-m_{l}-m_{n}$, that is the log
odds ratio does not depend on the predictors.

  * •

With a _main and association effect structure_ we obtain

$(\bm{v}_{k}+\bm{v}_{o}-\bm{v}_{l}-\bm{v}_{n})=\bm{B}_{z}^{\prime}(\bm{z}_{k}+\bm{z}_{o}-\bm{z}_{l}-\bm{z}_{n})=(z_{kr}z_{kr^{\prime}}+z_{or}z_{or^{\prime}}-z_{lr}z_{lr^{\prime}}-z_{nr}z_{nr^{\prime}})\sum_{s}b_{rr^{\prime}s}^{(z)}=\sum_{s}b_{rr^{\prime}s}^{(z)},$

so that, when $x_{p}$ increases with a unit, the log odds ratio
$\theta_{r,r^{\prime}}(\bm{x}_{i})$ goes up by
$\sum_{s}b_{ps}^{(x)}b_{rr^{\prime}s}^{(z)}$.

###  3.4 Relationship with loglinear models and logistic regression

When all predictor variables are categorical a loglinear analysis would be
feasible. In a loglinear analysis, we can specify for each (sub)set of
variables whether they are associated or not. Often, loglinear models are
constrained to be hierarchical, in the sense that when an association is
included in the model all lower order associations between and main effects of
the variables are also included.

Our extended stereotype model is equivalent to certain loglinear models.
First, in the extended stereotype model, we do not make any assumption about
the relationships among the predictor variables, so in the corresponding
loglinear model we need all associations among the predictor variables.
Second, the associations among the response variables are defined through the
design matrix $\bm{W}$, so the corresponding loglinear model needs to impose
the same structure on the responses. Finally, the core of the extended
stereotype model is in the relationships between predictors and responses, and
is captured by the coefficients $\bm{B}_{x}$ and $\bm{B}_{z}$. Note that in
our model we do not have the possibility to define precisely which predictor
variable has an influence on which outcome variable. When a predictor variable
is included in the model, it will have an effect on all responses as defined
in $\bm{Z}$. Therefore, not all loglinear models can be formulated as an
extended stereotype model. In the next Section, we will delve further into
these relationships.

Another difference is that the loglinear model operates with so-called grouped
data, whereas our extended stereotype model works with ungrouped data (see
Agresti,, 2013, Section 5.2.3). The deviance for those two forms of data is
different, even when they represent the same information. For grouped data,
the deviance can be used as a model fit statistic, which is not the case for
ungrouped data. We need to take this into account when comparing the results
of loglinear models with the results of the extended stereotype model.
Differences between deviances of two models represent statistical tests of
model parameters. These differences between deviances are the same for the
loglinear model and the extended stereotype model.

When there is only a single response variable, the extended stereotype model
is equal to a logistic regression model with the same predictors. Therefore,
both logistic regression as well as loglinear modelling are special cases of
our extended stereotype model. Of course, if the data allows for direct
loglinear or logistic regression modelling, a researcher better uses specific
software designed for that purpose.

##  4 Maximum Likelihood Estimation and Model Selection

###  4.1 MM algorithm

In this section, we present an MM algorithm (Heiser,, 1995; Hunter and Lange,,
2004; Nguyen,, 2017) for the estimation of the model parameters. The idea of
MM for finding a minimum of the function $\mathcal{L}(\bm{\theta})$, where
$\bm{\theta}$ is a vector of parameters, is to define an auxiliary function,
called a _majorization function_ , $\mathcal{M}(\bm{\theta}|\bm{\vartheta})$
with two characteristics

$\mathcal{L}(\bm{\vartheta})=\mathcal{M}(\bm{\vartheta}|\bm{\vartheta})\\\ $

where $\bm{\vartheta}$ is a supporting point, and

$\mathcal{L}(\bm{\theta})\leq\mathcal{M}(\bm{\theta}|\bm{\vartheta}).$

The two equations tell us that $\mathcal{M}(\bm{\theta}|\bm{\vartheta})$ is a
function that lies above (i.e., majorizes) the original function and touches
the original function at the support point. The support point is defined by
the current parameter estimates.

A convergent algorithm can be constructed because

$\mathcal{L}(\bm{\theta}^{+})\leq\mathcal{M}(\bm{\theta}^{+}|\bm{\vartheta})\leq\mathcal{M}(\bm{\vartheta}|\bm{\vartheta})=\mathcal{L}(\bm{\vartheta}),$

where $\bm{\theta}^{+}$ is

$\bm{\theta}^{+}=\mathrm{argmin}_{\bm{\theta}}\
\mathcal{M}(\bm{\theta}|\bm{\vartheta}),$

the updated parameters. An advantage of MM algorithms is that they always
converge monotonically to a (local) minimum. The challenge is to find a
parametrized function family, $\mathcal{M}(\bm{\theta}|\bm{\vartheta})$, that
can be used in every step.

For our multinomial models, we like to minimize the negative log-likelihood,
that is

$\mathcal{L}(\bm{\theta})=\sum_{i}\mathcal{L}_{i}(\bm{\theta}_{i})=-\sum_{i}\sum_{k}g_{ik}\log\left[\frac{\exp(\theta_{ik})}{\sum_{k^{\prime}=1}^{K}\exp(\theta_{ik^{\prime}})}\right],$

where $\bm{\theta}_{i}=[\theta_{i1},\ldots,\theta_{iK}]^{\prime}$ and $g_{ik}$
are the elements of the response profile indicator matrix ($\mathbf{G}$), rows
of which are indicated by $\bm{g}_{i}$. For finding a majorization function we
can work with $\mathcal{L}_{i}(\bm{\theta}_{i})$ and afterwards sum the
contributions.

We majorize the negative log-likelihood using a quadratic upper bound (Hunter
and Lange,, 2004, Section 3.4) which is based on the following inequality

$\mathcal{L}_{i}(\bm{\theta}_{i})\leq\mathcal{L}_{i}(\bm{\vartheta}_{i})+\frac{\partial\mathcal{L}_{i}(\bm{\theta}_{i})}{\partial\bm{\theta}_{i}}(\bm{\theta}_{i}-\bm{\vartheta}_{i})+\frac{1}{2}(\bm{\theta}_{i}-\bm{\vartheta}_{i})^{\prime}\bm{\Omega}(\bm{\theta}_{i}-\bm{\vartheta}_{i}),$

for any matrix $\bm{\Omega}$ such that
$\bm{\Omega}-\partial^{2}\mathcal{L}_{i}(\bm{\theta}_{i})$ is positive semi
definite. Evans, (2014) showed that the matrix
$\bm{\Omega}=\frac{1}{4}\mathbf{I}$ defines such a matrix. The first
derivative is of $\mathcal{L}_{i}(\bm{\theta}_{i})$ is

$\frac{\partial\mathcal{L}_{i}(\bm{\theta}_{i})}{\partial\bm{\theta}_{i}}=-(\mathbf{g}_{i}-\bm{\pi}_{i}),$

such that

$\mathcal{L}_{i}(\bm{\theta}_{i})\leq\mathcal{L}_{i}(\bm{\vartheta}_{i})-(\bm{\theta}_{i}-\bm{\vartheta}_{i})^{\prime}(\mathbf{g}_{i}-\bm{\pi}_{i})+\frac{1}{8}(\bm{\theta}_{i}-\bm{\vartheta}_{i})^{\prime}(\bm{\theta}_{i}-\bm{\vartheta}_{i}).$

Working out this expression and defining
$\bm{h}_{i}=\bm{\vartheta}_{i}+4(\mathbf{g}_{i}-\bm{\pi}_{i})$, we obtain
after some simplification

$\mathcal{L}_{i}(\bm{\theta}_{i})\leq\frac{1}{8}\bm{\theta}_{i}^{\prime}\bm{\theta}_{i}-\frac{1}{8}\cdot
2\bm{\theta}_{i}^{\prime}\bm{h}_{i}+\mathrm{constant}$

(see De Leeuw,, 2006; Groenen et al.,, 2003). Note that $\bm{h}_{i}$ is
independent of $\bm{\theta}_{i}$, so we can add
$\frac{1}{8}\bm{h}_{i}^{\prime}\bm{h}_{i}$ to obtain

$\mathcal{L}_{i}(\bm{\theta}_{i})\leq\frac{1}{8}\bm{\theta}_{i}^{\prime}\bm{\theta}_{i}-\frac{1}{8}\cdot
2\bm{\theta}_{i}^{\prime}\bm{h}_{i}+\frac{1}{8}\bm{h}_{i}^{\prime}\bm{h}_{i}+\mathrm{constant},$

which can be rewritten as

$\mathcal{L}_{i}(\bm{\theta}_{i})\leq\frac{1}{8}||\bm{h}_{i}-\bm{\theta}_{i}||^{2}+\mathrm{constant}=\mathcal{M}_{i}(\bm{\theta}_{i}|\bm{h}_{i})+\mathrm{constant}.$

By using the summation property, we therefore have the following majorization
function

$\mathcal{M}(\bm{\theta}|\bm{\delta})=\sum_{i}\mathcal{M}_{i}(\bm{\theta}_{i}|\bm{h}_{i})=\sum_{i}\sum_{k}(h_{ik}-\theta_{ik})^{2}=\|\bm{H}-\bm{\Theta}\|^{2},$

where $\bm{H}=\\{h_{ik}\\}$ with $h_{ik}=\vartheta_{ik}+4(g_{ik}-\pi_{ik})$,
so called _working responses_. The matrix $\bm{\Theta}$ has elements
$\theta_{ik}$ and is parametrized as

$\bm{\Theta}=\bm{1m}^{\prime}+\bm{XB}_{x}\bm{B}^{\prime}_{z}\bm{Z}^{\prime}.$

Within each iteration, we have to minimize

$\|\bm{H}-\bm{1m}^{\prime}-\bm{XB}_{x}\bm{B}^{\prime}_{z}\bm{Z}^{\prime}\|^{2},$

which we do in an alternating fashion. The update for $\bm{b}_{w}$ we first
define the auxiliary unconstrained vector

$\displaystyle\tilde{\bm{m}}=\frac{1}{N}(\bm{H}-\bm{XB}_{x}\bm{B}^{\prime}_{z}\bm{Z}^{\prime})^{\prime}\bm{1},$

from which we find the update as

$\displaystyle\bm{b}_{w}^{+}=\left(\bm{W}^{\prime}\bm{W}\right)^{-1}\bm{W}^{\prime}\tilde{\bm{m}},$
(3)

and subsequently compute $\bm{m}=\bm{W}\bm{b}_{w}$.

To find updates for $\bm{B}_{x}$ and $\bm{B}_{z}$, we first compute

$\bm{H}_{c}=\bm{H}-\bm{1m}^{\prime}.$

to obtain the loss function

$\|\bm{H}_{c}-\bm{XB}_{x}\bm{B}^{\prime}_{z}\bm{Z}^{\prime}\|^{2},$

which is the so-called CANDELINC problem (Carroll et al.,, 1980). Updates of
$\bm{B}_{x}$ and $\bm{B}_{z}$ can be obtained from a generalized singular
value decomposition of the matrix

$(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{H}_{c}\bm{Z}(\bm{Z}^{\prime}\bm{Z})^{-1}$

in the metrics $(\bm{X}^{\prime}\bm{X})$ and $(\bm{Z}^{\prime}\bm{Z})$
(Takane,, 2013). Therefore, take the SVD of

$(\bm{X}^{\prime}\bm{X})^{-\frac{1}{2}}\bm{X}^{\prime}\bm{H}_{c}\bm{Z}(\bm{Z}^{\prime}\bm{Z})^{-\frac{1}{2}}=\bm{P}\bm{\Phi}\bm{Q}^{\prime}$

and define the updates

$\displaystyle\bm{B}_{x}^{+}$ $\displaystyle=$
$\displaystyle\sqrt{N}(\bm{X}^{\prime}\bm{X})^{-\frac{1}{2}}\bm{P}_{S}$ (4)
$\displaystyle\bm{B}_{z}^{+}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}(\bm{Z}^{\prime}\bm{Z})^{-\frac{1}{2}}\bm{Q}_{S}\bm{\Phi}_{S},$
(5)

where $\bm{P}_{S}$ and $\bm{Q}_{S}$ are the first $S$ left and right singular
vectors, and $\bm{\Phi}_{S}$ is the $S\times S$ diagonal matrix with the $S$
largest singular values.

In short, in every iteration of the MM algorithm, we compute the matrix of
working responses $\bm{H}$. Then update $\bm{b}_{w}$ by Equation 3 and update
$\bm{B}_{x}$ and $\bm{B}_{z}$ by Equations 4 and 5. Having these updates, new
working responses are computed. This sequence of steps monotonically minimizes
the negative log-likelihood to a global minimum.

###  4.2 Software

We developed R-functions for estimation of the model parameters. The functions
are made public on the github-page of the author.

###  4.3 Model Selection

For selecting a model there are four issues:

  * •

Selecting the dimensionality/rank $S$;

  * •

Selecting the structure on the profile scores ($\bm{v}_{k}$), that is,
determine $\bm{Z}$;

  * •

Selecting the linear predictor of the object scores ($\bm{u}_{i}$), that is,
determine $\bm{X}$;

  * •

Determining the residual association structure, that is, determine $\bm{W}$.

As we have a maximum likelihood method, we can rely on information criteria
and likelihood ratio statistics. Another approach would be to compute the
standard errors. For maximum likelihood estimation, standard errors are given
by the Hessian matrix. However, as we use an MM algorithm that does not
optimize the likelihood itself this Hessian is not a by-product of the
algorithm. This need not be a problem, confidence intervals for the parameters
can be obtained from a bootstrap procedure (Efron,, 1979; Efron and
Tibshirani,, 1986).

As there are many choices to make, we suggest to use a step wise approach.
Starting with a full matrices $\bm{Z}$, $\bm{X}$, and $\bm{W}$ first select
the dimensionality $S$. Second, using the just selected $S$, find the optimal
profile structure (i.e., find the optimal $\bm{Z}$). With the optimal $S$ and
profile structure, find which predictor variables contribute to the prediction
of the responses (i.e., find the optimal $\bm{X}$). Finally, determine which
residual associations are necessary (i.e., find the optimal $\bm{W}$).

For selecting the dimensionality, the likelihood ratio statistic does not
follow a chi-square distribution (Takane et al.,, 2003). Therefore, we use
Akaike’s Information Criterion (AIC) for selecting an optimal $S$. Having an
optimal dimensionality, for selecting the appropriate structure on the profile
scores the likelihood ratio test or an information criterion can be used. For
selecting the appropriate structure on the object scores, again likelihood
ratio statistics or an information criterion can be used. Finally, we have to
select the residual associations, that is, which effects are presented in
$\bm{W}$ taking into account that the matrix $\bm{Z}$ should represent a
subset. We favor the use of AICs for all these tasks, as it makes less
assumptions

For the AIC, or another information criterion, the degrees of freedom are
needed. For reduced rank models, Mukherjee et al., (2015) argued that the
naive estimator of the degrees of freedom is not optimal and derived another
unbiased estimator. This unbiased estimator, however, is only available for
Gaussian reduced rank models. Therefore, we will use the naive one, where
$df=T+S(P+Q-S)$.

##  5 Empirical Examples

In this section, we will show three examples. The first analysis pertains to
alcohol, cigarette and marijuana use among high school students. The analysis
focuses on the relationship between our extended stereotype model and
loglinear analysis. The second analysis pertains to respiratory problems among
coal miners. We include this analysis to contrast our approach with marginal
models, and reiterate the conclusion from McCullagh and Nelder, (1989) about
the differences between marginal and conditional models. The third example
pertains to the analysis of depressive and anxiety disorders and discusses
model selection and interpretation for a data set having five response
variables and eight predictors.

###  5.1 Alcohol, cigarette and marijuana use among high school seniors

Agresti (2013, Table 10.1) shows data of 2276 senior high school students and
their use of alcohol (A), cigarettes (C), and marijuana (M), all binary yes or
no. Also available are two predictor variables, race (R; white/other) and
gender (G; female/male). This data set with three responses and two predictor
variables has been analyzed by Agresti using loglinear models. In this
Section, we re-analyze the data using our extended stereotype models. We would
like to show relationships between the two approaches. Table 1 shows the main
results.

|  | LLA | ESM
---|---|---|---
| Model | $G^{2}$ | $df$ | Deviance | npar
1 | Mutual independende + GR | 1325.14 | 25 | 7900.19 | 3
2 | Homogeneous Association | 15.34 | 16 | 6590.38 | 12
3 | All three-factor terms | 5.27 | 6 | - | -
4a | (2) - AC | 201.20 | 17 | 6776.24 | 11
4b | (2) - AM | 106.96 | 17 | 6681.99 | 11
4c | (2) - CM | 513.47 | 17 | 7088.51 | 11
4d | (2) - AG | 18.72 | 17 | - | -
4e | (2) - AR | 20.32 | 17 | - | -
4f | (2) - CG | 16.32 | 17 | - | -
4g | (2) - CR | 15.78 | 17 | - | -
4h | (2) - GM | 25.16 | 17 | - | -
4i | (2) - GR | 18.93 | 17 | - | -
5 | (AC, AN, CM, AG, AR, GM, MR, GR) | 16.74 | 18 | 6601.55 | 10
6 | (AC, AN, CM, AG, AR, GM, GR) | 19.91 | 19 | - | -
7 | (AC, AN, CM, AG, AR, GR) | 28.81 | 20 | 6603.84 | 8
Table 1: Results from loglinear analyses (LLA) and the extended stereotype
model (ESM) for the alcohol, cigarettes, and marijuana data.

The left column of Table 1 shows all loglinear models as fitted by Agresti.
The middle block shows $G^{2}$-statistics and degrees of freedom as reported
in Table 10.2 of Agresti. Note that the contingency table has 32 cells, the
number of parameters of the models is therefore, $32-df$. The last two columns
show the results of corresponding extended stereotype models. When no numbers
are presented there is no equivalent extended stereotype model.

Model 1 can be estimated through a ESM by setting $S=0$ and defining $\bm{W}$
to only include the main effects. The number of parameters is 3, the length of
$\bm{b}_{w}$. Note that the four parameters representing the association
structure among the predictors are not counted as such in the ESM. In the
loglinear model there are 7 parameters (32 - 25).

Model 2, the model with all pairwise associations among the five variables can
be obtained as an ESM. The definition of $\bm{W}$ should exclude the three-
variable association among the three response variables, but include all three
pairwise associations and the three main effects. The design matrix $\bm{Z}$
should include only the main effects of the three responses and $\bm{X}$
should have only the two main effects (not an interaction between race and
gender). The number of fitted parameters is 12, nine more than in Model 1. The
same difference can be found in the $df$ of the two models for the loglinear
analysis. Note also that 1325.14 - 15.34 = 7900.19 - 5590.38 = 1309.8.

There is no extended stereotype equivalent of the all three-factor terms
loglinear model. In the extended stereotype models, model complexity and
effects can be defined through the three design matrices $\bm{X}$, $\bm{Z}$,
and $\bm{W}$. We need to include an interaction term between race and gender
in the matrix $\bm{X}$ to model the three variable associations RGA, RGC, and
RGM. For the three variable associations between RAC, RAM, RCM, and GAC, GAM,
GCM, however, we need to define the main effect and pairwise associations in
$\bm{Z}$. When including an interaction effect through $\bm{X}$ and pairwise
associations in $\bm{Z}$, the four way associations RGAC, RGCM, and RGAM are
automatically included in the model. Similar arguments play a role, for all
other LLA models that do not have an ESM equivalent.

Differences in $G^{2}$ for the loglinear models equal differences in deviances
for the extended stereotype models. Differences in degrees of freedom match
differences in number of parameters. So, tests of specific (sets of) effects
through the likelihood ratio statistic are the same for the two approaches.

###  5.2 Analysis of Coalminers data

McCullagh and Nelder, (1989) show data collected by Ashford and Sowden in 1970
about respiratory problems among coal miners. Each participant was classified
as either suffering from Breathlessnees ($B$) or not and suffering from Wheeze
($W$) or not. Also available is the age of the participant.

As our Extended Stereotype Model needs person level data, we transformed the
data in two matrices $\bm{X}$ and $\bm{Y}$ each having 18282 rows. Following
McCullagh and Nelder (pp 229 - 235), we defined $X_{1}=(age-42)/5$. The matrix
$\bm{Y}$ has two columns corresponding to Breathlessness and Wheeze. The
corresponding profile matrix $\bm{G}$ has four columns.

For this data, we would like to have a main and association effect structure
on the response profiles and therefore define the matrix $\bm{Z}$ with four
rows ($K=4$) and three columns ($Q=3$) as

Profile | B | W | B:W
---|---|---|---
00 | $-\frac{1}{2}$ | $-\frac{1}{2}$ | $\frac{1}{4}$
01 | $-\frac{1}{2}$ | $\frac{1}{2}$ | $-\frac{1}{4}$
10 | $\frac{1}{2}$ | $-\frac{1}{2}$ | $-\frac{1}{4}$
11 | $\frac{1}{2}$ | $\frac{1}{2}$ | $\frac{1}{4}$

.

A model in maximum dimensionality ($S=1$) is fitted. The estimated
coefficients are shown in Table 2. The effect of age is substantial, with
coefficients 0.43, 0.12, and -0.16. As discussed in Section 3, interpretation
should focus on the highest order terms. Therefore, with every unit increase
in age the log odds ratio between breathlessness and wheeze goes down by 0.16.

| ESM | MLM
---|---|---
| B | W | B:W | B | W | B:W
1 | -1.94 | -0.49 | 3.03 | -2.26 | -1.49 | 3.02
age | 0.43 | 0.12 | -0.16 | 0.52 | 0.33 | -0.13
Table 2: Fitted coefficients for the coalminer data set. In the rows we have
the intercept, the effect of age, and the effect of age squared (age2). In the
columns, we have the effect of those predictors on the log-odds of
Breathlessness (B), Wheeze (W), and their interaction (B:W). First three
columns are coefficients of the extended stereotype model (ESM), last three
columns of the marginal logit model (MLM) as fitted by McCullagh, (1980).

McCullagh and Nelder, (1989) analysed the data using a marginal logit model.
The coefficients of their model are shown in the last three columns of Table 2
and show that the coefficients for breathlessness and the association are
quite similar, the coefficients for wheeze, however, are quit a bit smaller.
Where the coefficients of McCullagh and Nelder, (1989) refer to _marginal
effects_ our coefficients refer to _average conditional effects_. It is well
known that these effects can differ. In this case, the effects on
Breathlessness and the association are similar. The effect on Wheeze, however,
is much smaller as was also pointed out by McCullagh and Nelder, (1989).

###  5.3 Depressive and Anxiety disorders

Here, we will analyse a data set from the Netherlands Study on Depression and
Anxiety (NESDA; Penninx et al.,, 2008), where we have for 786 participants
data available on 5 psychiatric diagnosis (dysthymia (D), major depressive
disorder (M), generalized anxiety disorder (G), social phobia (S), and panic
disorder (P)). Psychiatric disorders commonly co-occur, a phenomenon called
comorbidity. In this particular data set, 272 participants have one disorder,
235 have two disorder, 147 have three, 96 have four, and 36 participants have
all five disorders. Research focusses on the link between personality
characteristics and these psychiatric disorders (Spinhoven et al.,, 2009).
Measurements of the big-five personality characteristics neuroticism (N),
extraversion (E), openness to experience (O), agreeableness (A), and
conscientiousness (C) are available and will be used as predictors. Besides
those five personality characteristics, sex, age and education are available
as control variables, which will be included in all models. The variable sex
is dichotomous with code 1 for females, the other predictor variables are
centered and scaled to have standard deviation equal to one. With five
response variables, $2^{5}=32$ response profiles are possible of which 28
occur in the data.

The results of the model selection procedure are shown in Table 3. Using a
profile structure of the order 3, that is all main effects plus all two- and
three-way associations, we first selected the optimal dimensionality ($S$). It
can be seen that the two-dimensional solution has the lowest AIC. Next, we
selected the profile structure ($\bm{Z}$), where we first compare the model
with all main effects and two- and three-way association to the model with the
model with all main effects and two-way associations, and the model with all
main effects. The latter has the lowest AIC. Than we checked whether any of
the two-way associations improves model fit over this order 1 (i.e., main
effects) model. Adding the two-way association between dysthymia and
generalized anxiety disorders and between major depressive disorder and
generalized anxiety disorder resulted in a lower AIC. Adding both, however,
lead to an increase in AIC, so we only include the association between
dysthymia and generalized anxiety disorders. Next, we checked whether a
personality characteristic ($\bm{X}$) can be left out of the model and found
that both Openness and Agreeableness can be left out of the model. In the
final step, we verified which residual associations are needed ($\bm{W}$) ,
where we see that all four and three-way associations do not contribute
significantly to the model fit.

$S$ | $Z$ | $X$ | $W$ | AIC
---|---|---|---|---
5 | 3 | N + E + O + A + C | 5 | 4449.2
4 | 3 | N + E + O + A + C | 5 | 4421.8
3 | 3 | N + E + O + A + C | 5 | 4404.6
2 | 3 | N + E + O + A + C | 5 | _4392.2_
1 | 3 | N + E + O + A + C | 5 | 4392.5
2 | 2 | N + E + O + A + C | 5 | 4379.0
2 | 1 | N + E + O + A + C | 5 | _4363.3_
2 |  1 \+ D:M | N + E + O + A + C | 5 | 4366.7
2 |  1 \+ D:G | N + E + O + A + C | 5 | _4361.7_
2 |  1 \+ D:S | N + E + O + A + C | 5 | 4366.9
2 |  1 \+ D:P | N + E + O + A + C | 5 | 4364.7
2 |  1 \+ M:G | N + E + O + A + C | 5 | 4362.5
2 |  1 \+ M:S | N + E + O + A + C | 5 | 4365.2
2 |  1 \+ M:P | N + E + O + A + C | 5 | 4367.1
2 |  1 \+ G:S | N + E + O + A + C | 5 | 4366.4
2 |  1 \+ G:P | N + E + O + A + C | 5 | 4367.1
2 |  1 \+ S:P | N + E + O + A + C | 5 | 4365.2
2 |  1 \+ D:G | E + O + A + C | 5 | 4430.2
2 |  1 \+ D:G | N + O + A + C | 5 | 4378.7
2 |  1 \+ D:G | N + E + A + C | 5 | _4358.0_
2 |  1 \+ D:G | N + E + O + C | 5 | 4359.2
2 |  1 \+ D:G | N + E + O + A | 5 | 4363.9
2 |  1 \+ D:G | N + E + C | 5 | _4355.7_
2 |  1 \+ D:G | N + E + C | 4 | 4353.7
2 |  1 \+ D:G | N + E + C | 3 | 4347.0
2 |  1 \+ D:G | N + E + C | 2 | 4336.0
Table 3: Model Selection for Nesda data. $S$ indicates the dimensionality or
rank. The column labeled $Z$ refers to the structure on the profiles scores,
where a 3 indicates all main effects plus all two- and three-way associations,
2 indicates all main effects and all two-way associations, and 1 indicates all
main effects. 1 + D:M indicates all main effects plus the two-way association
between D and M (Dysthymia and Major depressive disorder). Under $X$ it is
indicated which predictor variables are in the model. The column labeled $W$
refers to the structure on the associations among the response variables,
where the same notation is used as for the profile scores (column $Z$). The
final column gives the AIC of the model.

The coefficients of the final model are given in Table 4. We see, for example,
that with every standard deviation increase of Neuroticism the log-odds of
Major depressive disorder goes up by 0.53, the log-odds of social phobia goes
up by 0.55, and the log-odds of panic disorder goes up by 0.17. The log-odds
ratio between dysthymia and generalized anxiety disorder goes down by 0.38
with every standard deviation increase in Neuroticism. Similar interpretations
can be given to the other coefficients.

| D | M | G | S | P | D:G
---|---|---|---|---|---|---
1 | -1.61 | 1.19 | -0.82 | -0.48 | 0.03 | 1.05
Gender | -0.07 | 0.02 | -0.01 | -0.04 | 0.20 | -0.02
Age | 0.04 | 0.05 | 0.02 | 0.05 | 0.00 | -0.03
Edu | -0.04 | -0.22 | -0.05 | -0.16 | -0.32 | 0.17
N | 0.36 | 0.53 | 0.18 | 0.55 | 0.17 | -0.38
E | -0.25 | -0.25 | -0.10 | -0.31 | 0.11 | 0.16
C | -0.10 | 0.01 | -0.01 | -0.06 | 0.24 | -0.02
Table 4: Estimated regression coefficients for the NESDA data. The
coefficients for the intercept are given by $\bm{b}_{w}$, the others by
$\bm{A}=\bm{B}_{x}\bm{B}^{\prime}_{z}$. Abbreviations: D = dysthymia, M =
major depressive disorder, G = generalized anxiety disorder, S = social
phobia, P = panic disorder, Edu = number of years of education, N =
Neuroticism, E= Extraversion, C = Conscientiousness.

The estimated pairwise residual associations among the response variables are
1.01 for DG, -0.01 for DS, 0.33 for DP, 0.61 for MG, -0.14 for MS, 0.16 for
MP, 0.45 for GS, 0.62 for GP, and 1.23 for SP.

##  6 Discussion

In this paper, we developed a new regression method for multivariate binary
outcome variables. The methods can deal with both numerical and categorical
predictors, as opposed to traditional loglinear models that are only defined
for categorical variables. We focused on binary response variables, but the
model can, in theory, be extended to other types of categorical response
variables. The method is a full maximum likelihood method, assuming a
multinomial distribution for the response profiles. In the model, we can
specify effects of predictor variables on the response variables, but also on
associations between response variables. In the first case, the log odds is a
(linear) function of the predictor variables, in the latter case, the log odds
ratio between two response variables is a (linear) function of one or more
predictor variables. A researcher can specify, by appropriately defining the
design matrices $\bm{X}$ and $\bm{Z}$, on which effects to focus. A limitation
of the model is that it is not possible to define different linear predictors
for the different responses. More specifically, if $\bm{X}$ includes two
predictors, than these two predictors have an effect on each main effect and
every association as defined by $\bm{Z}$. We elaborated on this issue with the
first empirical example (Section 5.1).

We developed an MM algorithm for maximum likelihood estimation of the model
parameters. In MM algorithms a majorization function is defined in every
iteration which is subsequently minimized. MM algorithms often require many
iterations, but when the majorization function is suitable chosen the
minimization of this function is easy, as in our case where in every iteration
a singular value decomposition has to be computed of a (small) matrix.

The extended stereotype model is a so-called conditional model, that is, the
effects on response variable $r$ should be interpreted conditional on the
other response variables. This type of model can be contrasted with marginal
models (McCullagh and Nelder,, 1989; Liang et al.,, 1992; Glonek and
McCullagh,, 1995) where the effects of predictors on the responses are
marginal effects. We showed the differences between the conditional and
marginal approach in our second empirical example (Section 5.2). Estimation of
marginal models by maximum likelihood is, however, awkward ((Agresti,, 2013).
With only categorical variables an algorithm has been proposed by Bergsma et
al., (2009). With numerical predictors and having more than two response
variables researchers commonly use generalized estimating equations (Liang and
Zeger,, 1986; Liang et al.,, 1992; Hardin and Hilbe,, 2002). A marginal
approach can also be specified by reduced rank logistic models (Yee and
Hastie,, 2003; De Rooij,, 2022), where conditional on the reduced rank the
responses are assumed to be independent. Such an approach, however, does not
allow to also model the relationships between responses.

For model selection, likelihood ratio tests and information criteria can be
used. In the example in Section 5.3, we used Akaike’s Information criterion.
We used a naive approach for determining the effective degrees of freedom. For
Gaussian reduced rank models there are unbiased estimators for the effective
degrees of freedom (Mukherjee et al.,, 2015), but these have not been
generalized for other reduced rank models. Mukherjee et al., (2015) point out
that having a better estimator of the effective degrees of freedom influences
the final model choice. Furtehr research is needed. We suggested a step wise
procedure for model selection, starting with the dimensionality, then the
selection of the profile score structure, followed by the selection of the
object score structure (which predictor variables), and finally determining
the residual association structure. The step wise approach might not be
optimal in all circumstances, further research is needed.

## Statements and Declarations

#### Data availability:

This paper uses three empirical data sets as illustrations.

  1. 1.

The alcohol, cigarette, and marijuana data are available in the R-package
CatDataAnalysis

  2. 2.

The coalminers data are available in the R-package VGAM

  3. 3.

The depressive and anxiety disorder data set is available upon request from
the nesda consortium, see https://www.nesda.nl/nesda-english/

R-code for the analyses of these data sets is available on the github-page of
the author.

#### Funding:

The author declares that no funds, grants, or other support were received
during the preparation of this manuscript.

#### Competing Interests:

The authors have no relevant financial or non-financial interests to disclose.

#### Author Contributions:

As this is a single author paper, everything was done by the main author.

## References

  * Afifi and Elashoff, (1969) Afifi, A. A. and Elashoff, R. (1969).  Multivariate two sample tests with dichotomous and continuous variables. i. the location model.  The Annals of Mathematical Statistics, 40(1):290–298.
  * Agresti, (2002) Agresti, A. (2002).  Categorical data analysis.  John Wiley & Sons, second edition.
  * Agresti, (2013) Agresti, A. (2013).  Categorical data analysis.  John Wiley & Sons, third edition.
  * Anderson, (1984) Anderson, J. A. (1984).  Regression and ordered categorical variables.  Journal of the Royal Statistical Society: Series B (Methodological), 46(1):1–22.
  * Ashford and Sowden, (1970) Ashford, J. and Sowden, R. (1970).  Multi-variate probit analysis.  Biometrics, 26:535–546.
  * Bergsma et al., (2009) Bergsma, W., Croon, M., and Hagenaars, J. (2009).  Marginal Models: For Dependent, Clustered, and Longitudinal Categorical Data.  Springer New York.
  * Bishop et al., (1975) Bishop, Y. M., Fienberg, S. E., and Holland, P. W. (1975).  Discrete multivariate analysis: theory and practice.  MIT press.
  * Carroll et al., (1980) Carroll, D. J., Pruzansky, S., and Kruskal, J. B. (1980).  Candelinc: A general approach to multidimensional analysis of many-way arrays with linear constraints on parameters.  Psychometrika, 45:3–24.
  * De Leeuw, (2006) De Leeuw, J. (2006).  Principal component analysis of binary data by iterated singular value decomposition.  Computational statistics & data analysis, 50(1):21–39.
  * De Rooij, (2022) De Rooij, M. (2022).  A new algorithm and a discussion about visualization for logistic reduced rank regression.  Submitted paper.
  * Efron, (1975) Efron, B. (1975).  The efficiency of logistic regression compared to normal discriminant analysis.  Journal of the American Statistical Association, 70(352):892–898.
  * Efron, (1979) Efron, B. (1979).  Bootstrap methods: Another look at the jackknife.  Annals of Statistics, 7:1–26.
  * Efron and Tibshirani, (1986) Efron, B. and Tibshirani, R. (1986).  Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy.  Statistical science, 1:54–75.
  * Evans, (2014) Evans, G. W. (2014).  Logistic Gifi: A Logistic Distance Association Model for Exploratory Analysis of Categorical Data.  PhD thesis, UCLA.
  * Glonek and McCullagh, (1995) Glonek, G. F. and McCullagh, P. (1995).  Multivariate logistic models.  Journal of the Royal Statistical Society: Series B (Methodological), 57(3):533–546.
  * Groenen et al., (2003) Groenen, P. J. F., Giaquinto, P., and Kiers, H. A. L. (2003).  Weighted majorization algorithms for weighted least squares decomposition models.  Econometric Institute Research Papers EI 2003-09, Erasmus University Rotterdam.
  * Hardin and Hilbe, (2002) Hardin, J. W. and Hilbe, J. M. (2002).  Generalized estimating equations.  chapman and hall/CRC.
  * Heiser, (1995) Heiser, W. J. (1995).  Convergent computation by iterative majorization: Theory and applications in multidimensional data analysis.  In Krzanowski, W. J., editor, Recent advances in descriptive multivariate analysis, pages 157–189. Clarendon Press.
  * Hunter and Lange, (2004) Hunter, D. R. and Lange, K. (2004).  A tutorial on MM algorithms.  The American Statistician, 58(1):30–37.
  * Izenman, (1975) Izenman, A. J. (1975).  Reduced-rank regression for the multivariate linear model.  Journal of multivariate analysis, 5(2):248–264.
  * Krzanowski, (1976) Krzanowski, W. (1976).  Canonical representation of the location model for discrimination or classification.  Journal of the American Statistical Association, 71(356):845–848.
  * Krzanowski, (1993) Krzanowski, W. (1993).  The location model for mixtures of categorical and continuous variables.  Journal of Classification, 10:25–49.
  * Krzanowski, (1975) Krzanowski, W. J. (1975).  Discrimination and classification using both binary and continuous variables.  Journal of the American statistical association, 70(352):782–790.
  * Liang and Zeger, (1986) Liang, K.-Y. and Zeger, S. L. (1986).  Longitudinal data analysis using generalized linear models.  Biometrika, 73(1):13–22.
  * Liang et al., (1992) Liang, K.-Y., Zeger, S. L., and Qaqish, B. (1992).  Multivariate regression analyses for categorical data.  Journal of the Royal Statistical Society: Series B (Methodological), 54(1):3–24.
  * McCullagh, (1980) McCullagh, P. (1980).  Regression models for ordinal data.  Journal of the Royal Statistical Society: Series B (Methodological), 42(2):109–127.
  * McCullagh and Nelder, (1989) McCullagh, P. and Nelder, J. A. (1989).  Generalized Linear Models.  Chapman & Hall / CRC.
  * Mukherjee et al., (2015) Mukherjee, A., Chen, K., Wang, N., and Zhu, J. (2015).  On the degrees of freedom of reduced-rank estimators in multivariate regression.  Biometrika, 102(2):457–477.
  * Ng and Jordan, (2001) Ng, A. and Jordan, M. (2001).  On discriminative vs. generative classifiers: A comparison of logistic regression and naive bayes.  Advances in neural information processing systems, 14.
  * Nguyen, (2017) Nguyen, H. D. (2017).  An introduction to majorization-minimization algorithms for machine learning and statistical estimation.  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 7(2):e1198.
  * Olkin and Tate, (1961) Olkin, I. and Tate, R. F. (1961).  Multivariate correlation models with mixed discrete and continuous variables.  The Annals of Mathematical Statistics, pages 448–465.
  * Penninx et al., (2008) Penninx, B. W., Beekman, A. T., Smit, J. H., Zitman, F. G., Nolen, W. A., Spinhoven, P., Cuijpers, P., De Jong, P. J., Van Marwijk, H. W., Assendelft, W. J., et al. (2008).  The Netherlands study of depression and anxiety (NESDA): rationale, objectives and methods.  International journal of methods in psychiatric research, 17(3):121–140.
  * Reinsel et al., (2022) Reinsel, G. C., Velu, R. P., and Chen, K. (2022).  Multivariate Reduced-Rank Regression: Theory, Methods and Applications, volume 225.  Springer Nature.
  * Rubinstein et al., (1997) Rubinstein, Y. D., Hastie, T., et al. (1997).  Discriminative vs informative learning.  In KDD, volume 5, pages 49–53.
  * Spinhoven et al., (2009) Spinhoven, P., De Rooij, M., Heiser, W., Smit, J. H., and Penninx, B. W. (2009).  The role of personality in comorbidity among anxiety and depressive disorders in primary care and specialty care: a cross-sectional analysis.  General hospital psychiatry, 31(5):470–477.
  * Takane, (2013) Takane, Y. (2013).  Constrained principal component analysis and related techniques.  CRC Press.
  * Takane et al., (2003) Takane, Y., van der Heijden, P. G., and Browne, M. W. (2003).  On likelihood ratio tests for dimensionality selection.  In Higuchi, T., Iba, Y., and Ishiguro, M., editors, Proceedings of science of modeling: The 30th anniversary meeting of the information criterion (AIC), pages 348–349. The Institute of Statistical Mathematics Tokyo.
  * Wickens, (1989) Wickens, T. D. (1989).  Multiway contingency tables analysis for the social sciences.  Lawrence Erlbaum Associates, Inc.
  * Yee and Hastie, (2003) Yee, T. W. and Hastie, T. J. (2003).  Reduced-rank vector generalized linear models.  Statistical modelling, 3(1):15–41.


# Epidemic oscillations induced by social network control: the discontinuous
case

Daniele De Martino  Biofisika Institute (CSIC,UPV-EHU) and Ikerbasque Basque
Foundation for Science, Bilbao 48013, Spain  Fabio Caccioli  Department of
Computer Science, University College London, Gower Street WC1E 6EA London (UK)
Systemic Risk Centre, London School of Economics and Political Sciences, WC2A
2AE, London (UK)  London Mathematical Laboratory, 8 Margravine Gardens, London
WC 8RH (UK).

###### Abstract

Epidemic spreading can be suppressed by the introduction of containment
measures such as social distancing and lock downs. Yet, when such measures are
relaxed, new epidemic waves and infection cycles may occur. Here we explore
this issue in compartmentalized epidemic models on graphs in presence of a
feedback between the infection state of the population and the structure of
its social network for the case of discontinuous control. We show that in
random graphs the effect of containment measures is simply captured by a
renormalization of the effective infection rate that accounts for the change
in the branching ratio of the network. In our simple setting, a piece-wise
mean-field approximations can be used to derive analytical formulae for the
number of epidemic waves and their length. A variant of the model with
imperfect information is used to model data of the recent covid-19 epidemics
in the Basque Country and Lombardy, where we estimate the extent of social
network disruption during lock downs and characterize the dynamical
attractors.

## Introduction

The onset of oscillations in a system as a consequence of feedback has been
highlighted since the inception of control theory [1, 2]. Traditionally,
feedback-induced oscillations have been studied in engineering artificial
tools like thermostats and steering devices [3]. More recently, research
focused also on its application to natural systems [4], in particular
homeostasis and its disruption in biological systems, a classical example
being glycemic control and diabetes in human metabolism [5].

Feedback-induced oscillations are currently emerging as governments are trying
to control the evolution of the covid-19 pandemic crisis with containment
measures such as social distancing, lock downs and quarantine.

The modeling of containement measures in compartmentalized epidemic models [6]
is thus under the focus of intense research [7, 8]. It has been very recently
rigorously demonstrated that compartimentalized epidemic models display
oscillations in presence of feedback between infection rate and infection
states [9],and that in general a feedback between order and control parameters
in large interacting systems subject to phase transitions triggers self
oscillations [10, 11], where an Andronov-Hopf bifurcation takes over the usual
phase transition.

As infection and recovery rates are changed, epidemic models on networks
display out-of-equilibrium phase transitions between a phase where a disease
is prevented from spreading and a phase where a macroscopic finite fraction of
the population becomes infected [12].

In this article, we will study the SIS and SIR models in a full microscopic
settings on random networks in presence of a feedback that changes the
structure of the underlying social network, and we will show that such
feedback triggers self-oscillations along the theory proposed in [10], where
suitably defined connectivity properties play the role of the control
parameter. In order to mimic the occurrence of lock downs, we will focus on a
a simple discontinuous feedback control, where a certain fraction of links is
deleted if the fraction of infections exceeds a given threshold $I_{2}$. The
same links are then be reinstated once the fraction of infections has been
reduced below a second threshold value $I_{1}<I_{2}$.

Oscillations in epidemic spreading have been studied mainly from the point of
view of seasonal effects that act as an external driving forces, while in the
case studied here oscillations are autonomously driven by an internal
feedback. The resulting models are described by time-independent equations and
parameters, and such oscillations can be considered as emerging self-
oscillations [13, 14, 15].

The article is organized as follows: In the first section we define the model
and illustrate its behavior with results from numerical simulations on an
instance of a real social network. In the second section we then study the
model on mean-field uncorrelated networks, where we will show that the overall
effect of lock downs on the dynamics is captured by a renormalization of the
effective infection rate through a change of the network branching ratio.

This finding is then exploited in the third section, where we analyze simple
piece-wise well-mixed models with point transformation techniques, leading to
analytical formulae for the number of waves and their length. In the last
section, we consider the realistic case of imperfect information on the
infection state, and we infer parameters from data on the current evolution of
the covid-19 pandemic, for which we estimate the extent of social network
disruption and characterize its limit cycle attractors.

## Results

### Model definition

We consider compartmentalized epidemic models on random networks, specifically
the SIS and SIR models (see [16] for a review), where individual agents are
represented as the nodes of a social network and can be in different states:
Infected, Susceptible and Recovered. Infected individuals recover with a
Poissonian rate $\gamma$, which is a parameter of the model, becoming either
susceptible (in the SIS model) or recovered (in the SIR model), and they
infect neighboring susceptible nodes with a Poissonian rate $\beta$, which is
the second parameter of the model. In order to model containment measures and
their relaxation, we consider a feedback between the network structure and the
infected state and its history in the following terms:

  * •

Starting from a state with few infections (whose relative number we will
indicate with $I$) that fastly spread, if the spreading overcomes a certain
threshold $I>I_{2}$ a central authority decides to disrupt the network
structure by randomly removing a macroscopic fraction $q$ of the links. This
will eventually revert back the spreading.

  * •

Starting from a regressing infected state in a disrupted network, when the
infection state is reverted to an acceptably low value $I<I_{1}$, the network
structure is restored back to its initial conditions.

Figure 1: $A\to B$: In a dense social network few infected highly contagious
people give rise to an epidemic spreading. $B\to C$: during the epidemic
outbreak a centralized authority decides for containment measures by severing
the network. $C\to D$ Under confinement, the epidemic regresses to few cases.
$D\to A$ Once the epidemic is supposedly under control containment measures
are withdrawn and the social network is restored. Red dots represent infected
individuals, while white dots represent susceptible ones. Light grey links in
panels C and D represent links that have been removed because of the
containment measure.

For the SIS model, this will eventually lead to an infection cycle as
illustrated in Fig. 1, where we show results of simulations on a school
friendship network reconstructed in [17] (number of nodes $N=134$).

In Fig. 2, we show instead a simulation of the SIR model on the same network
with and without feedback control, thus illustrating the effect of enforcing
containment measures. The control successfully reduces the spreading of the
infection, but for this to occur a series of lock downs have to be put in
place.

Figure 2: Fraction of infected and subsequently recovered individuals as a
function of time (that can be measured in days) from epidemic simulations of
the SIR model in a school friendship network [17] (number of nodes $N=134$)
with (left) and without (right) lock down measures in place, with parameters
$\beta=0.05$ (bare infection rate), $\gamma=0.07$ (recovery rate). The lock
down for the feedback case is enforced if the number of infections is above
$13$ and it is relaxed if they are below $2$, and it consists in a dilution of
the network links by a factor $q=0.99$. There are three lock downs occurring
at times $t=25$, $t=80$ and $t=130$.

We will consider in the next section the SIS and SIR models for the case of
large uncorrelated random networks, where it is possible to characterize
analytically general features of the dynamics.

### Networks

We consider here the case of large annealed uncorrelated random graphs with
degree distribution $P(k)$. At odds with static networks, in annealed networks
we assume that links are randomly rewired over a faster time scale of the
spreading process, while the assumption of uncorrelated networks implies there
are no correlations between the degrees of neighboring nodes. The locally
treee-like structure of these networks makes it possible to make analytical
progress in the study of dynamical processes taking place on them, since it
allows to recur to well-controlled approximations for the factorization of
probability states. In particular, we will consider here the heterogeneous
mean-field approximation, where nodes are grouped in classes according to
their degree.

Let us start from the SIS model. In absence of feedback, the rate equation for
the fraction $I_{k}$ of infected individuals of degree $k$ can be written as
follows [18]

$\dot{I}_{k}(t)=\beta\left(1-I_{k}(t)\right)k\Theta(t)-\gamma I_{k}(t),$ (1)

where $\Theta(t)=\sum_{k}\frac{kP(k)}{\langle k\rangle}I_{k}(t)$ is the
probability that a randomly selected neighbor of a node of degree $k$ is
infected, and we denote by $\langle k\rangle$ the average degree of the
network.

Here we consider the case in which, when the fraction of infected individuals
$I(t)=\sum_{k}P(k)I_{k}(t)$ exceeds a given threshold $I_{2}$, a lock down
measure is implemented that removes a fraction $q$ of links, which are then
reinstated once the condition $I(t)<I_{1}$ is satisfied. The equations of the
model in presence of this feedback mechanism can therefore be written in terms
of a state-dependent infection rate as follows

$\dot{I}_{k}(t)=\tilde{\beta}(I,\dot{I})\left(1-I_{k}(t)\right)k\Theta(t)-\gamma
I_{k}(t),$ (2)

where

$\tilde{\beta}(I,\dot{I})=\begin{cases}\beta,{\rm if}\leavevmode\nobreak\
I(t)<I_{1}\leavevmode\nobreak\ {\rm or}\leavevmode\nobreak\
\begin{cases}I_{1}\leq I(t)\leq I_{2}\\\ \dot{I}(t)>0\end{cases}\\\
(1-q)\beta,\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
if}\leavevmode\nobreak\ I(t)>I_{2}\leavevmode\nobreak\ {\rm
or}\begin{cases}I_{1}\leq I(t)\leq I_{2}\\\
\dot{I}(t)<0\end{cases}.\end{cases}$ (3)

We note that we can express $\tilde{\beta}$ as a function of the fraction of
infected population and its derivative because we are considering
deterministic rate equations. A more general representation for the discrete
stochastic case would require the introduction of a binary state variable to
denote the occurrence or absence of a lock down.

In Figure 3, we compare the result of numerical simulations with the numerical
solution of equation (2) for Erdős-Rényi and scale-free random networks. In
both cases we clearly see the emergence of oscillations due to the feedback.

Figure 3: Fraction of infected individuals as a function of time in an
homogeneous and heterogeneous network (both with $N=10^{5}$ nodes ). Dots
refer to 100 simulations of the SIS model (parameters $\beta=0.03$
$\gamma=0.08$ $I_{1}=0.05$ $I_{2}=0.2$ $q=0.95$), solid lines refer to the
numerical integration of the mean-field equations . Left: Erdős-Rényi random
network of average degree $\langle k\rangle=8$. Right: scale-free network of
minimum degree $k_{min}=2$ and exponent $\alpha=2.5$.

In Figure 4, we visualize the feedback-induced oscillations by means of a
phase portrait, where we plot the fraction of infected individuals vs. the
fraction of new positives. The figure clearly shows the emergence of a limit
cycle as an attractor of the dynamic. We also note from both Figures 3 and 4
that the dynamics on scale-free networks display bigger sample-to-sample
fluctuations than that on Erdős-Rényi networks.

Figure 4: Left panel: Phase portrait of an SIS model in an Erdős-Rényi network
with $N=10^{5}$ and average degree $\langle k\rangle=8$. Right panel: Phase
portrait of an SIS model in a scale-free network with $N=10^{5}$, $\gamma=2.5$
$k_{\rm min}=2$. In both panels, blue dots refer to 100 simulations of the SIS
model (parameters $\beta=0.03$ $\gamma=0.08$ $I_{1}=0.05$ $I_{2}=0.2$
$q=0.95$). Red dots refer to the numerical integration of the mean-field
equations.

Figure 5: Fraction of infected and recovered individuals as a function of time
for the SIR model with feedback. Left panel: Erdős-Rényi network with
$N=10^{5}$ nodes. Data refer to 100 simulations with parameters $\beta=0.03$
$\gamma=0.08$ $I_{1}=0.2$ $I_{2}=0.01$ $q=0.95$. Right panel: scale-free
network with $N=10^{5}$ nodes, $\gamma=2.5$ and $k_{m}in=2$. Data refer to 100
simulations with parameters $\beta=0.03$ $\gamma=0.08$ $I_{1}=0.2$
$I_{2}=0.05$ $q=0.95$. Black points: fraction of infected individuals from
numerical simulations. Blue points: fraction of recovered individuals from
numerical simulations. Solid red line: fraction of infected individuals from
numerical integration of the mean-field equations. Dashed green line: fraction
of recovered individuals from numerical integration of the mean-field
equations.

The same feedback mechanism can be considered for the SIR model as well. If we
now define $R_{k}$ as the fraction of recovered nodes with degree $k$, the
rate equations that describe the dynamic of the SIR model are

$\displaystyle\dot{I}_{k}(t)$ $\displaystyle=$
$\displaystyle\tilde{\beta}\left(1-I_{k}(t)-R_{k}(t)\right)k\Theta(t)-\gamma
I_{k}(t)$ (4) $\displaystyle\dot{R}_{k}(t)$ $\displaystyle=$
$\displaystyle\gamma I_{k}(t),$ (5)

where, as before, $\Theta(t)=\sum_{k}\frac{kP(k)}{\langle
k\rangle}I_{k}(t)$111We are considering here the case of annealed networks,
where links are randomly rewired at each time step. In the case of a static
network, the factor $k$ in the definition of $\Theta$ would be replaced by a
factor $k-1$. and $\tilde{\beta}(I,\dot{I})$ is given as before by equation
(3).

In figure 5 we show the evolution over time of the fraction of infected and
recovered individuals for the case of Erdős-Rényi and scale-free networks. We
see that the introduction of the feedback can lead to oscillations
corresponding to multiple infection waves. Clearly, in contrast to the case of
the SIS model, these oscillations will eventually come to an end once a large
enough fraction of the population has been infected. The number of infection
waves depends on the parameters of the model. In the next section we provide
an analytical estimation for a well-mixed population in the limit when $q$ is
close to $1$.

### Mixed approximation

The analysis reported in the previous section, where the effect of network
embedding and the feedback enacting on it is captured by a simple
renormalization of the bare infection rate, suggests to study the model on
well-mixed populations to get analytical insights. In this section $I$, $S$
and $R$ will indicate the fraction of infected, susceptible and recovered
individuals in the total population, respectively. The feedback law mimicking
containment, with parameters $\beta_{1}>\beta_{0}$ and $I_{2}>I_{1}$, is given
by

$\tilde{\beta}(I,\dot{I})=\begin{cases}\beta_{1},{\rm if}\leavevmode\nobreak\
I(t)<I_{1}\leavevmode\nobreak\ {\rm or}\leavevmode\nobreak\
\begin{cases}I_{1}\leq I(t)\leq I_{2}\\\ \dot{I}(t)>0\end{cases}\\\
\beta_{0},\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
if}\leavevmode\nobreak\ I(t)>I_{2}\leavevmode\nobreak\ {\rm
or}\begin{cases}I_{1}\leq I(t)\leq I_{2}\\\
\dot{I}(t)<0\end{cases}.\end{cases}$ (6)

For this form of the function $\tilde{\beta}(I,\dot{I})$, mixed epidemic
models can be readily analytically solved piece-wise in each sector, where the
infection rate is constant, and the solutions can be joined at the boundaries
(see Andronov [2]).

For instance, in the mixed approximation the SIS model with no feedback and
infection rate $\beta$ has the simple solution ($S+I=1$)

$I_{\beta}(t)=\frac{1-\gamma/\beta}{\left(\frac{1-\gamma/\beta}{I(0)}-1\right)e^{-(\beta-\gamma)t}+1}$
(7)

In presence of the feedback, the piece-wise constructed solution shows that
for $\beta_{1}>\gamma>\beta_{0}$, $I_{2}<1-\gamma/\beta_{1}$ the dynamics
settles into a limit cycle, and the periods of the quiescent epidemic
spreading ($t_{1}$) and of the recovery under lock downs ($t_{2}$) are given
by the following analytical formulae

$\displaystyle
t_{1}=\frac{1}{\beta_{1}-\gamma}\log\left(\frac{\frac{1-\gamma/\beta_{1}}{I_{1}}-1}{\frac{1-\gamma/\beta_{1}}{I_{2}}-1}\right)$
(8) $\displaystyle
t_{2}=\frac{1}{\beta_{0}-\gamma}\log\left(\frac{\frac{1-\gamma/\beta_{0}}{I_{2}}-1}{\frac{1-\gamma/\beta_{0}}{I_{1}}-1}\right)$
(9)

For a swift and resolute population lock down, we can approximate
$I_{2}<<1-\gamma/\beta_{1}$, $\beta_{0}<<\gamma$ and obtain for the total
duration of an epidemic wave

$T=t_{1}+t_{2}\sim\frac{\log(I_{2}/I_{1})}{\gamma}\frac{R_{0}}{R_{0}-1}.$ (10)

where $R_{0}=\beta_{1}/\gamma$. For instance, from the values $R_{0}\sim 3$, ,
$I_{2}/I_{1}\sim 100$ and $1/\gamma\sim 2$ weeks we can calculate $T\sim 3$
months.

For the SIR model, the solution in each interval – starting from initial
conditions $S_{i},R_{i},I_{i}$ at time $t_{i}$ – reads

$\displaystyle I+R+S$ $\displaystyle=$ $\displaystyle 1$ (11) $\displaystyle
S$ $\displaystyle=$ $\displaystyle S_{i}e^{\beta/\gamma(R_{i}-R)}$ (12)
$\displaystyle t-t_{i}$ $\displaystyle=$
$\displaystyle\int_{R_{i}}^{R}\frac{dr}{1-r-S_{i}e^{\beta/\gamma(R_{i}-r)}}.$
(13)

The total number of lock downs can be worked out analytically by joining
solutions piece-wisely, and a first order expansion in
$(\beta_{0}/\beta_{1},\beta_{0}/\gamma)$ (see the appendix) gives the formula

$n^{*}(\beta_{0})\sim n^{*}(\beta_{0}=0)/x,$ (14)

where

$n^{*}(0)=\frac{1-(1+\log R_{0})/R_{0}}{I_{2}-I_{1}}$ (15)

and

$1-\beta_{0}/\beta_{1}\leq x\leq 1+\beta_{0}/\beta_{1}+\beta_{0}/\gamma.$ (16)

For instance for the values $R_{0}\sim 3$, $I_{2}-I_{1}\sim 0.1\pm 0.05$ we
get $n^{*}\sim 3\pm 2$. In the next section we will illustrate data modeling
applications of our framework.

### Imperfect information and data modeling

In this section, we illustrate our framework in the context of modeling
epidemic data of the covid-19 infection in 2020 in the Italian region of
Lombardy and the Spanish region of the Basque Country. Data includes daily
reports of new infections and active cases 222From https://github.com/pcm-
dpc/COVID-19 (Lombardy) and
https://opendata.euskadi.eus/catalogo-datos (Basque country), plus a single
prevalence estimate. For the purpose of data analysis, we consider a variant
of the model with imperfect information by splitting the total number of
infections into detected and undetected cases, whose numbers we denote by
$I_{d}$ and $I_{u}$ respectively. We assume the existence of a detection
process by which undetected infected individuals are spotted with rate $r$ and
then put in isolation (which is equivalent to removing them). The mean-field
rate equations are as follows

$\displaystyle\dot{I}_{u}$ $\displaystyle=$
$\displaystyle\tilde{\beta}(I_{d},\dot{I_{d}})I_{u}S-(\gamma+r)I_{u}$ (17)
$\displaystyle\dot{I}_{d}$ $\displaystyle=$ $\displaystyle rI_{u}-\gamma
I_{d}$ (18)

where $\tilde{\beta}(I_{d},\dot{I_{d}})$ is the piece-wise constant function
defined in equation (6) (parametrized by $\beta_{1},\beta_{0},I_{1},I_{2}$),
and analogous equations for $S$ and/or $R$, depending on whether we consider
the SIS and/or the SIR model.

The inference of the model parameters has been performed by approximately
solving Bayes equations under the hypothesis of Gaussian noise by means of
Monte Carlo methods and dynamical system numerical simulations (see the
appendix).

Figure 6: New daily cases vs time for the first epidemic wave of COVID-19:
data from Lombardy (left) and Basque country (right) against the inferred
model (maximum likelihood estimate).

Figure 7: Scatter plots of new daily and active cases for Lombardy (left) and
Basque country (right) for the first epidemic wave of covid-19. In both
panels, the solid blue line refers to the model inferred maximizing the
likelihood, red dots refer to daily data and green stars refer to weekly data.

In Fig. 6 we show the time series of new daily cases for the two analyzed
cases against the model with maximum likelihood parameters, the latter showing
a clear cusp peak and piece-wise exponential trends corresponding to the lock
down event. The model suggests that a closed trajectory must be observed in
the plane of new and active daily cases. This is shown in Fig. 7, where we
report data smoothened by a 7-days moving average as well. This closed
trajectory in the ideal case of a pure SIS model with feedback would be a
limit cycle attractor of the dynamics.

## Conclusions

Connectivity plays a crucial role in the definition of the parameters that
control the collective behavior of a system. This finding has striking
consequences, like the absence of an epidemic threshold in epidemic spreading
models defined on scale-free networks [19], and it suggests that it is
possible to control the spreading by acting on the network of social
interactions.

In this article we have shown that feedback control at the level of the social
network in epidemic models triggers self-oscillations along the theory
proposed in [10].

We have investigated self-oscillations induced by a simple discontinuous
feedback control mimicking lock down events in classical compartmentalized
epidemic models (SIS and SIR) on networks.

On random graphs, for Erdős-Rényi as well as scale-free networks with naive
populations, we have shown that the effect of lock downs simply amounts at
renormalizing the effective infection rate to account for the reduction in the
network branching ratio.

This led to simple piece-wise mean field approximations that we solved
analytically by means of transformation point methods, recovering formulae for
the number of waves and their extent in terms of the model parameters. These
formulae can be in principle tested against data, once a certain amount of
evidence accumulates on the number of lock downs and their lengths, prevalence
estimates, and basic infection numbers region by region.

A problem related to data collection during the covid-19 epidemic outbreak was
the fact that many positives were undetected. In order to bring the model to
data, we have therefore extended it by assuming the existence of a fraction of
undetected positives, who can then be detected at a given rate (for instance
through testing). We have applied our extended framework to analyze data from
the first epidemic wave of covid-19 in Lombardy and Basque country, where
parameters have been inferred leading to a characterization of the dynamical
attractors in the phase space. Apart from applications to predictive modeling
– which would require more extensive data analysis [20] and methods of system
identification- - we do point out here briefly some potentially interesting
theoretical problems stemming from this work.

First, the issue of optimal scheduling [21] in the control of the social
network, leading to continuous feedback and potentially smoother oscillations.
Current qualitative evidence from the second epidemic wave of covid-19 seems
indeed to show in some regions smoother trends and oscillations around the
phase transition point ($R_{t}\sim 1$), due to attempts of finer control like
partial restrictions and selected closures taken in due course. Recently
proposed analytical frameworks [7] can be very useful in this respect.

Second, within the framework proposed here, where epidemics can be regarded as
self-oscillators, it comes naturally the question of coupling and
synchronization [22, 23] of epidemic waves running on different networks that
are weakly connected, e.g. by migration processes.

Finally, another interesting issue concerns the impact of periodically
external drive on oscillators: analogously to well-known forced double well
oscillator [24], the combined effect of seasonal changes and feedback could
potentially lead to chaotic oscillations in strange attractors, an aspect that
adds to the problem of predictability of such systems, and that we leave for
future investigations.

##  Appendix A Number of epidemic waves: perturbative expansion

In this section we derive formula (14-16). In the SIR model at fixed infection
rate $\beta$ (as in a simple model without feedback, or in a given interval
for the piece-wise feedback model) the peak value of the infected fraction is
(when $\dot{I_{p}}=0$ and $S_{p}=\gamma/\beta$)

$I_{p}=I_{i}+S_{i}-\frac{\gamma}{\beta}\log
S_{i}-\gamma/\beta(1-\log(\gamma/\beta)).$ (19)

This value is not achieved if it is greater than the one triggering the lock
down, i.e. when $I_{p}>I2$. Thus we will assume as halting condition that
$I_{p}\leq I_{2}$, since in this case no lock down takes place and the system
proceeds towards herd immunity. We will now work out a series for fraction of
susceptible individuals at the various stages of the epidemic waves,
exploiting the piece-wise analytical solutions in each interval.

Suppose we are at the beginning of a wave $I=I_{1}$ with given susceptible
fraction $S=S_{n,-}$, the system (with $\beta=\beta_{1}$) will evolve towards
$I_{2}$ and a given $S_{n,+}$ that satisfies

$S_{n,+}-\frac{\gamma}{\beta_{1}}\log
S_{n,+}=S_{n,-}-\frac{\gamma}{\beta_{1}}\log S_{n,-}-\Delta I,$ (20)

where $\Delta I=I_{2}-I_{1}$. Then we have the lock down $\beta=\beta_{0}$,
and the system will evolve towards $I_{1}$ with a given $S_{n+1,-}$ that
satisfies

$S_{n+1,-}-\frac{\gamma}{\beta_{0}}\log
S_{n+1,-}=S_{n,+}-\frac{\gamma}{\beta_{0}}\log S_{n,+}+\Delta I$ (21)

These equations define a series eventually halting when $I_{p}\leq I_{2}$.

### Vanishing $\beta_{0}$

Suppose $\beta_{0}=0$. In this case $S_{n+1,-}=S_{n,+}\equiv S_{n}$ and we
have

$\displaystyle S_{n+1}-\frac{\gamma}{\beta_{1}}\log S_{n+1}$ $\displaystyle=$
$\displaystyle S_{n}-\frac{\gamma}{\beta_{1}}\log S_{n}-\Delta I$ (22)
$\displaystyle=$ $\displaystyle S_{n-1}-\frac{\gamma}{\beta_{1}}\log
S_{n-1}-2\Delta I$ (23) $\displaystyle=$ $\displaystyle\dots$ (24)
$\displaystyle=$ $\displaystyle S_{0}-\frac{\gamma}{\beta_{1}}\log
S_{0}-n\Delta I$ (25)

If we start from $S_{0}\approx 1$, from the halting condition we find that the
number $n^{*}$ of lockdowns is

$n^{*}=\left\lfloor\frac{1-(1+\log R_{0})/R_{0}}{\Delta I}\right\rfloor,$ (26)

where we denote by $\left\lfloor x\right\rfloor$ the integer part of $x$.

### First order expansion in $\beta_{0}$

A first order expansion in $\beta_{0}$ leads to

$S_{n+1,-}\sim S_{n,+}(1-\beta_{0}/\gamma\Delta I).$ (27)

Then defining

$F_{n}=S_{n,-}-\gamma/\beta_{1}\log S_{n,-}$ (28)

we have the recursion relation

$F_{n+1}=F_{n}-x_{n}\Delta I$ (29)

where

$x_{n}=1-\beta_{0}/\beta_{1}+\beta_{0}/\gamma S_{n+1,-}$ (30)

that can be bounded by $\mathcal{O}(\beta_{0})$ terms (given that $0\leq
S_{n,-}\leq 1$)

$1-\beta_{0}/\beta_{1}\leq x_{n}\leq 1-\beta_{0}/\beta_{1}+\beta_{0}/\gamma$
(31)

and the halting criterion leads to

$\displaystyle n^{*}(\beta_{0})\sim n^{*}(\beta_{0}=0)/x$ (32) $\displaystyle
1-\beta_{0}/\beta_{1}\leq x\leq 1-\beta_{0}/\beta_{1}+\beta_{0}/\gamma$ (33)

In Figure 8 we show the agreement between numerical simulations on an Erdős-
Rényi network and the analytical prediction given by equation (26) for
$\beta_{0}=0$.

Figure 8: Number of lock downs: comparison between analytical and numerical
results obtained for an Erdős-Rényi random network with $N=10^{5}$ nodes,
average degree $8$, $\beta_{0}=0$, $I_{1}=0.2$, $I_{2}=0.05$. The four points
correspond to $\beta_{1}\in\\{0.02,0.03,0.06,0.09\\}$.

##  Appendix B Inference of model parameters from epidemic data

Here we consider the task of fitting against epidemic data a model including
imperfect information on the state of the system. This is an instance of a
system identification problem [25], which we solved along the following lines:
We consider the time series of observed new daily and active cases
$(N_{t}^{o},A_{t}^{o})$ ($t=1\dots T$ is the temporal index in days, starting
from the 1st of March, $T=T_{l}=155$ for Lombardy, and $T=T_{b}=100$ for the
Basque country), and we assume it as coming from an instance of the model plus
a noise term

$\displaystyle N_{t}^{o}=rI_{u}(t,\wp)+\delta_{N}$ (34) $\displaystyle
A_{t}^{o}=I_{d}(t,\wp)+\delta_{A}$ (35)
$\displaystyle\wp=\left\\{\gamma,\beta_{1},\beta_{0},rI_{1},I_{2},I_{u}(0),I_{d}(0)\right\\},$
(36)

where we highlighted the dependence of the model trajectory by the dynamical
parameters and boundary values.

We assume shot-noise of the form

$\displaystyle\langle\delta_{N}\rangle=\langle\delta_{A}\rangle=0$ (37)
$\displaystyle\langle\delta_{N}^{2}\rangle=N\quad\langle\delta_{A}^{2}\rangle=A,$
(38)

which for large numbers we assume to be distributed normally. Upon assuming an
uniform prior, we have the following formula for the log-likelihood of the
parameters

$\mathcal{L}(\wp)=\sum_{t}\frac{(N_{t}^{o}-rI_{u}(t,\wp))^{2}}{2N_{t}^{o}}+\frac{(A_{t}^{o}-I_{d}(t,\wp))^{2}}{2A_{t}^{o}}+\textrm{const.}$
(39)

From the Bayes formula, the posterior probability distribution of parameters
$P(\wp)\propto e^{-\mathcal{L}(\wp)}$ has been sampled by a Metropolis
Montecarlo rule, where the evaluation of $\mathcal{L}(\wp)$ has been done by
numerical integration of the model equations. More explicitly we have been
following the following flowchart:

  * •

Start from some value of the parameters $\wp_{0}$: a warm start has been
provided by fitting the curve of new daily cases alone in linear approximation
($S<<N$).

  * •

Propose a change for the parameters $\wp_{n}\to\wp_{n+1}$: we used independent
geometrical random walks of stepsize $10^{-3}$.

  * •

Numerically integrate the model equations with the new proposed parameters to
evaluate their log-likelihood $\mathcal{L}(\wp_{n+1})$. We used the standard
Verlet algorithm.

  * •

Accept the proposed new parameters with probability
$\textrm{min}(1,\exp{\mathcal{L}(\wp_{n+1})-\mathcal{L}(\wp_{n})})$
(Metropolis rule), otherwise keep the old parameters.

This defines a series that asymptotically uniformly samples the posterior
probability for the parameters, whose peak values have been used for the
results showed in Fig. 6 and 7. Finally, we do point out that when we
numerically integrate the model equations we rejected solutions that i) do not
include at least one lock down event and ii) do not agree with prevalence
estimate from serological data. The latter gives a lower bound for the total
number $K$ of infected individuals at a certain time, $K\geq K_{o}$ that can
be reformulated as an inequality between the model parameters (in particular
$r$, $\beta_{1}$ and $\beta_{0}$) as follows: $K$ can decomposed in detected
and undetected cases $K=K_{d}+K_{u}$, where $K_{d}=\sum\Delta I_{d}^{+}$ is
given by the data, while for the latter we have

$\displaystyle K_{u}=\sum_{t}\Delta I_{u,t}^{+}=\sum\beta(I_{d,t})I_{u,t}=$
(40) $\displaystyle=\frac{1}{r}\sum\beta(I_{d,t})\Delta I_{d}^{+}=$ (41)
$\displaystyle=\frac{\beta_{1}}{r}\sum_{incr}\Delta
I_{d}^{+}+\frac{\beta_{0}}{r}\sum_{decr}\Delta I_{d}^{+}=$ (42)
$\displaystyle=\frac{K_{d}}{r}(\beta_{1}(1-x)+\beta_{0}x)$ (43)

In the first and second lines we have used the model hypothesis, in the third
we have decomposed the sum in terms of the increasing and decreasing part of
the wave, and finally in the last line $x$ is the fraction of the detected
infections during the decreasing part of the epidemic wave. We have finally
the inequality

$r(K_{o}/K_{d}-1)\leq\beta_{1}(1-x)+\beta_{0}x$ (44)

Given the simplicity of the model employed – in particular with respect to the
hypothesis of a constant detection rate $r$ – we obtain a fairly high value of
the $\chi^{2}\sim 20$, with a concomitant acceptably low average relative
error of $\epsilon\sim 20\%$ across data points that can be considered apt for
a qualitative description of the data. We report in the following table the
maximum likelihood inferred parameters with their standard deviation.

Region | $R_{0}=\beta_{1}/\gamma$ | $1-p=\beta_{0}/\beta_{1}$ |  $r$(day${}^{-1})$
---|---|---|---
Basque country | $2.9\pm 0.1$ | $1.3\pm 0.4\cdot 10^{-2}$ | $3\pm 1\cdot 10^{-3}$
Lombardy | $5.2\pm 0.2$ | $4\pm 1\cdot 10^{-2}$ | $8.4\pm 0.6\cdot 10^{-3}$

We do point out an anomalously high $R_{0}$ for Lombardy, due to a very low
average recovery rate of approximately $1$ month (while the one from Basque
country is around two weeks in agreement with WHO estimates). This is apparent
upon looking at the much slower decay of the infection curves during the lock
down, and it is probably due to a biased over-sampling of critically ill cases
whose average recovery is typically longer.

## References

  * [1] James Clerk Maxwell.  I. on governors.  Proceedings of the Royal Society of London, (16):270–283, 1868\.
  * [2] Aleksandr Aleksandrovich Andronov, Aleksandr Adol’fovich Vitt, and Semen Emmanuilovich Khaikin.  Theory of Oscillators.  Dover, 1966.
  * [3] Karl Johan Åström and Richard M Murray.  Feedback systems: an introduction for scientists and engineers.  Princeton university press, 2010.
  * [4] Simone Pigolotti, Sandeep Krishna, and Mogens H Jensen.  Oscillation patterns in negative feedback loops.  Proceedings of the National Academy of Sciences, 104(16):6533–6537, 2007.
  * [5] Keith N Frayn.  Metabolic regulation: a human perspective.  John Wiley & Sons, 2009.
  * [6] James D Murray.  Mathematical biology: I. An introduction, volume 17.  Springer Science & Business Media, 2007.
  * [7] Ginestra Bianconi and Pavel L Krapivsky.  Epidemics with containment measures.  arXiv preprint arXiv:2004.03934, 2020.
  * [8] Francesco Di Lauro, István Z. Kiss, Daniela Rus, and Cosimo Della Santina.  Covid-19 and flattening the curve: a feedback control perspective.  arXiv preprint arXiv:2020.3039322, 2020.
  * [9] Yunxiu Zhou, Simon A Levin, and Naomi E Leonard.  Active control and sustained oscillations in actsis epidemic dynamics.  arXiv preprint arXiv:2007.01424, 2020.
  * [10] Daniele De Martino.  Feedback-induced self oscillations in large interacting systems subjected to phase transitions.  Journal of Physics A: Mathematical and Theoretical, 2018.
  * [11] Daniele De Martino and Andre C Barato.  Oscillations in feedback-driven systems: Thermodynamics and noise.  Physical Review E, 100(6):062123, 2019.
  * [12] Haye Hinrichsen.  Non-equilibrium critical phenomena and phase transitions into absorbing states.  Advances in physics, 49(7):815–958, 2000.
  * [13] E De Lauro, S De Martino, M Falanga, and M Palo.  Self-sustained vibrations in volcanic areas extracted by independent component analysis: a review and new results.  Nonlinear Processes in Geophysics, 18(6):925–940, 2011.
  * [14] G Buccheri, E De Lauro, S De Martino, and M Falanga.  Experimental study of self-oscillations of the trachea–larynx tract by laser doppler vibrometry.  Biomedical Physics & Engineering Express, 2(5):055009, 2016.
  * [15] Alejandro Jenkins.  Self-oscillation.  Physics Reports, 525(2):167–222, 2013.
  * [16] Mark Newman.  Networks.  Oxford university press, 2018.
  * [17] Rossana Mastrandrea, Julie Fournet, and Alain Barrat.  Contact patterns in a high school: a comparison between data collected using wearable sensors, contact diaries and friendship surveys.  PloS one, 10(9):e0136497, 2015.
  * [18] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani.  Epidemic processes in complex networks.  Reviews of modern physics, 87(3):925, 2015.
  * [19] Romualdo Pastor-Satorras and Alessandro Vespignani.  Epidemic spreading in scale-free networks.  Physical review letters, 86(14):3200, 2001.
  * [20] Jose M. G. Vilar and Leonor Saiz.  The evolving worldwide dynamic state of the covid-19 outbreak.  medrXiv preprint 2020.11.26.20239434, 2020.
  * [21] Donald E Kirk.  Optimal control theory: an introduction.  Courier Corporation, 2004.
  * [22] Francesco Guerra.  Coupled self-oscillating systems: Theory and applications.  International Journal of Modern Physics B, 23(28n29):5505–5514, 2009.
  * [23] Francesca Di Patti, Duccio Fanelli, Filippo Miele, and Timoteo Carletti.  Ginzburg-landau approximation for self-sustained oscillators weakly coupled on complex directed graphs.  Communications in Nonlinear Science and Numerical Simulation, 56:447–456, 2018.
  * [24] John Guckenheimer and Philip Holmes.  Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42.  Springer Science & Business Media, 2013.
  * [25] Karel J Keesman.  System identification: an introduction.  Springer Science & Business Media, 2011.