# Inferring High-Dimensional Dynamic Networks Changing with Multiple
Covariates

Louis Dijkstra  Leibniz Institute for Prevention Research & Epidemiology –
BIPS, Bremen, Germany  Arne Godt  Leibniz Institute for Prevention Research &
Epidemiology – BIPS, Bremen, Germany  Ronja Foraita111Corresponding author.
E-mail<EMAIL_ADDRESS>Leibniz Institute for Prevention Research &
Epidemiology – BIPS, Bremen, Germany

###### Abstract

High-dimensional networks play a key role in understanding complex
relationships. These relationships are often dynamic in nature and can change
with multiple external factors (e.g., time and groups). Methods for estimating
graphical models are often restricted to static graphs or graphs that can
change with a single covariate (e.g., time). We propose a novel class of
graphical models, the covariate-varying network (CVN), that can change with
multiple external covariates.

In order to introduce sparsity, we apply a $L_{1}$-penalty to the precision
matrices of $m\geq 2$ graphs we want to estimate. These graphs often show a
level of similarity. In order to model this ‘smoothness’, we introduce the
concept of a ‘meta-graph’ where each node in the meta-graph corresponds to an
individual graph in the CVN. The (weighted) adjacency matrix of the meta-graph
represents the strength with which similarity is enforced between the $m$
graphs.

The resulting optimization problem is solved by employing an alternating
direction method of multipliers. We test our method using a simulation study
and we show its applicability by applying it to a real-world data set, the
gene expression networks from the study ‘German Cancer in childhood and
molecular-epidemiology’ (KiKme). An implementation of the algorithm in R is
publicly available under github.com/bips-hb/cvn.
Keywords: alternating direction method of multipliers $\cdot$ covariate-
varying networks $\cdot$ Gaussian graphical model $\cdot$ graphical LASSO

##  1 Introduction

In many scientific fields, network models are essential for describing and
understanding complex relationships (Newman,, 2018). In genetic and molecular
epidemiology, networks elucidate interactions between molecular regulators
governing gene expression such as gene regulatory networks (Barabasi and
Oltvai,, 2004) and how microbial communities in the microbiome interact with
each other and the host (e.g. Curtis et al.,, 2012; de Vos and de Vos,, 2012).
In neuroscience, network models map intricate brain connections, aiding in
understanding communication and coordination between different regions
(Badhwar et al.,, 2017).

Extensive research has been devoted to inferring high-dimensional networks
from data (Friedman et al.,, 2008; Cai et al.,, 2011; Danaher et al.,, 2014),
primarily on static networks where edges remain fixed and unaffected by
external influences. In this paper, we use the terms ‘graph’ and ‘network’
interchangeably. While these static networks provide valuable insights, it is
crucial to also understanding how network structures evolve with respect to
external covariates, such as time or group differences. For example, long-term
exposure to ultraviolet (UV) radiation can alter gene networks, impairing DNA
repair mechanisms and potentially leading to skin cancer (Brozyna et al.,,
2007). Microbiome networks differ between healthy individuals and those with
obesity (de la Cuesta-Zuluaga et al.,, 2018). In neuroepidemiology, changes in
brain connectome interactions due to cholesterol levels and aging indicate
Alzheimer’s disease (Badhwar et al.,, 2017). Observing how dynamic networks
change with external factors – such as UV radiation, treatment regimes, or
cholesterol levels – is fundamental to understanding the underlying biological
processes.

Current literature provides, to be best of our knowledge, only methods for
analyzing networks that change with a single discrete external covariate, such
as time (Monti et al.,, 2014; Hallac et al.,, 2017; Gibberd and Nelson,, 2017)
or groups (Danaher et al.,, 2014; Hallac et al.,, 2017; Vinciotti et al.,,
2024). Here, we propose extending these methods to multiple, possibly
intricately related, discrete covariates. We present a very general framework
for defining a dynamical graphical model of this kind. We call this the
_covariate-varying network_ (CVN) model. We present an alternating direction
method of multipliers (ADMM) algorithm that can be used to estimate the model
on the basis of (potentially high-dimensional) data.

The paper is organized as follows: Section 2 introduces static Gaussian
graphical models in which the (in)dependency structure is unaffected by
external covariates; we simultaneously present much of the notation used
throughout the paper. In Section 3, we formally define the CVN graphical
model, including the concept of a _meta-graph_ to represent the similarities
between the various individual graphs. This section also discusses the
estimator considered throughout the paper. Additionally, we highlight several
known models from the literature that are special cases of the CVN framework.

We solve the estimation problem for the CVN using an ADMM algorithm that we
describe in detail in Section 4. The computational complexity of the algorithm
is determined in Section 5. We show how one could interpolate a graph given an
estimated CVN in Section 6. The two tuning parameters of the model, one
regulating sparsity and the other controling the level of similarity between
the graphs, need to be chosen. In Section 7, we show how this could be done
using either the Akaike or Bayesian information criterion (AIC and BIC). We
propose an alternative tuning parameterization in Section 8 that leads to more
stable choices that are barely influenced by the introduction of more
variables and/or covariates.

We assess the performance of the method in a simulation study in which we
explore to what extent the proposed method can reconstruct a CVN. We present
the simulator in Section 9 and the simulation set-up itself in Section 10. The
measures used to assess the performance are presented in the latter section as
well.

In addition, we apply the method to a real data set from the KiKme study
(Marron et al.,, 2021), see Section 11. The results of both the simulation and
this case study are shown in Section 12. We end with our conclusions and a
discussion in Section 13.

The ADMM algorithm has been implemented as an R package and is publicly
accessible at www.github.com/bips-hb/CVN. Additionally, the CVN data simulator
is available as an R package at www.github.com/bips-hb/CVNSim, and the code
for the simulation study can be found at www.github.com/bips-hb/CVNStudy. All
results of the simulation study are hosted at cvn.bips.eu.

##  2 Static Networks

We represent a static network as a probabilistic undirected graphical model
with the tuple $\\{\mathbf{X},f(\mathbf{X}),G=(V,E)\\}$ where
$\mathbf{X}=(X_{1},X_{2},\ldots,X_{p})^{\top}$ is a real-valued
$p$-dimensional random vector with joint density function $f(\mathbf{X})$. The
nodes of the graph $G$, $V=\\{1,2,\ldots,p\\}$, correspond to the variables
$X_{1},X_{2},\ldots,X_{p}$. The edges $E\subseteq V\times V$ reflect the
(in)dependence structure between the variables in accordance with the density
function $f(\mathbf{X})$. The presence or absence of an edge represents
whether the two corresponding variables are conditionally independent, i.e.,
the edge $\\{s,t\\}$ is in $E$ if $X_{s}\not\\!\perp\\!\\!\\!\perp
X_{t}\mid\mathbf{X}_{V\setminus\\{s,t\\}}$, where
$\mathbf{X}_{V\setminus\\{s,t\\}}$ are all variables in $\mathbf{X}$ except
$X_{s}$ and $X_{t}$. Let $\mathbf{A}=(a_{st})_{p\times p}$ be the graph’s
adjacency matrix where $a_{st}=\mathds{1}(\\{s,t\\}\in E)$ and
$\mathds{1}(\cdot)$ is the indicator function.

We assume the random vector $\mathbf{X}$ to follow a multivariate normal
distribution with mean $\bm{\mu}$ and covariance matrix $\mathbf{\Sigma}$.
Under this assumption, the entries of the precision matrix
$\mathbf{\Theta}=\left(\theta_{st}\right)_{p\times p}=\mathbf{\Sigma}^{-1}$
are proportional to the partial correlations between the variables.
Determining the network structure, i.e., the edge set $E$, is, therefore,
equivalent to determining which entries in the precision matrix
$\mathbf{\Theta}$ are zero and non-zero. In terms of the graph’s adjacency
matrix, $a_{st}=\mathds{1}(\theta_{st}\neq 0)$ for all $\\{s,t\\}\in V\times
V$ (see Lauritzen,, 1996). Without loss of generality, we assume that the data
are centered around the mean $\bm{\mu}=\mathbf{0}$ for the remainder of the
paper.

Suppose we observe $n$ realizations $\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{n}$
of the random vector $\mathbf{X}$. The maximum likelihood estimator of the
precision matrix is given by

$\widehat{\mathbf{\Theta}}=\operatorname*{arg\,min}_{\mathbf{\Theta}\succ
0}\ell(\mathbf{\Theta})=\operatorname*{arg\,min}_{\mathbf{\Theta}\succ
0}\leavevmode\nobreak\
\frac{n}{2}\left[\text{trace}\left(\widehat{\mathbf{\Sigma}}\mathbf{\Theta}\right)-\log\text{det}\left(\mathbf{\Sigma}\right)\right],$
(1)

where $\mathbf{\Theta}\succ 0$ denotes positive definiteness,
$\ell(\mathbf{\Theta})$ is the log-likelihood function, ‘trace’ and ‘det’
denote the trace and determinant of a matrix, and
$\widehat{\mathbf{\Sigma}}=n^{-1}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{\top}$ is
the empirical covariance matrix (see, for example Uhler,, 2018). However, this
estimator does lead to dense graphs. Moreover, in high-dimensional settings,
the empirical covariance matrix $\widehat{\mathbf{\Sigma}}$ is rank deficient
and cannot be used to estimate the precision matrix
$\widehat{\mathbf{\Theta}}$. One can estimate $\widehat{\mathbf{\Theta}}$ by
imposing a sparse graph structure and, hence, a sparse precision matrix. A
popular method to estimate sparse networks is the Graphical Least Absolute
Selection and Shrinkage Operator (GLASSO; Banjeree and Ghaoui, (2008) and
Friedman et al., (2008)) which places a $L_{1}$-norm penalty on the precision
matrix in order to ‘shrink’ its off-diagonal entries to zero. It involves
solving the following convex optimization problem:

$\widehat{\mathbf{\Theta}}_{\text{GLASSO}}=\operatorname*{arg\,min}_{\mathbf{\Theta}\succ
0}\ell(\mathbf{\Theta})+\lambda\left\lVert\mathbf{\Theta}\right\rVert_{1}^{\text{off}}.$
(2)

Here, $\lambda>0$ is the tuning parameter and
$\left\lVert\mathbf{\Theta}\right\rVert_{1}^{\text{off}}=\sum_{s\neq
t}|\theta_{st}|$ is the off-diagonal $L_{1}$-norm of the precision matrix.

##  3 The Covariate-Varying Network Model

We represent a covariate-varying network (CVN) graphical model with the
quintuple

$\text{CVN}=\\{\mathbf{X},\mathbf{U},\mathcal{U},f(\mathbf{X}\mid\mathbf{U}),\\{G(u)=(V,E(u))\\}_{u\in\mathcal{U}}\\},$
(3)

where $\mathbf{X}=(X_{1},X_{2}\ldots,X_{p})^{\top}$ is a $p$-dimensional
random vector and $\mathbf{U}=(U_{1},U_{2},\ldots,U_{q})^{\top}$ is a random
vector representing $q$ external discrete covariates, i.e., variables that are
not included in the graph with $(K_{1},\ldots,K_{q})^{\top}$ categories. The
vector $\mathbf{U}$ lies in the discrete space $\mathcal{U}$ with cardinality
$m\leq\prod_{k=1}^{q}K_{k}$. The joint density function of $\mathbf{X}$
conditioned on $\mathbf{U}$ is $f(\mathbf{X}\mid\mathbf{U})$. The fifth
element of the CVN is a set of $m$ graphs, one for each value $u$ in
$\mathcal{U}$. The vertices of $G(u)$, $V=\\{1,2,\ldots,p\\}$, correspond to
the variables $X_{1},X_{2},\ldots,X_{p}$ and do not change with $\mathbf{U}$.
The (in)dependence structure between the variables is captured by the edge set
$E(u)$ which can change with $\mathbf{U}$. Specifically, the edge $\\{s,t\\}$
is in $E(u)$ if $X_{s}\not\\!\perp\\!\\!\\!\perp X_{t}\mid U=u\ \text{and
}\mathbf{X}_{V\setminus\\{s,t\\}}$. We denote the adjacency matrix of graph
$G(u)$ as $\mathbf{A}(u)=(a_{st}(u))_{p\times p}$ where
$a_{st}(u)=\mathds{1}(\\{s,t\\}\in E(u))$. In this paper, we are interested in
estimating the graphs $\\{G(u)=(V,E(u))\\}_{u\in\mathcal{U}}$.

We assume throughout that $\mathbf{X}\mid\mathbf{U}=u$ follows a multivariate
normal distribution with mean $\bm{\mu}(u)=\mathbf{0}$ and covariance matrix
$\mathbf{\Sigma}(u)$. Under the normality assumption, the entries of the
precision matrices correspond with the edge sets of the different graphs. Let
$\mathbf{\Theta}(u)=\mathbf{\Sigma}(u)^{-1}=(\theta_{st}(u))_{p\times p}$ be
the precision matrix, then the edge $\\{s,t\\}$ is in $E(u)$ if
$\theta_{st}(u)\neq 0$. Hence, estimating the structure of the graphs
$\\{G(u)\\}_{u\in\mathcal{U}}$, corresponds to determining the zero and non-
zero entries in the precision matrices
$\\{\mathbf{\Theta}(u)\\}_{u\in\mathcal{U}}$.

We index each element in $\mathcal{U}$ by assigning them a unique element from
the index set $\mathcal{I}=\\{1,2,\ldots,m\\}$. If $i\in\mathcal{I}$ is the
index of $u\in\mathcal{U}$, we write $G_{i}=G(u)$, $E_{i}=E(u)$,
$\mathbf{\Sigma}_{i}=\mathbf{\Sigma}(u)$,
$\mathbf{\Theta}_{i}=\mathbf{\Theta}(u)$ and $\mathbf{A}_{i}=\mathbf{A}(u)$.
The $(s,t)$ entry in the $i$-th precision matrix is written as
$\theta_{st}^{(i)}$. The collection of all precision matrices is denoted by
$\mathbf{\Theta}_{\text{CVN}}=\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}$.

For each value of $u\in\mathcal{U}$, we obtain $n_{i}>0$ observations:
$\left(\bm{x}_{1}^{(i)},\bm{x}_{2}^{(i)},\ldots,\bm{x}_{n_{i}}^{(i)}\right)$
where $\bm{x}_{j}^{(i)}\in\mathds{R}^{p}$ and $j=1,\ldots,n_{j}$. Similarly to
eq. (1), we can express the log-likelihood function for precision matrix
$\mathbf{\Theta}_{i}$ as

$\ell_{i}(\mathbf{\Theta}_{i})=\frac{n_{i}}{2}\cdot\left[\text{trace}\left(\widehat{\mathbf{\Sigma}}_{i}\mathbf{\Theta}_{i}\right)-\log\text{det}\left(\mathbf{\Sigma}_{i}\right)\right],$
(4)

where
$\widehat{\mathbf{\Sigma}}_{i}=n{{}_{i}}^{-1}\sum_{j=1}^{n_{i}}\bm{x}_{j}^{(i)}(\bm{x}_{j}^{(i)})^{\top}$
is the empirical covariance matrix. Overall, the log-likelihood for the CVN
model itself is simply the sum over all graphs, i.e.,
$\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})$.

To introduce sparsity, we apply, just as for the GLASSO in eq. (2), the off-
diagonal $L_{1}$-norm penalty for each individual graph:

$\operatorname*{arg\,min}_{\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}}\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}},$

where $\lambda_{1}>0$ is the GLASSO tuning parameter. Larger values of
$\lambda_{1}$ increase the penalty and result in a sparser estimate of the
precision matrices.

In many applications, we expect some similarities between the graphs. For
example, two networks observed on two consecutive time points tend to have
similar edge sets, while networks observed on far removed time points show
less similarity. We refer to enforcing similarity between certain graphs as
_smoothing_. For introducing smoothness across the $m$ graphs, we first define
a _meta-graph_. Let $\mathcal{G}=(\mathcal{V},\mathcal{E},w)$ be a weighted,
undirected graph where each node in the node set
$\mathcal{V}=\\{1,2,\ldots,m\\}$ corresponds to the graphs
$G_{1},G_{2},\ldots,G_{m}$,
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$, and
$w:\mathcal{E}\rightarrow[0,1]$ is the weight function that assigns a value
between $0$ and $1$ to each edge $e\in\mathcal{E}$. We denote the symmetric
weighted adjacency matrix of the meta-graph by $\mathbf{W}=(w_{ij})_{m\times
m}$ where $w_{ij}=w_{ji}=w(\\{i,j\\})$.

The weight matrix reflects some prior knowledge about the similarity of how
the $m$ graphs might change with the external covariates. For a consistent
structure estimation, we must impose an additional smoothness constraint on
$\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}$ to detect graph-wise changes in the
underlying (in)dependence structure. We define the CVN estimator of
$\mathbf{\Theta}_{\text{CVN}}=\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}$ as

$\widehat{\mathbf{\Theta}}_{\text{CVN}}=\operatorname*{arg\,min}_{\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}}\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{i<j}w_{ij}\cdot\left\lVert\mathbf{\Theta}_{i}-\mathbf{\Theta}_{j}\right\rVert_{1}^{\text{off}},$
(5)

where $\lambda_{2}\geq 0$ is the _smoothing tuning parameter_. The term
$w_{ij}\left\lVert\mathbf{\Theta}_{i}-\mathbf{\Theta}_{j}\right\rVert_{1}^{\text{off}}$
penalizes the differences between the edge sets of graph $G_{i}$ and graph
$G_{j}$. The value $w_{ij}$ encourages a similar structure between graphs
$G_{i}$ and $G_{j}$ whereas $\lambda_{2}$ regulates how strongly differences
between all $\binom{m}{2}$ pairs of graphs are penalized.

###  3.1 Special Cases of Covariate-Varying Networks

By choosing the space $\mathcal{U}$, the indexation $\mathcal{I}$ and the
weighted adjacency matrix $\mathbf{W}$ of the meta-graph appropriately, we can
express various models from the literature as special cases of the CVN
graphical model.

#### Time-Varying Graphical LASSO (TVGL)

Suppose we observe graphs at $T$ time points, i.e.,
$\mathcal{U}=\\{t_{1},t_{2},\ldots,T_{T}\\}$. We define the index set
$\mathcal{I}=\\{1,2,\ldots,T\\}$ where $i\in\mathcal{I}$ corresponds to time
point $t_{i}$. Consequently, the meta-graph $\mathcal{G}_{\text{TVGL}}$ has
$T$ nodes, each corresponding to a single graph, and only graphs on
consecutive time-points are smoothed. We define the $(T\times T)$-dimensional
adjacency matrix of $\mathcal{G}_{\text{TVGL}}$ as

$\mathbf{W}_{\text{TVGL}}=\left(w_{ij}^{\text{TVGL}}\right)_{T\times
T}\qquad\text{where }w_{ij}^{\text{TVGL}}=\begin{cases}1&\text{if }|i-j|=1\\\
0&\text{otherwise}.\end{cases}$ (6)

The time-varying graphical LASSO
$\\{\widehat{\mathbf{\Theta}}_{t}^{\text{TVGL}}\\}_{t=1}^{T}$ can then be
estimated by solving

$\\{\widehat{\mathbf{\Theta}}_{t}^{\text{TVGL}}\\}_{t=1}^{T}=\operatorname*{arg\,min}_{\\{\mathbf{\Theta}_{t}\\}_{t=1}^{T}}\sum_{t=1}^{T}\ell_{t}(\mathbf{\Theta}_{t})+\lambda_{1}\sum_{t=1}^{T}\left\lVert\mathbf{\Theta}_{t}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{t=2}^{T}\left\lVert\mathbf{\Theta}_{t}-\mathbf{\Theta}_{t-1}\right\rVert_{1}^{\text{off}}$
(7)

which is equivalent to the definition in Hallac et al., (2017) and Monti et
al., (2014). Figure 1 shows an example of a TVGL in which each of the $T$
graphs has five nodes corresponding to $X_{1},X_{2},\ldots,X_{5}$.

Figure 1: The Time-Varying Graphical LASSO (TVGL) as a special case of the CVN
model, see eq. (6) and (7). A total of $T$ graphs $G_{i}$ with $p=5$ variables
$X_{1},X_{2},\ldots,X_{5}$ are observed at time points
$t_{1},t_{2},\ldots,t_{T}$. The circle around each graph represents a node in
the meta-graph $\mathcal{G}_{\text{TVGL}}$; the links between the graphs are
the meta-graph’s edges. Graphs on consecutive time points are smoothed with
weight $1$.

#### Fused Graphical LASSO (FGL)

Figure 2: The joint graphical LASSO with the fused graphical LASSO (FGL)
penalty term as a special case of the CVN model, see eq. (8) and (9). A total
of $m$ graphs with $p=5$ variables are observed. Just as in Figure 1, each
graph represents a node in the meta-graph $\mathcal{G}_{\text{FGL}}$ and links
between the graphs are the meta-graph’s edges.

When there is only one external covariate, the CVN reduces in some cases to
the fused graphical LASSO (FGL) from Danaher et al., (2014). The FGL takes
advantage of structural similarities across different groups to estimate $m$
graphical models that share common edges. Let $m$ be the number of categories
of the single external variable $U$, and for each category there are
$n_{1},n_{2},\ldots,n_{m}>0$ observations. As before, we choose the index set
$\mathcal{I}=\\{1,2,\ldots,m\\}$ where the graph of $i$-th category is denoted
with $G_{i}$. The CVN reduces to the FGL with the fused graphical LASSO
penalty term if the $(m\times m)$-dimensional weight matrix is defined as

$\mathbf{W}_{\text{FGL}}=(w_{ij}^{\text{FGL}})_{m\times m}\qquad\text{where
}w_{ij}^{\text{FGL}}=\begin{cases}1&\text{if }i\neq j\\\
0&\text{otherwise}.\end{cases}$ (8)

This leads to the following minimization problem:

$\\{\widehat{\mathbf{\Theta}}^{\text{FGL}}_{i}\\}_{i=1}^{m}=\operatorname*{arg\,min}_{\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}}\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{i<j}\left\lVert\mathbf{\Theta}_{i}-\mathbf{\Theta}_{j}\right\rVert_{1}^{\text{off}},$
(9)

which is the FGL (Danaher et al.,, 2014). Expressed as CVN this means that the
meta-graph is fully connected with weights of $1$. The FGL penalty enforces
that every graph is smoothed with every other graph. See Figure 2 for a visual
representation of this model.

#### Multiple Static Graphical Models

A trivial case is when no structural constraints are imposed on the $m$
graphs. The weight matrix is hence defined as $\mathbf{W}_{0}=\mathbf{0}$. The
CVN reduces then to $m$ independent GLASSO problems as in eq. (2):

$\\{\widehat{\mathbf{\Theta}}_{i}^{0}\\}_{i=1}^{m}=\operatorname*{arg\,min}_{\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}}\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}=\left\\{\operatorname*{arg\,min}_{\mathbf{\Theta}_{i}}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}\right\\}_{i=1}^{m}.$

Figure 3 shows a visual representation.

##  4 An Alternating Direction Method of Multipliers Algorithm for the
Covariate-Varying Network Problem

Estimating the CVN graphical model defined in eq. (3) requires us to solve the
optimization problem in eq. (5). We propose to do this using an Alternating
Direction Method of Multipliers (ADMM) algorithm (Boyd and Vandenberghe,,
2004). By doing so, we can split the original problem into three subproblems,
which in turn can be solved either analytically or numerically.

Figure 3: A trivial case with $\mathbf{W}_{0}=\mathbf{0}$ without smoothing
the $m$ graphs.

First, let us introduce a set of _consensus variables_
$\mathbf{Z}=\\{\mathbf{Z}_{i}\\}_{i=1}^{m}$ where
$\mathbf{Z}_{i}=\mathbf{\Theta}_{i}$ for all $i$. In addition, we write
$\mathbf{\Theta}=\\{\mathbf{\Theta}_{i}\\}_{i=1}^{m}$. We can then express eq.
(5) as the equality-constrained convex optimization problem:

$\begin{gathered}\operatorname*{arg\,min}_{\mathbf{\Theta},\mathbf{Z}}\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{Z}_{i}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{i<j}w_{ij}\left\lVert\mathbf{Z}_{i}-\mathbf{Z}_{j}\right\rVert_{1}^{\text{off}}\\\\[4.30554pt]
\text{subject to }\mathbf{\Theta}_{i}-\mathbf{Z}_{i}=\mathbf{0}\text{ and
}\mathbf{\Theta}_{i},\mathbf{Z}_{i}\succ 0\text{ for
}i=1,2,\ldots,m.\end{gathered}$ (10)

Thus, the problem does not have to be solved jointly but can be solved for
$\mathbf{Z}$ and $\mathbf{\Theta}$ separately. The corresponding augmented
Lagrangian is

$\begin{split}\mathcal{L}_{\rho}(\mathbf{\Theta},\mathbf{Z},\mathbf{Y})&=\sum_{i=1}^{m}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{Z}_{i}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{i<j}w_{ij}\left\lVert\mathbf{Z}_{i}-\mathbf{Z}_{j}\right\rVert_{1}^{\text{off}}\\\
&+(\rho/2)\sum_{i=1}^{m}\left[\left\lVert\mathbf{\Theta}_{i}-\mathbf{Z}_{i}+\mathbf{Y}_{i}\right\rVert_{F}^{2}-\left\lVert\mathbf{Y}_{i}\right\rVert_{F}^{2}\right],\end{split}$
(11)

where $\mathbf{Y}=\\{\mathbf{Y}_{i}\\}_{i=1}^{m}$ is a set of _scaled dual
variables_ ,
$\left\lVert\mathbf{A}\right\rVert_{F}^{2}=\left(\sum_{s,t}a_{st}^{2}\right)^{1/2}$
is the Frobenius norm, and $\rho>0$ is the ADMM’s _dual update step length_
(Boyd and Vandenberghe,, 2004; Danaher et al.,, 2014; Hallac et al.,, 2017).
Minimizing the function in eq. (11) is equivalent to minimizing the original
problem in eq. (5). The advantage, however, of this alternative definition, is
that eq. (10) can now be solved by iteratively updating $\mathbf{\Theta}$,
$\mathbf{Z}$ and $\mathbf{Y}$ (Boyd and Vandenberghe,, 2004). The update steps
for the ADMM are

Likelihood update:
$\displaystyle\mathbf{\Theta}(k+1)=\operatorname*{arg\,min}_{\mathbf{\Theta}}\leavevmode\nobreak\
\mathcal{L}_{\rho}(\mathbf{\Theta},\mathbf{Z}(k),\mathbf{Y}(k)),$ (12)
Constraint update:
$\displaystyle\mathbf{Z}(k+1)=\operatorname*{arg\,min}_{\mathbf{Z}}\leavevmode\nobreak\
\mathcal{L}_{\rho}(\mathbf{\Theta}(k+1),\mathbf{Z},\mathbf{Y}(k))\text{ and }$
(13) Dual update:
$\displaystyle\mathbf{Y}(k+1)=\operatorname*{arg\,min}_{\mathbf{Y}}\leavevmode\nobreak\
\mathcal{L}_{\rho}(\mathbf{\Theta}(k+1),\mathbf{Z}(k+1),\mathbf{Y}),$

where $k$ is the iteration step count. We find that the update step for
$\mathbf{Y}(k+1)$ has the analytical solution

$\mathbf{Y}_{i}(k+1)=\mathbf{Y}_{i}(k)+\left[\mathbf{\Theta}_{i}(k+1)-\mathbf{Z}_{i}(k+1)\right]\quad\text{for
}i=1,2,\ldots,m.$ (14)

The update steps for $\mathbf{\Theta}(k+1)$ and $\mathbf{Z}(k+1)$ have to be
solved numerically. The first has a known solution (Danaher et al.,, 2014).
The latter is more involved and is discussed in Section 4.2.

The algorithm iterates until convergence is reached. To guarantee convergence
to an optimal solution, the constraint
$\mathbf{\Theta}_{i}-\mathbf{Z}_{i}=\mathbf{0}$ must be satisfied while the
augmented Lagrangian forms (12) and (13) have to be minimized. In each
iteration, the minimization constraint on the dual variable is satisfied,
which ensures that the constraint update is satisfied by construction (Boyd et
al.,, 2010; Monti et al.,, 2014). The likelihood update step holds
asymptotically as $k\rightarrow\infty$ using the following stopping criterion
(Danaher et al.,, 2014; Wu et al.,, 2019):

$\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}(k+1)-\mathbf{\Theta}_{i}(k)\right\rVert_{1}/\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}(k)\right\rVert_{1}<\epsilon.$
(15)

As a default, we choose $\epsilon=10^{-5}$ as our convergence threshold and
control the step size by setting $\rho=1$.

###  4.1 Update Step for $\mathbf{\Theta}$

The update step for the precision matrices $\mathbf{\Theta}(k+1)$ in eq. (12)
can be divided into $m$ distinct problems, one for each precision matrix:

$\begin{split}\mathbf{\Theta}_{i}(k+1)&=\operatorname*{arg\,min}_{\mathbf{\Theta}_{i}}\leavevmode\nobreak\
\mathcal{L}_{\rho}(\mathbf{\Theta}_{i},\mathbf{Z}_{i}(k),\mathbf{Y}_{i}(k))\\\
&=\operatorname*{arg\,min}_{\mathbf{\Theta}_{i}}\leavevmode\nobreak\
\frac{n_{i}}{2}\left[\text{trace}\left(\widehat{\mathbf{\Sigma}}_{i}\mathbf{\Theta}_{i}\right)-\log\text{det}\left(\mathbf{\Sigma}_{i}\right)\right]+\\\
&\qquad(\rho/2)\left\lVert\mathbf{Y}_{i}(k)+\left[\mathbf{\Theta}_{i}-\mathbf{Z}_{i}(k)\right]\right\rVert_{F}^{2}\end{split}$
(16)

for $i=1,2,\ldots,m$. The precision matrix $\mathbf{\Theta}_{i}(k+1)$ can be
determined numerically. Due to $\rho>0$, one must choose $\mathbf{\Theta}_{i}$
in such a way that $\mathbf{\Theta}_{i}$ both minimizes the log-likelihood and
is in the proximity of $\mathbf{Z}_{i}(k)$. The degree to which this must be
enforced depends on both $\rho$ and $\mathbf{Y}_{i}(k)$. Others have shown
that eq. (16) can be solved numerically by rescaling an eigenvalue
decomposition (see Witten and Tibshirani, (2009) and Monti et al., (2014)).
More precisely, let $\mathbf{Q}_{i}\mathbf{\Gamma}_{i}\mathbf{Q}_{i}^{\top}$
be the eigenvalue decomposition of
$\widehat{\mathbf{\Sigma}}_{i}-(\rho/n_{i})\left[\mathbf{Z}_{i}(k)-\mathbf{Y}_{i}(k)\right]$,
where $\mathbf{Q}_{i}$ the columns are the eigenvectors and
$\mathbf{\Gamma}_{i}=\text{diag}(\gamma_{1},\gamma_{2},\ldots,\gamma_{p})$ is
a diagonal matrix with the corresponding eigenvalues on the diagonal. We
rescale this matrix such that
$\widetilde{\mathbf{\Gamma}}_{i}=\text{diag}(\widetilde{\gamma}_{i1},\widetilde{\gamma}_{i2},\ldots,\widetilde{\gamma}_{ip})$
where
$\widetilde{\gamma}_{ij}=\dfrac{n_{i}}{2\rho}\left((\gamma_{ij}^{2}+4\rho/n_{i})^{\tfrac{1}{2}}-\gamma_{ij}\right)$
for $j=1,2,\ldots,p$. The update step from eq. (16) then equals

$\mathbf{\Theta}_{i}(k+1)=\mathbf{Q}_{i}\widetilde{\mathbf{\Gamma}}_{i}\mathbf{Q}_{i}^{\top}.$
(17)

Note that each $\mathbf{\Theta}_{i}(k+1)$ is symmetric. Furthermore, all
$\widetilde{\gamma}_{i}$ are strictly positive, ensuring that each
$\mathbf{\Theta}_{i}(k+1)$ is positive definite.

The computationally most challenging part of the update step is the eigenvalue
decomposition of a $(p\times p)$-matrix with a complexity of
$\mathcal{O}(p^{3})$. Since this is needed for each of the $m$ graphs, the
$\mathbf{\Theta}$-update step’s complexity is $\mathcal{O}(mp^{3})$.

###  4.2 Update Step for $\mathbf{Z}$

The minimization of the augmented Lagragian with respect to $\mathbf{Z}$ in
eq. (13) is equal to

$\begin{split}\mathbf{Z}(k+1)&=\operatorname*{arg\,min}_{\mathbf{Z}}\mathcal{L}_{\rho}(\mathbf{\Theta}(k+1),\mathbf{Z},\mathbf{Y}(k+1))\\\
&=\operatorname*{arg\,min}_{\mathbf{Z}}\leavevmode\nobreak\
(\rho/2)\sum_{i=1}^{m}\left\lVert\mathbf{Z}_{i}-\mathbf{\Theta}_{i}(k+1)-\mathbf{Y}_{i}(k)\right\rVert_{F}^{2}+\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{Z}_{i}\right\rVert_{1}^{\text{off}}\\\
&\qquad\qquad\qquad+\lambda_{2}\sum_{i<j}w_{ij}\left\lVert\mathbf{Z}_{i}-\mathbf{Z}_{j}\right\rVert_{1}^{\text{off}}.\end{split}$
(18)

In contrast to $\mathbf{\Theta}(k+1)$, for which we could divide the update
step into $m$ distinct problems, we need to solve the problem for
$\mathbf{Z}(k+1)$ jointly due to the smoothness term
$\left\lVert\mathbf{Z}_{i}-\mathbf{Z}_{j}\right\rVert_{1}^{\text{off}}$.
However, we can divide the problem into multiple, independent optimization
problems, one for each node pair $\\{s,t\\}$, corresponding to a potential
edge in the graphs. Let $\\{s,t\\}$ be the potential edge between nodes
$X_{s}$ and $X_{t}$. Eq. (18) can then be written as

$\begin{split}\mathbf{Z}(k+1)&=\operatorname*{arg\,min}_{\mathbf{Z}}\leavevmode\nobreak\
(\rho/2)\sum_{i=1}^{m}\sum_{s,t}\left[z_{st}^{(i)}-\theta_{st}^{(i)}(k+1)-y_{st}^{(i)}(k)\right]^{2}\\\
&\qquad\qquad\qquad+\lambda_{1}\sum_{i=1}^{m}\sum_{s\neq
t}|z_{st}^{(i)}|+\lambda_{2}\sum_{i<j}\sum_{s\neq
t}w_{ij}|z^{(i)}_{st}-z^{(j)}_{st}|,\end{split}$

which is completely separable with respect to each node pairs $\\{s,t\\}$. Let

$\displaystyle\bm{\theta}_{st}$
$\displaystyle=\left(\theta_{st}^{(1)},\theta_{st}^{(2)},\ldots,\theta_{st}^{(m)}\right)^{\top},$
$\displaystyle\bm{z}_{st}$
$\displaystyle=\left(z_{st}^{(1)},z_{st}^{(2)},\ldots,z_{st}^{(m)}\right)^{\top},\quad\text{and}$
$\displaystyle\bm{y}_{st}$
$\displaystyle=\left(y_{st}^{(1)},y_{st}^{(2)},\ldots,y_{st}^{(m)}\right)^{\top}$

be $m$-dimensional vectors with the $(s,t)$ entries of the matrices
$\mathbf{\Theta}$, $\mathbf{Z}$ and $\mathbf{Y}$, respectively. There is an
analytical solution for the diagonal elements of $\bm{Z}(k+1)$, specifically:

$\bm{z}_{ss}(k+1)=\operatorname*{arg\,min}_{\bm{z}_{ss}\in\mathds{R}^{m}}\leavevmode\nobreak\
\left\lVert\bm{z}_{ss}-\bm{\theta}_{ss}(k+1)-\bm{y}_{ss}(k)\right\rVert_{2}^{2}=\bm{\theta}_{ss}(k+1)+\bm{y}_{ss}(k).$

For $s\neq t$, we can update each $(s,t)$ entry by solving

$\begin{split}\bm{z}_{st}(k+1)&=\operatorname*{arg\,min}_{\bm{z}_{st}\in\mathds{R}^{m}}\leavevmode\nobreak\
(\rho/2)\left\lVert\bm{z}_{st}-\bm{\theta}_{st}(k+1)-\bm{y}_{st}(k)\right\rVert_{2}^{2}\\\
&\qquad\qquad+\lambda_{1}\left\lVert\bm{z}_{st}\right\rVert_{1}+\lambda_{2}\sum_{i<j}w_{ij}|z_{st}^{(i)}-z_{st}^{(j)}|.\end{split}$

We can express this problem as a _weighted Fused LASSO Signal Approximator_
(wFLSA). For ease of notation, let us define $\bm{\beta}=\bm{z}_{st}$ and
${\bm{y}=\bm{\theta}_{st}(k+1)+\bm{y}_{st}(k)}$; we then find that our problem
can be written as

$\bm{\beta}(k+1)=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{m}}\leavevmode\nobreak\
\frac{1}{2}\left\lVert\bm{y}-\bm{\beta}\right\rVert_{2}^{2}+\eta_{1}\left\lVert\bm{\beta}\right\rVert_{1}+\eta_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|$
(19)

where $\eta_{1}=\lambda_{1}/\rho$ and $\eta_{2}=\lambda_{2}/\rho$. Since the
wFLSA was not solved before, we developed an ADMM algorithm specifically to
this end, see Dijkstra et al., (2024). The complexity of this algorithm is in
the order of $\mathcal{O}(m^{2})$.

Figure 4: A graphical representation of the algorithm for estimating the CVN
model, see eq. (5). The input is the raw data set and the weight matrix; the
output is the estimated CVN model. The ADMM consists of three updates steps
for $\mathbf{\Theta}$, $\mathbf{Z}$ and $\mathbf{Y}$ that are repeated till
convergence. Update step $\mathbf{\Theta}$ and $\mathbf{Y}$ can be subdivided
into $m$ independent step. During the $\mathbf{Z}$-update step, $\binom{p}{2}$
wFLSA problems need to be solved, one for each potential edge.

###  4.3 Pseudo-algorithm for the CVN model

Combining the results presented in the previous section, leads us to the
following procedure for estimating CVN models. A graphical representation can
be found in Figure 4.

  1. 1.

Compute the empirical covariance matrices
$\widetilde{\mathbf{\Sigma}}_{i}=(\sigma_{st}^{(i)})_{p\times
p}=n_{i}^{-1}\sum_{j=1}^{n_{i}}\bm{x}_{j}^{(i)}\left(\bm{x}_{j}^{(i)}\right)^{\top}$
with $i=1,2,\ldots,m$.

  2. 2.

Initialize
$\mathbf{\Theta}_{i}(1)=\text{diag}(\widetilde{\mathbf{\Sigma}}_{i})^{-1}$ and
$\mathbf{Z}_{i}(1)=\mathbf{Y}_{i}=\mathbf{0}$. Alternatively, one can use a
warmstart where
$\mathbf{\Theta}_{i}(1)=\operatorname*{arg\,min}_{\mathbf{\Theta}_{i}}\ell_{i}(\mathbf{\Theta}_{i})+\lambda_{1}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}$
is the GLASSO estimate for the individual graphs $i=1,2,\ldots,m$.

  3. 3.

Update $\mathbf{\Theta}$. For each graph $i=1,2,\ldots,m$, do:

    1. (a)

Compute the eigenvalue decomposition
$\mathbf{Q}_{i}\mathbf{\Gamma}_{i}\mathbf{Q_{i}}^{\top}$ of
$\widetilde{\mathbf{\Sigma}}_{i}-(\rho/n_{i})\left[\mathbf{Z}_{i}(k)-\mathbf{Y}_{i}(k)\right]$;

    2. (b)

Let
$\widetilde{\mathbf{\Gamma}}_{i}=\text{diag}(\widetilde{\gamma}_{i1},\widetilde{\gamma}_{i2},\ldots,\widetilde{\gamma}_{ip})$
where
$\widetilde{\gamma}_{ij}=\frac{n_{i}}{2\rho}\left(\sqrt{\gamma_{ij}^{2}+4\rho/n_{i}}-\gamma_{ij}\right)$;

    3. (c)

Then
$\mathbf{\Theta}_{i}(k+1)=\mathbf{Q}_{i}\widetilde{\mathbf{\Gamma}}_{i}\mathbf{Q}_{i}^{\top}$.

  4. 4.

Update $\mathbf{Z}$ by going over all potential edges $\\{s,t\\}$:

    1. (a)

Let
$\bm{\theta}_{st}=\left(\theta_{st}^{(1)},\theta_{st}^{(2)},\ldots,\theta_{st}^{(m)}\right)^{\top}$,
$\bm{z}_{st}=\left(z_{st}^{(1)},z_{st}^{(2)},\ldots,z_{st}^{(m)}\right)^{\top}$
and
$\bm{y}_{st}=\left(y_{st}^{(1)},y_{st}^{(2)},\ldots,y_{st}^{(m)}\right)^{\top}$;

    2. (b)

Update the diagonal: $\bm{z}_{ss}(k+1)=\bm{\theta}_{ss}(k+1)+\bm{y}_{ss}(k)$
for $s=1,2,\ldots,p$;

    3. (c)

Update the off-diagonal entries $\bm{z}_{st}(k+1)$ and $\bm{z}_{ts}(k+1)$ for
all $s<t$. Let $\bm{\beta}=\bm{z}_{st}(k+1)$ and
$\bm{y}=\bm{\theta}_{st}(k+1)+\bm{y}_{st}(k)$, then

$\bm{z}_{st}(k+1)=\operatorname*{arg\,min}_{\bm{\beta}\in\mathds{R}^{m}}\leavevmode\nobreak\
\frac{1}{2}\left\lVert\bm{y}-\bm{\beta}\right\rVert_{2}^{2}+\eta_{1}\left\lVert\bm{\beta}\right\rVert_{1}+\eta_{2}\sum_{i<j}w_{ij}|\beta_{i}-\beta_{j}|$

is the corresponding wFLSA problem, where $\eta_{1}=\lambda_{1}/\rho$ and
$\eta_{2}=\lambda_{2}/\rho$. Solve this using the algorithm presented in
Dijkstra et al., (2024).

  5. 5.

Update $\mathbf{Y}$:
$\mathbf{Y}_{i}(k+1)=\mathbf{\Theta}_{i}(k+1)-\mathbf{Z}_{i}(k+1)+\mathbf{Y}_{i}(k)$
for $i=1,2,\ldots,m$.

  6. 6.

In case the stopping criterion
$\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}(k+1)-\mathbf{\Theta}_{i}(k)\right\rVert_{1}/\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}(k)\right\rVert_{1}<\epsilon$
has been met, return the estimate of the CVN:
$\widehat{\mathbf{\Theta}}_{\text{CVN}}=\\{\mathbf{\Theta}_{i}(k+1)\\}_{i=1}^{m}$.
Otherwise, repeat steps 3 to 5.

An implementation of this algorithm is publicly available as an R package at
www.github.com/bips-hb/CVN.

##  5 Algorithmic Complexity

Here, we determine the computational complexity of the CVN algorithm presented
in the previous section. We first consider each of the update steps separately
and combine the results at the end.

The $\mathbf{\Theta}(k+1)$-update step can be decomposed into $m$ distinct
steps, each associated with a specific graph $i=1,2,\ldots,m$, as illustrated
in Figure 4. The most computationally demanding aspect of the
$\mathbf{\Theta}_{i}(k+1)$ update step, see eq. (12), involves performing the
eigendecomposition of a $(p\times p)$-dimensional matrix, a task known to be
on the order of $\mathcal{O}(p^{3})$. The entire update step requires, thus,
in the order of $\mathcal{O}(mp^{3})$ operations.

In the case of the $\mathbf{Y}(k+1)$ update step, see eq. 14, the task entails
summing three $(p\times p)$-dimensional matrices for each graph
$i=1,2,\ldots,m$. This operation, thus, requires $\mathcal{O}(mp^{2})$
operations.

The $\mathbf{Z}(k+1)$-update step requires one to solve a wFLSA problem for
each of the $\binom{p}{2}$ potential edges in the graph. We know that the
complexity of solving a single wFLSA problem is $\mathcal{O}(rm^{2})$, where
$r$ is the (average) number of iterations needed for the algorithm to converge
(Dijkstra et al.,, 2024). The overall complexity of the
$\mathbf{Z}(k+1)$-update step is, thus, $\mathcal{O}\left(\binom{p}{2}\cdot
rm^{2}\right)=\mathcal{O}(rm^{2}p^{2})$.

Combining these results, yields the computational complexity of the ADMM
algorithm for solving CVN models as

$\mathcal{O}\Big{(}R\cdot\big{(}\underbrace{mp^{3}}_{\mathbf{\Theta}(k+1)\text{-update}}+\underbrace{mp^{2}}_{\mathbf{Y}(k+1)\text{-update}}+\underbrace{rm^{2}p^{2}}_{\mathbf{Z}(k+1)\text{-update}}\big{)}\Big{)}=\mathcal{O}\left(Rrm^{2}p^{3}\right),$

where $R$ is the number of iterations needed for the ADMM to satisfy the
stopping criterion in eq. (15), $r$ denotes the average iterations needed for
the wFLSA algorithm to converge (Dijkstra et al.,, 2024), $m$ is the number of
graphs, and $p$ is the number of variables.

##  6 Interpolation

In this section, we propose a method for interpolating a graph for which there
are no observations based on an estimated CVN model with $m$ graphs. We denote
the graph we want to interpolate by $G_{m+1}$. We are concerned solely with
the existence of edges and not values of the entries in the precision matrix
$\mathbf{\Theta}_{m+1}$. Approaches in the literature often rely on the so-
called ‘AND’ and ‘OR’ rules (see, for example, Hallac et al., (2017)). In case
of the former, the interpolated graph has the edge $\\{s,t\\}$ if and only if
_all_ $m$ graphs in the CVN model have the edge $\\{s,t\\}$. In case of the
‘OR’ rule, the interpolated graph has the edge $\\{s,t\\}$ if at least one of
the other estimated graphs has that edge. Note that with this approach, it is
difficult to control the level of sparsity and take different levels of
similarity between the graphs into account.

We interpolate graph $G_{m+1}$ given the estimated CVN model
$\widehat{\mathbf{\Theta}}_{\text{CVN}}$ by solving

$\widehat{\mathbf{\Theta}}_{m+1}=\operatorname*{arg\,min}_{\mathbf{\Theta}_{m+1}\succ
0}\ell_{m+1}(\mathbf{\Theta}_{m+1})+\lambda_{1}\left\lVert\mathbf{\Theta}_{m+1}\right\rVert_{1}^{\text{off}}+\lambda_{2}\sum_{i=1}^{m}\omega_{i}\left\lVert\widehat{\mathbf{\Theta}}_{i}-\mathbf{\Theta}_{m+1}\right\rVert_{1}^{\text{off}}$

where $\bm{\omega}=(\omega_{1},\omega_{2},\ldots,\omega_{m})^{\top}$ are the
_smoothing coefficients_ which reflect the level of similarity between the
individual graphs of the CVN and the interpolated graph. Note that, since
there are no observations for the graph $G_{m+1}$, the corresponding log-
likelihood is equal to $\ell_{m+1}(\mathbf{\Theta}_{m+1})=0$, see eq. (4).
Similar to the $\mathbf{Z}$-update step, we can divide the problem into
$\binom{p}{2}$ separate optimization problems, one for each potential edge
$\\{s,t\\}$. In fact, we can write the function as

$\begin{split}\widehat{\mathbf{\Theta}}_{m+1}&=\operatorname*{arg\,min}_{\mathbf{\Theta}_{m+1}\succ
0}\
\lambda_{1}\sum_{s,t}|\theta_{st}^{(m+1)}|+\lambda_{2}\sum_{i=1}^{m}\omega_{i}\sum_{s,t}|\widehat{\theta}_{st}^{(i)}-\theta_{st}^{(m+1)}|\\\
&=\operatorname*{arg\,min}_{\mathbf{\Theta}_{m+1}\succ 0}\
\sum_{s,t}\left[\lambda_{1}|\theta_{st}^{(m+1)}|+\lambda_{2}\sum_{i=1}^{m}w_{i}|\widehat{\theta}_{st}^{(i)}-\theta_{st}^{(m+1)}|\right].\end{split}$

We, thus, need to solve for each potential edge $\\{s,t\\}$:

$\widehat{\theta}_{st}^{(m+1)}=\operatorname*{arg\,min}_{\theta_{st}^{(m+1)}\in\mathds{R}}\lambda_{1}|\theta_{st}^{(m+1)}|+\lambda_{2}\lVert\bm{\omega}^{\top}\widehat{\bm{\theta}}_{st}-\theta_{st}^{(m+1)}\rVert_{1}.$

This function is convex (see Appendix A) and can be solved using any
derivative-free search method. In this context, we employ Brent’s algorithm
(Brent,, 2013). The sign of $\widehat{\theta}_{st}^{(m+1)}$ is
inconsequential, as our primary objective is to establish the presence or
absence of the edge $\\{s,t\\}$ in the interpolated graph, i.e.,
$\widehat{\theta}_{st}^{(m+1)}\neq 0$. Note that the smoothing coefficients
$\omega_{1},\omega_{2},\ldots,\omega_{m}$ need to be chosen a priori, which
requires a certain level of prior knowledge.

##  7 Tuning Parameter Selection

After fitting the CVN model, the next step involves selecting suitable values
for the tuning parameters: $\lambda_{1}$, which governs the sparsity of the
graphs, and $\lambda_{2}$, responsible for regulating the smoothness or
similarity between the graphs. Since having prior knowledge about appropriate
values for the tuning parameters is rare, one often opts for a data-driven
approach. Numerous methods exist for tuning parameter selection, and in this
context, we explore two frequently used measures: the Akaike Information
Criterion (AIC) and Bayesian Information Criterion (BIC). However, it is
important to note that in our scenario, these criteria can only be
approximated, see Danaher et al., (2014).

Let $\widehat{\mathbf{\Theta}}_{i}(\lambda_{1},\lambda_{2})$ denote the
estimator of the precision matrix for the $i$-th graph for specific values of
$\lambda_{1}$ and $\lambda_{2}$. The AIC can be approximated as

$\text{AIC}(\lambda_{1},\lambda_{2})=\sum_{i=1}^{m}\ell_{i}(\widehat{\mathbf{\Theta}}_{i}(\lambda_{1},\lambda_{2}))+2\lVert\widehat{\mathbf{\Theta}}_{i}(\lambda_{1},\lambda_{2})\rVert_{0},$
(20)

where $\left\lVert\cdot\right\rVert_{0}$ denotes the $L_{0}$-norm, equivalent
to the number of non-zero entries in the matrix (Danaher et al.,, 2014). On
the other hand, the BIC, taking into account the number of observations for
the graphs, tends to result in sparser models than the AIC:

$\text{BIC}(\lambda_{1},\lambda_{2})=\sum_{i=1}^{m}\ell_{i}(\widehat{\mathbf{\Theta}}_{i}(\lambda_{1},\lambda_{2}))+2\log(n_{i})\lVert\widehat{\mathbf{\Theta}}_{i}(\lambda_{1},\lambda_{2})\rVert_{0},$
(21)

where $n_{i}$ is the number of observations for the $i$-th graph. The optimal
values for $\lambda_{1}$ and $\lambda_{2}$ are those that minimize either the
approximated AIC or BIC. In our simulation study, see Section 9, we consider
both information criteria. Note that while AIC and BIC are commonly used, they
are not the sole metrics for selecting the optimal tuning parameters based on
the data. Other methods are cross-validation, bootstrap (Hastie et al.,,
2009), and alternative information criteria like the extended BIC (Chen and
Chen,, 2008) and Predictive Information Criteria (Ando,, 2007). Here, we opt
for AIC and BIC due to the computational efficiency, which is advantageous for
a computationally intensive algorithm such as the ADMM propsed here.

##  8 An Alternative Tuning Parameterization

Due to the computational complexity of the algorithm exploring an extensive
range of values for the tuning parameters $\lambda_{1}$ and $\lambda_{2}$ is
impractical, and, as mentioned above, obtaining expert or prior knowledge
about ‘optimal’ values for these parameters is challenging.

We introduce an alternative parameterization for the CVN model. Instead of
using $\lambda_{1}$ and $\lambda_{2}$, we propose employing two different
parameters, $\gamma_{1}$ and $\gamma_{2}$, which we define in this section.
The benefit of this new parameterization is that the interpretation of these
values do _not_ heavily depend on the number of graphs ($m$) and/or the number
of variables ($p$).

Suppose a model has been fitted and the best values for the tuning parameters
$\lambda_{1}$ and $\lambda_{2}$ have already been selected based on some
criterion, see the previous section. Now, consider a scenario where the same
data set is used, but a number of variables are added or removed (thereby
changing $p$), and/or other external covariates are added or removed (changing
the number of graphs $m$). Due to the way the model is defined, see eq. (5),
the interpretation of the penalty terms $\lambda_{1}$ and $\lambda_{2}$
change, i.e., the optimal values for the first model provide little
information about the best values for the tuning parameters for the new model.

Note that this challenge is not unique to the CVN model; it extends to other
penalized regression problems as well. However, this challenge is typically
overlooked for two main reasons. Firstly, it is in other situations uncommon
for the considered variables to change. Secondly, the associated algorithm is
usually computationally light, making it straightforward to reapply to new
datasets. Nevertheless, we would like to point out that the alternative tuning
parameterization proposed here may be useful for other LASSO-related methods
as well.

In the existing parameterization, penalties are imposed on the sum of the
$L_{1}$-norm of the precision matrices, governed by $\lambda_{1}$, and the sum
of the $L_{1}$-norm of the differences between the precision matrices,
regulated by $\lambda_{2}$, see eq. (5). We propose a modification where we
penalize each individual edge in the CVN model to regulate sparsity, and each
pair of edges in the CVN model to control the level of smoothness. We denote
the new parameters as $\gamma_{1}$, regulating sparsity, and $\gamma_{2}$,
managing smoothness.

Recall that the total number of potential edges in the CVN model is $mp(p-1)$,
and the total number of pairs of potential edges is $m(m-1)p(p-1)$. We define
$\gamma_{1}$ and $\gamma_{2}$ as follows:

$\gamma_{1}=2\frac{\lambda_{1}}{mp(p-1)}\qquad\text{and}\qquad\gamma_{2}=4\frac{\lambda_{2}}{m(m-1)p(p-1)}.$

Consequently, the sparsity penalty term can be expressed as

$\lambda_{1}\sum_{i=1}^{m}\left\lVert\mathbf{\Theta}_{i}\right\rVert_{1}^{\text{off}}=\sum_{i=1}^{m}\sum_{s<t}\gamma_{1}|\theta_{st}^{(i)}|.$

In other words, $\gamma_{1}$ regulates how much an individual edge in a graph
is penalized. Consequently, alterations in the number of edges have a
relatively modest impact on this tuning parameter in comparison to
$\lambda_{1}$. Similarly, the smoothness penalty term of the CVN model can be
written in terms of $\gamma_{2}$:

$\lambda_{2}\sum_{i<j}\left\lVert\mathbf{\Theta}_{i}-\mathbf{\Theta}_{j}\right\rVert_{1}^{\text{off}}=\sum_{i<j}w_{ij}\sum_{s,t}\gamma_{2}|\theta_{st}^{(i)}-\theta_{st}^{(j)}|.$

Here, $\gamma_{2}$ is used to penalize a single pair of edges. Once more,
changing the number of graphs and edges has little effect on the
interpretation of this tuning parameter.

The idea that optimal values for $\gamma_{1}$ and $\gamma_{2}$ are consistent
across various datasets hinges on the assumption that both the sparsity of the
graphs and their smoothness will remain relatively stable across those
datasets.

##  9 Simulation

In this simulation study, we consider a CVN model with two external
covariates, $U_{1}$ and $U_{2}$, each of which takes values from the set
$\\{1,2,3\\}$ featuring nine graphs. The graph corresponding to $U_{1}=i$ and
$U_{2}=j$ is denoted as $G_{ij}=G(U_{1}=i,U_{2}=j)$, and its associated
adjacency, precision, and covariance matrices are similarly indexed as
$\mathbf{A}_{ij}$, $\mathbf{\Theta}_{ij}$ and $\mathbf{\Sigma}_{ij}$,
respectively. Our simulation setup allows for varying the number of
observations ($n$), the number of variables ($p$) and the the general graph
structure, specifically Erdös-Rényi and Barabasi-Albert graphs (Newman,, 2018)
in our case. Additionally, we can control the extent to which the edges in the
graph change with $U_{1}$ and $U_{2}$. Adjusting $n$ and $p$ is
straightforward, as we will see later. Our focus is first on how to employ a
specific general graph structure and how to regulate the extent with which the
graphs change with changes in the two external covariates.

Figure 5: Illustration of the CNV model with two external covariates, $U_{1}$
and $U_{2}$, as used in the simulation study. The values of $U_{1}$ and
$U_{2}$ are shown on the horizontal and vertical axes, respectively. The
‘starting graph’ $G_{11}$ is located in the lower-left corner, followed by
graph $G_{21}$ to its right, and so on, ending with graph $G_{33}$ is in the
upper-right corner.

We can represent our case as a grid of graphs, as illustrated in Figure 5,
with $U_{1}$ and $U_{2}$ positioned along the horizontal and vertical axes,
respectively. Graph $G_{11}$ is located in the bottom-left corner, followed by
graph $G_{21}$ to its right, and so forth, concluding with graph $G_{33}$ in
the upper-right corner. We designate graph $G_{11}$ as the _starting graph_ ;
it is randomly generated according to a specific graph model. In this context,
we consider two types:

  1. 1.

The Erdös-Rényi (ER) graph, which has a single parameter $\pi$, representing
the probability of the existence of any given edge. The expected density of
the graph, defined as the number of edges divided by the total number of
possible edges (i.e., $p(p-1)/2$), is expressed as (Newman,, 2018)

$\mathds{E}_{\text{ER}}\left[|E|\right]=\pi.$

Here, $|E|$ denotes the cardinality of the edge set, and
$\mathds{E}_{\text{ER}}$ signifies the expectation under the ER model. The
relevance of the expected density of the graph becomes apparent later;

  2. 2.

The Barabasi-Albert (BA) graph, which is a scale-free graph. The degrees of
the graph follow a power law distribution with an exponent of 3 (Newman,,
2018). The expected density of this graph type is approximately given by
(Pósfai and Barabasi,, 2016)

$\mathds{E}_{\text{BA}}\left[|E|\right]\approx\frac{p-1}{2},$ (22)

where $\mathds{E}_{\text{BA}}$ represents the expectation under the BA model.

(a) Dividing the starting graph $G_{11}$

(b) Changing the first subgraph

(c) Changing the second subgraph

(d) Combining the results

Figure 6:  An example illustrating the generation of the nine graphs used in
the simulation study. (a) The starting graph is divided into two subgraphs
(orange and blue). (b) The edges in the orange subgraph change with respect to
$U_{1}$. (c) The edges in the blue subgraph change with respect to $U_{2}$.
(d) Combining both modifications in the orange and blue subgraphs yields the
nine graphs. Note that the edges between the two subgraphs (depicted in gray)
remain unaltered.

As previously mentioned, we aim to control the extent to which graphs change
with $U_{1}$ and $U_{2}$. To achieve this, we partition the starting graph
into two equal-sized subgraphs. See, for an example, Figure 6(a), where the
first subgraph is shown in orange and the second in blue. In cases where the
number of variables $p$ is odd, the first and second subgraphs comprise
$\lceil{p/2}\rceil$ and $\lfloor{p/2}\rfloor$ nodes, respectively, where
$\lceil\cdot\rceil$ and $\lfloor\cdot\rfloor$ denote the ceiling and floor
functions. The adjacency matrix of the starting graph can be expressed in
terms of these subgraphs as follows:

$\mathbf{A}_{11}=\begin{bmatrix}\mathbf{A}_{1}^{1}&\mathbf{A}\\\
\mathbf{A}^{\top}&\mathbf{A}_{1}^{2}\end{bmatrix}.$ (23)

Here, $\mathbf{A}_{1}^{1}$ represents the adjacency matrix for the first
subgraph, and $\mathbf{A}_{1}^{2}$ is the matrix for the second subgraph.
Similarly, we divide the remaining eight graphs into subgraphs using the same
structure:

$\mathbf{A}_{ij}=\begin{bmatrix}\mathbf{A}_{i}^{1}&\mathbf{A}\\\
\mathbf{A}^{\top}&\mathbf{A}_{j}^{2}\end{bmatrix}$

for $i=1,2,3$ and $j=1,2,3$.

When varying $U_{1}$, only the subgraph corresponding to $\mathbf{A}_{i}^{1}$
undergoes changes, while the remainder of the graph remains unaltered.
Similarly, when varying $U_{2}$, only the matrix $\mathbf{A}_{j}^{2}$
undergoes changes. It is important to note that the edges between the two
subgraphs, represented here as $\mathbf{A}$, remain constant throughout and
form the _core_ of the network.

After sampling the starting graph $G_{11}$ and obtaining $\mathbf{A}_{1}^{1}$
and $\mathbf{A}_{1}^{2}$, we derive $\mathbf{A}_{2}^{1}$ by randomly adding
and removing $b_{1}$ edges. If there are no $b_{1}$ edges in the subgraph, all
existing edges change. Subsequently, we obtain $\mathbf{A}_{3}^{1}$ by
altering $b_{1}$ edges in $\mathbf{A}_{2}^{1}$. A similar process is followed
for the second subgraph: we obtain $\mathbf{A}_{2}^{2}$ by randomly adding and
removing $b_{2}$ edges in the matrix $\mathbf{A}_{1}^{2}$. The matrix
$\mathbf{A}_{3}^{2}$ is then created by modifying $b_{2}$ edges in the matrix
$\mathbf{A}_{2}^{2}$.

We visually depict this process in Figure 6. Figure 6(a) displays the starting
graph $G_{11}$ with the first subgraph depicted in orange and the second in
blue. Changes occur exclusively in the edges of the orange subgraph as $U_{1}$
varies, as illustrated in Figure 6(b). Similarly, only edges in the blue
subgraph undergo changes with variations in $U_{2}$, as shown in Figure 6(c).
Combining the changes made in the orange and blue subgraphs results in the
nine graphs depicted in Figure 6(d).

Focusing on the number of edges $b_{1}$ and $b_{2}$ that change with varying
either $U_{1}$ or $U_{2}$ may be challenging to interpret, particularly for
different values of $p$ or the number of edges. Therefore, we find it more
meaningful to think in terms of the percentage of edges in the graph that
change with $U_{1}$ and $U_{2}$. These percentages are denoted as $\pi_{1}$
and $\pi_{2}$. The number of edges changing with $U_{1}$ and $U_{2}$ are then
given by

$b_{1}=\text{round}\left(\pi_{1}\mathds{E}[|E|]\right)\qquad\text{and}\qquad
b_{2}=\text{round}\left(\pi_{2}\mathds{E}[|E|]\right),$ (24)

where $\text{round}(\cdot)$ rounds to the nearest integer. In other words,
$\pi_{1}$ and $\pi_{2}$ represent the percentage of edges the graph is
expected to have ($\mathds{E}[|E|]$) that change with $U_{1}$ and $U_{2}$,
respectively.

We translate the obtained adjacency matrices into nine precision matrices.
Subsequently, these precision matrices are used to generate data using a
normal distribution. We adopt a methodology similar to the one outlined in
(Liu and Wang,, 2017). We define
$\widetilde{\mathbf{\Theta}}_{ij}=v\mathbf{A}_{ij}$ as the adjacency matrix
for $U_{1}=i$ and $U_{2}=j$, multiplied by the scalar $v>0$. Let $\kappa_{ij}$
represent the smallest eigenvalue of the matrix
$\widetilde{\mathbf{\Theta}}_{ij}$. To ensure the precision matrix is positive
definite, we ‘increase’ the diagonal, specifically,

$\mathbf{\Theta}_{ij}=\widetilde{\mathbf{\Theta}}_{ij}+\text{diag}\left(|\kappa_{ij}|+.1+u\right),$
(25)

where $u>0$ is a scalar as well. Throughout our simulation, we maintain $v=.4$
and $u=.1$, consistent with the values employed by Liu and Wang, (2017). The
corresponding covariance matrix is simply the inverse of the precision matrix,
i.e., $\mathbf{\Sigma}_{ij}=\mathbf{\Theta}_{ij}^{-1}$, and is determined
numerically. The data for each graph, $\mathbf{X}_{ij}$, consists of $n$
independent samples from the $p$-dimensional normal distribution
$\mathcal{N}(\mathbf{0},\mathbf{\Sigma}_{ij})$.

In conclusion, our simulator follows the steps:

  1. 1.

Specify the number of observations ($n$), parameters ($p$), initial graph type
(ER or BA model), and the percentage of edges in the graphs changing in the
$U_{1}$ and $U_{2}$ directions ($\pi_{1}$ and $\pi_{2}$);

  2. 2.

Determine the number of changing edges $b_{1}$ and $b_{2}$ based on the
percentages $\pi_{1}$ and $\pi_{2}$ using equation (24);

  3. 3.

Generate the starting graph $G_{11}$ with $p$ vertices following the selected
model from step 1.

  4. 4.

Divide the starting graph $G_{11}$ into two subgraphs, see eq. (23) and Figure
6(a);

  5. 5.

Generate the adjacency matrix $A_{2}^{1}$ by randomly adding and removing
$b_{1}$ edges from $A_{1}^{1}$, see Figure 6(b);

  6. 6.

Generate $A_{3}^{1}$ by adding and removing $b_{1}$ edges from $A_{2}^{1}$
obtained in the previous step;

  7. 7.

Repeat the last two steps for the adjacency matrices $A_{2}^{2}$ and
$A_{3}^{2}$ by changing $b_{2}$ edges, see Figure 6(c);

  8. 8.

Steps 3–7 result in nine adjacency matrices, see Figure 6(d). Determine the
corresponding nine precision matrices using eq. (25);

  9. 9.

Numerically determine the covariance matrices (inverse of the precision
matrices) from the previous step, and

  10. 10.

Generate nine data sets ${\mathbf{X}_{ij}}$, where $\mathbf{X}_{ij}$ comprises
$n$ independent draws from the $p$-dimensional normal distribution
$\mathcal{N}(\mathbf{0},\mathbf{\Sigma}_{ij}$).

The simulator is publicly available as an R package at www.github.com/bips-
hb/CVNSim.

##  10 Simulation Set-up

In our simulation study, we explore various parameter settings outlined in
Table 1. We set the number of variables, $p$, to 100 and vary the number of
observations, $n$, from 100 to 200. Both scenarios are high-dimensional, as
the number of edges in each graph is $\binom{p}{2}=4{,}950$. As previously
mentioned, we consider both the ER and BA models. For the ER model, the
density $\pi$ is set to $10\%$. The sparsity of the BA model is constant, as
given in eq. (22). We introduce variability in the degree to which graphs
change in the $U_{1}$ and $U_{2}$ directions by setting $(\pi_{1},\pi_{2})$ to
$(.1,.1)$, $(.1,.2)$, and $(.2,.2)$.

Table 1: Parameter settings for the simulation study Description | Parameter | Value
---|---|---
Number of variables | $p$ | 100
Number of observations | $n$ | 100, 200
Density ER graph | $\pi$ | .1
Percentage edges changing with $U_{1}$ | $\pi_{1}$ | .1, .2
Percentage edges changing with $U_{2}$ | $\pi_{2}$ | .1, .2

Using the CVN algorithm necessitates the specification of $\lambda_{1}$ and
$\lambda_{2}$, or equivalently, $\gamma_{1}$ and $\gamma_{2}$ (see Section 8),
along with the weight matrix $\mathbf{W}$. In our simulation, we opt for
$(\gamma_{1},\gamma_{2})\in\Gamma\times\Gamma$, where
$\Gamma=\\{10^{-5},5\cdot 10^{-5},10^{-4},5\cdot 10^{-4},10^{-3},5\cdot
10^{-3}\\}$. We explore three distinct weight matrices:

  1. 1.

$\mathbf{W}_{0}=\bm{0}_{m\times m}$, the zero matrix. This corresponds to
applying the GLASSO to each of the nine graphs individually, i.e., there is no
smoothing;

  2. 2.

$\mathbf{W}_{\text{grid}}$, defined as

$\mathbf{W}_{\text{grid}}=\begin{bmatrix}0&1&0&1&0&0&0&0&0\\\
1&0&1&0&1&0&0&0&0\\\ 0&1&0&0&0&1&0&0&0\\\ 1&0&0&0&1&0&1&0&0\\\
0&1&0&1&0&1&0&1&0\\\ 0&0&1&0&1&0&0&0&1\\\ 0&0&0&1&0&0&0&1&0\\\
0&0&0&0&1&0&1&0&1\\\ 0&0&0&0&0&1&0&1&0\\\ \end{bmatrix}$ (26)

representing the scenario where each graph is smoothed with its ‘adjacent’
graphs, and

  3. 3.

$\mathbf{W}_{\text{full}}=\mathbf{1}_{m\times m}-\text{diag}(\mathbf{1})$, a
matrix of ones with zeros on its diagonal. In this case, each graph is equally
smoothed with every other graph.

See Figure 7 for a visual representation of the meta-graphs implied by these
three weight matrices. Given the simulation process for the nine graphs, the
second weight matrix, $\mathbf{W}_{\text{grid}}$, mirrors the true underlying
structure closest. We repeat the simulation for each simulation setting and
weight matrix 20 times.

As performance measure, we consider the $F_{1}$-score which we define here.
Recall from the previous section that the true adjacency matrix is denoted as
$\mathbf{A}_{ij}=\left(a_{st}^{(ij)}\right)_{p\times p}$, and let
$\widehat{\mathbf{A}}_{ij}(\gamma_{1},\gamma_{2})=\left(\widehat{a}_{st}^{(ij)}(\gamma_{1},\gamma_{2})\right)_{p\times
p}$ be the estimated adjacency matrix for $U_{1}=i$, $U_{2}=j$ and tuning
parameters $\gamma_{1}$ and $\gamma_{2}$. The number of true positives (TP),
false positives (FP) and false negatives (FN) are

$\begin{split}\text{TP}(\gamma_{1},\gamma_{2})&=\sum_{i=1}^{3}\sum_{j=1}^{3}\left(\sum_{s<t}a_{st}^{(ij)}\cdot\widehat{a}_{st}^{(ij)}(\gamma_{1},\gamma_{2})\right),\\\
\text{FP}(\gamma_{1},\gamma_{2})&=\sum_{i=1}^{3}\sum_{j=1}^{3}\left(\sum_{s<t}\left(1-a_{st}^{(ij)}\right)\cdot\widehat{a}_{st}^{(ij)}(\gamma_{1},\gamma_{2})\right)\text{
and }\\\
\text{FN}(\gamma_{1},\gamma_{2})&=\sum_{i=1}^{3}\sum_{j=1}^{3}\left(\sum_{s<t}a_{st}^{(ij)}\cdot\left(1-\widehat{a}_{st}^{(ij)}(\gamma_{1},\gamma_{2})\right)\right).\end{split}$

Note that one needs to sum over all nine graphs. Using these definitions, we
can express precision as follows:

$\text{Precision}\left(\gamma_{1},\gamma_{2}\right)=\frac{\text{TP}\left(\gamma_{1},\gamma_{2}\right)}{\text{TP}\left(\gamma_{1},\gamma_{2}\right)+\text{FP}\left(\gamma_{1},\gamma_{2}\right)}.$

And recall as:

$\text{Recall}\left(\gamma_{1},\gamma_{2}\right)=\frac{\text{TP}\left(\gamma_{1},\gamma_{2}\right)}{\text{TP}\left(\gamma_{1},\gamma_{2}\right)+\text{FN}\left(\gamma_{1},\gamma_{2}\right)}.$

The $F_{1}$-score is then simply the harmonic mean of the precision and
recall:

$F_{1}\left(\gamma_{1},\gamma_{2}\right)=2\cdot\frac{\text{Precision}\left(\gamma_{1},\gamma_{2}\right)\cdot\text{Recall}\left(\gamma_{1},\gamma_{2}\right)}{\text{Precision}\left(\gamma_{1},\gamma_{2}\right)+\text{Recall}\left(\gamma_{1},\gamma_{2}\right)}.$

The $F_{1}$-score requires selecting the tuning parameters $\gamma_{1}$ and
$\gamma_{2}$ (or equivalently, $\lambda_{1}$ and $\lambda_{2}$). We consider
all values of $(\gamma_{1},\gamma_{2})\in\Gamma\times\Gamma$. If the selection
is based on the AIC or the BIC, see eq. (20) and (21), respectively, the
$F_{1}$-scores are given by

$F_{1}^{\text{AIC}}=\min_{\left(\gamma_{1},\gamma_{2}\right)\in\Gamma\times\Gamma}\text{AIC}\left(\gamma_{1},\gamma_{2}\right)\quad\text{and}\quad
F_{1}^{\text{BIC}}=\min_{\left(\gamma_{1},\gamma_{2}\right)\in\Gamma\times\Gamma}\text{BIC}\left(\gamma_{1},\gamma_{2}\right).$
(27)

i.e., the scores obtained for the $\gamma_{1}$ and $\gamma_{2}$ values that
minimize the AIC or the BIC. Additionally, we consider the oracle
$F_{1}$-score, which represents the highest $F_{1}$-score achievable using the
CVN algorithm if the true graphs are known. This score is given by

$F_{1}^{\text{oracle}}=\max_{\left(\gamma_{1},\gamma_{2}\right)\in\Gamma\times\Gamma}F_{1}\left(\gamma_{1},\gamma_{2}\right).$
(28)

The oracle score serves as a benchmark to determine the best possible
performance achievable when selecting the optimal tuning parameters.

(a) No smoothing, $\mathbf{W}_{0}$

(b) Grid, $\mathbf{W}_{\text{grid}}$

(c) Fully connected, $\mathbf{W}_{\text{full}}$

Figure 7: A visualization of the three meta-graphs corresponding to the three
weight matrices considered in the simulation study.

##  11 Real Data Analysis

High doses of ionizing radiation (2 Gray [Gy]) used in childhood cancer
treatment are known risk factors for acute myeloid leukemia in children and
for second primary neoplasms later in life. However, the effects of low doses
of radiation ($\leq.05$ Gy), typically used in medical diagnostics, on the
incidence of childhood cancer are less understood.

Our case study uses data from the German ‘Cancer in childhood and molecular-
epidemiology’ study (KiKme), which was collected between 2013 and 2019 from
591 former childhood cancer patients registered in the German Childhood Cancer
Registry and cancer-free controls (Marron et al.,, 2021). Within a sub-sample
of carefully matched 156 participants, a radiation experiment was conducted
where individual reactions of long non-coding ribonucleic acids (lncRNAs) to
different radiation exposures in normal somatic cells of the same person were
investigated. Primary fibroblasts were donated from long-term childhood cancer
survivors who had a first primary cancer during childhood (N1, $n=52$), with
at least one second primary neoplasm (N2+, $n=52$), and tumor-free controls
(N0, $n=52$) (Grandt et al.,, 2022). Cultured fibroblasts of each donor were
exposed to 0 Gy (sham-irradiation), .05 Gy, and 2 Gy X-rays to investigate
cellular reaction to ionizing radiation.

We used $p=191$ genes that are annotated with the GO-term ‘response to
ionizing radiation’ (GO:0010212). All gene expressions were standardized with
mean 0 and a standard deviation of 1 within each class. We estimated all nine
graphs based on the various combinations of the levels of the external
covariates _group_ (N0, N1, N2+) and _radiation_ (0, .05 and 2 Gy).

The tuning parameter to control the sparsity level was set to
$\gamma_{1}=10^{-5}$ to achieve very sparse networks and $\gamma_{2}=10^{-6}$
to enforce smoothness. We assumed a regular grid $\mathbf{W}_{\text{grid}}$,
see eq. (26), and a grid with weights accordingly to the radiation doses:

$\mathbf{W}_{\text{rad}}=\begin{bmatrix}0&.025&0&.5&0&0&0&0&0\\\
.025&0&.975&0&.5&0&0&0&0\\\ 0&.975&0&0&0&.5&0&0&0\\\ .5&0&0&0&.025&0&.5&0&0\\\
0&.5&0&.025&0&.975&0&.5&0\\\ 0&0&.5&0&.975&0&0&0&.5\\\
0&0&0&.5&0&0&0&.025&0\\\ 0&0&0&0&.5&0&.025&0&.975\\\ 0&0&0&0&0&.5&0&.975&0\\\
\end{bmatrix}$

##  12 Results

In this section, we present the results from both the simulation and the case
study.

###  12.1 Simulation Study

We consider 12 different parameter settings and repeat the simulation 20 times
for each setting, see Section 10. We evaluate the $F_{1}$-score using the AIC
and BIC criteria for selecting the tuning parameters $\gamma_{1}$ and
$\gamma_{2}$, see eq. (27). Additionally, we assess the performance in an
‘oracle’ scenario where the true graphs are known, see eq. (28).

Figures 8, 9, and 10 show results from the simulation study. These plots share
a common structure: the top row displays results for the ER graphs, while the
bottom row shows results for the BA graphs. Additionally, the left column
corresponds to the setting with $n=200$ observations, and the right column
corresponds to the setting with $n=100$ observations.

Figure 8 shows the $F_{1}$-scores obtained using three different methods for
selecting the ‘optimal’ tuning parameters: AIC, BIC, and the oracle. These
results pertain to the setting where $\pi_{1}=\pi_{2}=.1$; similar patterns
are observed for other values of $\pi_{1}$ and $\pi_{2}$. The AIC performs
reasonably well when $n=200$, but extremely poorly when $n=100$. The BIC
performs poorly in both cases, since it tends to be too stringent, resulting
in very sparse graphs. However, if the optimal values were known, one could
achieve rather good $F_{1}$-scores. This is shown by the oracle which performs
best since it considers the true underlying structure of the CVN. This at
least shows that, in principle, the CVN can be estimated rather well. Due to
the poor performance of both information criteria, we consider the oracle
$F_{1}$-score for the rest of this section.

Figure 9 shows the best $F_{1}$-score achievable by selecting the ‘optimal’
$\gamma_{1}$ and $\gamma_{2}$ for different weight matrices. Recall that, as
described in Section 9, the grid weight matrix, $\mathbf{W}_{\text{grid}}$,
aligns best with the data-generating process. As expected, using the grid
weight matrix yields the best results when $n=200$. However, with fewer
observations (right column), the fully connected weight matrix,
$\mathbf{W}_{\text{full}}$, performs slightly better. This might be due to the
fact that smoothing all nine graphs simultaneously increases the ‘evidence’
for the existence of edges.

Figure 10 shows the best $F_{1}$-score when the extent to which the graphs
change in the $U_{1}$ and $U_{2}$ directions varies. Whether 10% or 20% of the
expected number of edges change in these directions has little impact on the
algorithm’s performance. This pattern is consistent across other simulation
settings as well.

All the simulation results can be found and explored interactively at
cvn.bips.eu.

Figure 8: The $F_{1}$-scores when the tuning parameters are selected based on
the AIC, BIC, or an oracle. The top row shows the results for the ER graphs,
while the bottom row shows the results for the BA graphs. The left column
corresponds to $n=100$ observations, and the right column corresponds to
$n=200$ observations. The parameters $\pi_{1}=\pi_{2}=.1$ for all plots.
Figure 9: The $F_{1}$-scores when the tuning parameters $\gamma_{1}$ and
$\gamma_{2}$ are selected using an oracle, with the CVN model using either
$\mathbf{W}_{0}$ (no smoothing), $\mathbf{W}_{\text{grid}}$, or
$\mathbf{W}_{\text{full}}$ as the weight matrix. The top row displays results
for the ER graphs, and the bottom row shows results for the BA graphs. The
left column corresponds to $n=100$ observations, and the right column
corresponds to $n=200$ observations. The parameters $\pi_{1}=\pi_{2}=.1$ for
all plots. Figure 10: The $F_{1}$-scores when the tuning parameters
$\gamma_{1}$ and $\gamma_{2}$ are selected using an oracle, with the CVN model
using either $\mathbf{W}_{0}$ (no smoothing), $\mathbf{W}_{\text{grid}}$, or
$\mathbf{W}_{\text{full}}$ as the weight matrix. The top row displays results
for the ER graphs, and the bottom row shows results for the BA graphs. The
number of observations is $n=200$ for all plots. In the left column,
$\pi_{1}=.1$ and $\pi_{2}=.2$, while in the right column, $\pi_{1}=.2$ and
$\pi_{2}=.2$.

###  12.2 Case Study

Figure 11 shows the nine estimated graphs based on the KiKme dataset with
$\gamma_{1}=5\times 10^{-5}$ and $\gamma_{2}=5\times 10^{-6}$ and the weight
matrix $\mathbf{W}_{\text{rad}}$. The nodes, corresponding to the $p=191$ gene
expressions, are positioned on a circle, with lines indicating the edges. Red
lines denote $42$ edges present in _all_ nine graphs, representing the ‘core
graph’ unaffected by external covariates, while blue lines indicate edges
present in only a subset of the graphs, i.e., the edges that are influenced by
the external covariates.

(a) N2+, 0 Gy

(b) N2+, .05 Gy

(c) N2+, 2 Gy

(d) N1, 0 Gy

(e) N1, .05 Gy

(f) N1, 2 Gy

(g) N0, 0 Gy

(h) N0, .05 Gy

(i) N0, 2 Gy

Figure 11: Estimated graphs based on $\mathbf{W}_{\text{rad}}$,
$\gamma_{1}=5\times 10^{-5}$ and $\gamma_{2}=5\times 10^{-6}$. Blue lines
denote $42$ edges present in _all_ nine graphs.

Here, the focus is not just on the graphs themselves, but also on the level of
dis)similarity between them. Here, we consider the Hamming distances between
graphs. Let $\widehat{\mathbf{A}}_{ij}\left(\gamma_{1},\gamma_{2}\right)$ and
$\widehat{\mathbf{A}}_{kl}\left(\gamma_{1},\gamma_{2}\right)$ denote the
estimated adjacency matrices when $(U_{1}=i,U_{2}=j)$ and $(U_{1}=k,U_{2}=l)$,
respectively, and the tuning parameters are $\gamma_{1}$ and $\gamma_{2}$. The
Hamming distance between these two graphs is defined as the number of edges
that is present in one graph but absent in the other:

$\text{Hamming
distance}\left(\gamma_{1},\gamma_{2}\right)=\sum_{s<t}\mathds{1}\left(\widehat{a}_{st}^{(ij)}\left(\gamma_{1},\gamma_{2}\right)\neq\widehat{a}_{st}^{(kl)}\left(\gamma_{1},\gamma_{2}\right)\right).$

This distance reflects the number of edges that must be added and/or removed
to transform one graph into another.

Figure 12 shows a heatmap of the Hamming distances between the nine graphs.
The first three rows/columns represent the N0 control group, the fourth
through sixth rows/columns represent the N1 results, and the last three
rows/columns represent the N2+ results. Within these groups, the first
row/column corresponds to no radiation, the second to .05 Gy, and the third to
the highest dose of 2 Gy. The color gradient represents the Hamming distances.
Within a group, there is almost no difference between no radiation and low
radiation. However, these networks are clearly different from those with the
maximum radiation dose. The networks between cancer groups within an intensity
of irradiation are clearly different, regardless of the radiation dose.

Figure 12: The Hamming distances between the nine graphs of a CVN with $\mathbf{W}_{\text{rad}}$ based on the KiKme data set when $\gamma_{1}=5\times 10^{-5}$ and $\gamma_{2}=5\times 10^{-6}$. The rows/columns are grouped together based on the cancer group the individuals belong to. Group | Gy |  Number of nodes |  Number of edges |  Graph density |  Avg. shortest path length |  Avg. Degree |  Max. Degree and nodes
---|---|---|---|---|---|---|---
N0 | 0 | 31 | 75 | .16 | 2.38 | 4.84 | 21
N0 | .05 | 31 | 76 | .16 | 2.33 | 4.9 | 21
N0 | 2 | 27 | 68 | .19 | 2.15 | 5.04 | 21
N1 | 0 | 31 | 63 | .14 | 2.33 | 4.06 | 15
N1 | .05 | 30 | 61 | .14 | 2.34 | 4.07 | 14
N1 | 2 | 25 | 57 | .19 | 2.24 | 4.56 | 13
N2+ | 0 | 32 | 66 | .13 | 2.69 | 4.12 | 16
N2+ | .05 | 32 | 67 | .14 | 2.65 | 4.19 | 16
N2+ | 2 | 29 | 65 | .16 | 2.16 | 4.48 | 19
Table 2: Characteristics of the 9 estimated subgraphs. Nodes in each subgraph
have at least one adjacent node. Singletons were removed from the original
graph with $p=191$ nodes.

Table 2 demonstrates that within each CVN subgraph, only a subset of nodes are
interconnected. There are only few differences within a group between
radiation doses of 0 and .05 Gy, which may be due to random variation.
Notably, the size of the subsets decreases with increasing radiation dose,
whereas the density and the average node degree in each subgraph increases,
implying a higher level of inter-connectivity among the remaining genes. The
TP53 gene consistently exhibits the highest degree centrality, indicating its
importance in the network regardless of the cancer group or radiation dosage.
The TP53 gene encodes the tumor suppressor protein p53 that regulates DNA
repair, apoptosis, and cell cycle control. It plays a critical role in the
prevention of cancer, and mutations in this gene are common in human tumors.

Gene | Chr | Description | Found in |  N00 |  N01 |  N02 |  N10 |  N11 |  N12 |  N20 |  N21 |  N22
---|---|---|---|---|---|---|---|---|---|---|---|---
ATG4A | X | Autophagy Related 4a Cysteine Peptidase | N2+ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1
RAD9B | 12 | Rad9 Checkpoint Clamp Component B | N2 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0
STX6 | 1 | Syntaxin 6 | N2+ | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1
TRIM13 | 13 | Tripartite Motif Containing 13 | N2+ | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1
XRCC4 | 5 | X-Ray Repair Cross Complementing 4 | N2+ | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1
CES2 | 16 | Carboxylesterase 2 | 2 Gy | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1
SESN2 | 1 | Sestrin 2 | 2 Gy | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1
FBXO22 | 15 | F-Box Protein 22 | in cancer group(s) and for N0 with 2 Gy | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1
XRCC6 | 22 | X-Ray Repair Cross Complementing 6 | in cancer group(s) and for N0 with 2 Gy | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1
GRB2 | 17 | Growth Factor Receptor Bound Protein 2 | control group | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0
PRKAA1 | 5 | Protein Kinase Amp-Activated Catalytic Subunit Alpha 1 | control group | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0
TAF3 | 10 | Tata-Box Binding Protein Associated Factor 3 |  in N0 and N1 with $<$= .05 Gy | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0
NOX4 | 11 | Nadph Oxidase 4 | control group | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0
BAX | 19 | Bcl2 Associated X, Apoptosis Regulator | control group | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0
ASCC3 | 6 | Activating Signal Cointegrator 1 Complex Subunit 3 | control group | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0
TP53INP1 | 8 | Tumor Protein P53 Inducible Nuclear Protein 1 | control group | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0
SPIDR | 8 | Scaffold Protein Involved In Dna Repair | 0-2 Gy in control group, unclear pattern for cancer groups | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1
DNMT3A | 2 | Dna Methyltransferase 3 Alpha | 0-2 Gy in control group, unclear pattern for cancer groups | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1
NABP1 | 2 | Nucleic Acid Binding Protein 1 | 0-2 Gy in control group, unclear pattern for cancer groups | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0
GADD45A | 1 | Growth Arrest And Dna Damage Inducible Alpha | 0-2 Gy in control group, unclear pattern for cancer groups | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0
RHNO1 | 12 | Rad9-Hus1-Rad1 Interacting Nuclear Orphan 1 | 0-2 Gy in control group, unclear pattern for cancer groups | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0
RAD9A | 11 | Rad9 Checkpoint Clamp Component A | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
HMGA2 | 12 | High Mobility Group At-Hook 2 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
BLM | 15 | Blm Recq Like Helicase | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
XRCC5 | 2 | X-Ray Repair Cross Complementing 5 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
CASP3 | 4 | Caspase 3 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
ANKRA2 | 5 | Ankyrin Repeat Family A Member 2 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
SESN1 | 6 | Sestrin 1 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
TNFRSF10B | 8 | Tnf Receptor Superfamily Member 10b | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
ACER2 | 9 | Alkaline Ceramidase 2 | in all graphs | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1

Table 3: Genes ajdacent to TP53 in the estimated CVN.

Table 3 shows all adjacent genes of TP53 in each subgraph, nine of them being
shared by all of them. Interestingly, TP53 is only associated with BAX in the
control group for all radiation doses. The expression of BAX is upregulated in
the presence of the tumor suppressor p53 following significant cellular
damage, leading to the initiation of apoptotic (programmed) cell death. This
is beneficial for cancer prevention. If there is a mutation in either TP53 or
BAX, the process of TP53 initiating BAX may be disrupted and contribute to
cancer development.

##  13 Conclusions & Discussion

We introduced a new Gaussian graphical model called the covariate-varying
network (CVN). This model allows the graph structure to change with multiple
discrete external covariates, e.g., time and groups. To achieve this, we
applied a LASSO penalty to the entries of the precision matrix to introduce
sparsity. Additionally, we enforce similarities, or ‘smoothness’, between the
different graphs by penalizing the absolute differences between them. We
modeled smoothness using the concept of a meta graph with weighted edges, see
Section 3. This approach is the first to address the problem of a graph
changing with multiple external covariates. Many previously existing graphical
models are special cases within the CVN graphical model class, see Section
3.1.

We solved the resulting maximum likelihood optimization problem associated
with the CVN model, see eq. (5), using an Alternating Direction Method of
Multipliers (ADMM), see Section 4. In the second update step of the algorithm
(Section 4.2), we discovered that the problem could be subdivided into
individual optimization problems, one for each edge. This subproblem is
equivalent to a weighted Fused LASSO Signal Approximator (wFLSA), for which we
developed a specialized algorithm previously, see Dijkstra et al., (2024).

We determined that the complexity of our algorithm is cubic with respect to
the number of variables ($p^{3}$) and squared with respect to the number of
graphs ($m^{2}$), see Section 5. This is a downside, as the algorithm becomes
computationally expensive and challenging to run for very large instances. The
primary bottleneck is the second update step, where the wFLSA must be solved
for each of the $\binom{p}{2}$ edges.

We proposed a different type of tuning parameter parameterization as is common
for graphical models, see Section 8. Instead of using tuning parameters that
penalize the $L_{1}$-norm of the precision matrix ($\lambda_{1}$) and the
absolute differences between precision matrices ($\lambda_{2}$), our proposed
parameters penalize individual entries of a precision matrix ($\gamma_{1}$)
and individual differences between individual edges ($\gamma_{2}$). The
advantage of this approach is that once appropriate values for $\gamma_{1}$
and $\gamma_{2}$ are found, they are likely to change only slightly when the
dataset is extended with additional variables or graphs, i.e., when $p$ and
$m$ change. This stability would not be the case for $\lambda_{1}$ and
$\lambda_{2}$. Additionally, this parameterization is not only useful for the
CVN model but could also benefit other LASSO-type problems.

We conducted a simulation study (see Sections 9 and 10) using two types of
graphs: Erdős–Rényi and Barabasi–Albert. These graphs were influenced by two
external covariates, resulting in varied changes to the graph structure based
on each covariate. We found that the estimation of the CVN was generally
successful across different scenarios. However, selecting the tuning
parameters poses significant challenges. The extent to which the graphs change
with the external covariates did not appear to impact performance. The authors
anticipate that if higher percentages were considered, one would see a decline
in performance. It is important to note that comparing this method with others
is difficult since, to the best of our knowledge, there are currently no other
methods capable of handling the CVN graphical model.

To assess its applicability, we applied our method to a real dataset from the
KiKme study (see Section 11). We showed that the method is able to recover
differences and similarities of gene expression network structures that are
specific to three different case/control groups and three different radiation
dosages. This case study demonstrated that CVN is capable of estimating
plausible networks in a high-dimensional setting where the overall number of
samples ($n=156$) is smaller than the number of nodes ($p=191$). The number of
samples for each subgraph was therefore much smaller.

Moving forward, there are several avenues for future research. First,
extending the method to handle continuous covariates would be a promising
direction, broadening its applicability across different datasets. Currently,
dealing with continuous external covariates necessitates data discretization.
Second, based on our simulation study, we observed that the performance of AIC
and BIC for tuning parameter selection is notably inadequate. Exploring
alternative methods for parameter tuning could potentially improve
performance. Lastly, due to the high algorithmic complexity of our method,
exploring alternative optimization strategies, including heuristic approaches,
could help reduce computational costs and enhance efficiency.

The algorithm has been implemented as an R package and is publicly accessible
at www.github.com/bips-hb/CVN. Additionally, the CVN data simulator is
available as an R package at www.github.com/bips-hb/CVNSim, and the code for
the simulation study can be found at www.github.com/bips-hb/CVNStudy. All
results of the simulation study are hosted at cvn.bips.eu.

## Acknowledgments

We gratefully acknowledge the financial support of the German Research
Foundation (DFG – Project FO 1045/2-1).

## References

  * Ando, (2007) Ando, T. (2007).  Bayesian predictive information criterion for the evaluation of hierarchical bayesian and empirical bayes models.  Biometrika, 94(2):443–458.
  * Badhwar et al., (2017) Badhwar, A., Tam, A., Dansereau, C., Orban, P., Hoffstaedter, F., and Bellec, P. (2017).  Resting-state network dysfunction in Alzheimer’s disease: A systematic review and meta-analysis.  Alzheimer’s & Dementia: Diagnosis, Assessment & Disease Monitoring, 8:73–85.
  * Banjeree and Ghaoui, (2008) Banjeree, O. and Ghaoui, L. E. (2008).  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  Journal of Machine Learning Research, 9:458–516.
  * Barabasi and Oltvai, (2004) Barabasi, A.-L. and Oltvai, Z. N. (2004).  Network biology: Understanding the cell’s functional organization.  Nature Reviews Genetics, 5(2):101–113.
  * Boyd et al., (2010) Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2010).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends in Machine Learning, 3:1–122.
  * Boyd and Vandenberghe, (2004) Boyd, S. and Vandenberghe, L. (2004).  Convex optimization.  Cambridge University Press.
  * Brent, (2013) Brent, R. P. (2013).  Algorithms for minimization without derivatives.  Courier Corporation.
  * Brozyna et al., (2007) Brozyna, A., Zbytek, B., Granese, J., Carlson, J. A., Ross, J., and Slominski, A. (2007).  Mechanism of UV-related carcinogenesis and its contribution to nevi/melanoma.  Expert Review of Dermatology, 2(4):451–469.
  * Cai et al., (2011) Cai, T., Liu, W., and Luo, X. (2011).  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association, 106:594–607.
  * Chen and Chen, (2008) Chen, J. and Chen, Z. (2008).  Extended bayesian information criteria for model selection with large model spaces.  Biometrika, 95(3):759–771.
  * Curtis et al., (2012) Curtis, H., Blaser, M. J., Dirk, G., Kota, K. C., Rob, K., Liu, B., Wang, L., Sahar, A., White, J. R., Badger, J. H., et al. (2012).  Structure, function and diversity of the healthy human microbiome.  Nature, 486(7402):207–214.
  * Danaher et al., (2014) Danaher, P., Wang, P., and Witten, D. M. (2014).  The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society. Series B: Statistical Methodology, 76:373–397.
  * de la Cuesta-Zuluaga et al., (2018) de la Cuesta-Zuluaga, J., Corrales-Agudelo, V., Velásquez-Mejía, E. P., Carmona, J. A., Abad, J. M., and Escobar, J. S. (2018).  Gut microbiota is associated with obesity and cardiometabolic disease in a population in the midst of westernization.  Scientific Reports, 8(1):11356.
  * de Vos and de Vos, (2012) de Vos, W. M. and de Vos, E. A. (2012).  Role of the intestinal microbiome in health and disease: From correlation to causation.  Nutrition Reviews, 70:S45–S56.
  * Dijkstra et al., (2024) Dijkstra, L., Hanke, M., Koenen, N., and Foraita, R. (2024).  An alternating direction method of multipliers algorithm for the weighted fused lasso signal approximator.  Preprint.  arXiv:2407.18077.
  * Friedman et al., (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9:432–441.
  * Gibberd and Nelson, (2017) Gibberd, A. J. and Nelson, J. D. (2017).  Regularized estimation of piecewise constant Gaussian graphical models: The group-fused graphical Lasso.  Journal of Computational and Graphical Statistics, 26(3):623–634.
  * Grandt et al., (2022) Grandt, C. L., Brackmann, L. K., Poplawski, A., Schwarz, H., Hummel-Bartenschlager, W., Hankeln, T., Kraemer, C., Marini, F., Zahnreich, S., Schmitt, I., et al. (2022).  Radiation-response in primary fibroblasts of long-term survivors of childhood cancer with and without second primary neoplasms: The KiKme study.  Molecular Medicine, 28(1):105.
  * Hallac et al., (2017) Hallac, D., Park, Y., Boyd, S., and Leskovec, J. (2017).  Network inference via the time-varying graphical lasso.  Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Part F1296:205–213.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H. (2009).  The elements of statistical learning: Data mining, inference, and prediction, volume 2.  Springer.
  * Lauritzen, (1996) Lauritzen, S. (1996).  Graphical models.  Oxford Statistical Science Series. Clarendon Press.
  * Liu and Wang, (2017) Liu, H. and Wang, L. (2017).  TIGER: A tuning-insensitive approach for optimally estimating Gaussian graphical models.  Electronic Journal of Statistics, 11:241–294.
  * Marron et al., (2021) Marron, M., Brackmann, L. K., Schwarz, H., Hummel-Bartenschlager, W., Zahnreich, S., Galetzka, D., Schmitt, I., Grad, C., Drees, P., Hopf, J., et al. (2021).  Identification of genetic predispositions related to ionizing radiation in primary human skin fibroblasts from survivors of childhood and second primary cancer as well as cancer-free controls: Protocol for the nested case-control study kikme.  JMIR Research Protocols, 10(11):e32395.
  * Monti et al., (2014) Monti, R. P., Hellyer, P., Sharp, D., Leech, R., Anagnostopoulos, C., and Montana, G. (2014).  Estimating time-varying brain connectivity networks from functional MRI time series.  Neuroimage, 103:427–443.
  * Newman, (2018) Newman, M. (2018).  Networks.  Oxford University Press.
  * Pósfai and Barabasi, (2016) Pósfai, M. and Barabasi, A.-L. (2016).  Network science.  Citeseer.
  * Uhler, (2018) Uhler, C. (2018).  Gaussian graphical models.  In Handbook of graphical models, pages 217–238. CRC Press.
  * Vinciotti et al., (2024) Vinciotti, V., Wit, E. C., and Richter, F. (2024).  Random graphical model of microbiome interactions in related environments.  Journal of Agricultural, Biological and Environmental Statistics, pages 1–14.
  * Witten and Tibshirani, (2009) Witten, D. M. and Tibshirani, R. (2009).  Covariance-regularized regression and classification for high dimensional problems.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 71(3):615–636.
  * Wu et al., (2019) Wu, N., Huang, J., Zhang, X.-F., Ou-Yang, L., He, S., Zhu, Z., and Xie, W. (2019).  Weighted fused pathway graphical lasso for joint estimation of multiple gene networks.  Frontiers in Genetics, 10:1–12.

##  Appendix A Interpolation

As outlined in Section 6, interpolating a graph $G_{m+1}$ using an estimated
CVN model $\widehat{\mathbf{\Theta}}_{\text{CVN}}$ involves solving an
optimization problem for each possible edge $\\{s,t\\}$ in the graph:

$\widehat{\theta}_{st}^{(m+1)}=\operatorname*{arg\,min}_{\theta_{st}^{(m+1)}\in\mathds{R}}\lambda_{1}|\theta_{st}^{(m+1)}|+\lambda_{2}\lVert\bm{\omega}^{\top}\widehat{\bm{\theta}}_{st}-\theta_{st}^{(m+1)}\rVert_{1}.$

Here, $\widehat{\theta}_{st}^{(m+1)}$ represents the interpolated entry for
the edge $\\{s,t\\}$ in the precision matrix. The vector
$\bm{\omega}=\left(\omega_{1},\omega_{2},\ldots,\omega_{m}\right)^{\top}$
contains predetermined, non-negative weights, and
$\widehat{\bm{\theta}}_{st}=\left(\widehat{\theta}_{st}^{(1)},\widehat{\theta}_{st}^{(2)},\ldots,\widehat{\theta}_{st}^{(m)}\right)^{\top}$
are the estimated values in the precision matrices for edge $\\{s,t\\}$ of the
$m$ graphs. The tuning parameters are $\lambda_{1}>0$ and $\lambda_{2}\geq 0$.

Although this optimization problem lacks an analytical solution, it is convex,
allowing for solving it numerically using derivative-free search methods. We
employ Brent’s algorithm for this purpose (Brent,, 2013). For ease of
notation, let $x=\theta_{st}^{(m-1)}$ and $\bm{y}=\widehat{\bm{\theta}}_{st}$.
The goal is to solve:

$\operatorname*{arg\,min}_{x\in\mathds{R}}f(x;\bm{y},\bm{\omega})=\operatorname*{arg\,min}_{x\in\mathds{R}}\lambda_{1}|x|+\lambda_{2}\left\lVert\bm{\omega}^{\top}\bm{y}-x\right\rVert_{1}.$

The convexity of the function $f(x;\bm{y},\bm{\omega})$ in $x$ is evident
since

$\begin{split}f(\alpha x_{1}+(1-\alpha)x_{2};\bm{y},\bm{\omega})&\leq\alpha
f(x_{1};\bm{y},\bm{\omega})+(1-\alpha)f(x_{2};\bm{y},\bm{\omega})\end{split}$

for all $x_{1},x_{2}\in\mathds{R}$, $\alpha\in[0,1]$, $\lambda_{1}>0$,
$\lambda_{2}\geq 0$, $\bm{\omega}\in\mathds{R}_{+}^{m}$ and
$\bm{y}\in\mathds{R}^{m}$. This is easy to see for both

$|\alpha x_{1}+(1-\alpha)x_{2}|\leq\alpha|x_{1}|+(1-\alpha)|x_{2}|$

and

$\left\lVert\bm{\omega}^{\top}\bm{y}-(\alpha
x_{1}+(1-\alpha)x_{2})\right\rVert_{1}\leq\alpha\left\lVert\bm{\omega}^{\top}\bm{y}-x_{1}\right\rVert_{1}+(1-\alpha)\left\lVert\bm{\omega}^{\top}\bm{y}-x_{2}\right\rVert_{1}.$

Note that the sign of $x$ (i.e., $\widehat{\theta}_{st}^{(m+1)}$) is
irrelevant in this context, as our objective is solely to determine the
existence of the edge $\\{s,t\\}$ in the interpolated graph, i.e., whether
$\widehat{\theta}_{st}^{(m+1)}\neq 0$.