Covariate-Adjusted Inference for Differential Analysis of High-Dimensional Networks

To formalize our problem, we begin by introducing some notation. We compare networks between two groups, labeled by $g\in\{\mathrm{I},\mathrm{II}\}$. We obtain measurements of $p$ variables $X^{g}=\left(X^{g}_{1},\ldots,X_{p}^{g}\right)^{\top}$, corresponding to nodes in a graphical model [26], on $n^{\mathrm{I}}$ subjects in group $\mathrm{I}$ and $n^{\mathrm{II}}$ subjects in group $\mathrm{II}$. We define $\mathcal{X}\subseteq\mathds{R}^{p}$ as the sample space of $X^{g}$. Let $X^{g}_{i,j}$ denote the data for node $j$ for subject $i$ in group $g$, and let $\mathbf{X}_{j}^{g}=(X^{g}_{1,j},\ldots,X^{g}_{n^{g},j})^{\top}$ be an $n^{g}$-dimensional vector of measurements on node $j$ for group $g$.

Our objective is to determine whether the association between variables $X_{j}$ and $X_{k}$, conditional upon all other variables, differs by group. Our approach is to specify a model for $X^{g}$ such that the conditional dependence between any two nodes $X^{g}_{j}$ and $X^{g}_{k}$ can be represented by a single scalar parameter $\beta^{g,*}_{j,k}$. If the association between nodes $j$ and $k$ is the same in both groups $\mathrm{I}$ and $\mathrm{II}$, $\beta^{\mathrm{I},*}_{j,k}=\beta^{\mathrm{II},*}_{j,k}$. Conversely, if $\beta^{\mathrm{I},*}_{j,k}\neq\beta^{\mathrm{II},*}_{j,k}$, we say nodes $j$ and $k$ are differentially connected. We assess for differential connectivity by performing a test of the null hypothesis

$\displaystyle H^{0}_{j,k}:\beta^{\mathrm{I},*}_{j,k}=\beta^{\mathrm{II},*}_{j,k}.$

(1)

We consider a general class of exponential family pairwise interaction models. For $x=(x_{1},\ldots,x_{p})^{\top}$, we assume the density function for $X^{g}$ takes the form

$\displaystyle f^{g,*}(x)=\exp\left(\sum_{j=1}^{p}\mu_{j}(x_{j})+\sum_{j=1}^{p}\sum_{k=1}^{j}\beta^{g,*}_{j,k}\psi_{j,k}(x_{j},x_{k})-U\left(\boldsymbol{\beta}^{g,*}\right)\right),$

(2)

where $\psi_{j,k}$ and $\mu_{j}$ are fixed and known functions, $\boldsymbol{\beta}^{g,*}$ is a $p\times p$ matrix with elements $\beta^{g,*}_{j,k}$, and $U(\boldsymbol{\beta}^{g,*})$ is the log-partition function. The dependence between $X^{g}_{j}$ and $X^{g}_{k}$ is measured by $\beta^{g,*}_{j,k}$, and nodes $j$ and $k$ are conditionally independent in group $g$ if and only if $\beta^{g,*}_{j,k}=0$.

This class of exponential family distributions is rich and includes several models that have been studied previously in the graphical modeling literature. One such example is the Gaussian graphical model, perhaps the most widely-used graphical model for continuous data. For $x\in\mathds{R}^{p}$ the density function for mean-centered Gaussian random vectors can be expressed as

$\displaystyle f^{g,*}(x)\propto\exp\left(-\sum_{j=1}^{p}\sum_{k=1}^{j}\beta^{g,*}_{j,k}x_{j}x_{k}\right),$

(3)

and is thus a special case of (2) with $\psi_{j,k}=-x_{j}x_{k}$ and $\mu_{j}=0$. The non-negative Gaussian density, which takes the form of (3) with the constraint that $x$ takes values in $\mathds{R}^{p}_{+}$, also belongs to the exponential family class. Another canonical example is the Ising model, commonly used for studying conditional dependencies among binary random variables. For $x\in\{0,1\}^{p}$, the density function for the Ising model can be expressed as

$\displaystyle f(x)\propto\exp\left(\sum_{j=1}^{p}\sum_{k=1}^{j}\beta_{j,k}x_{j}x_{k}\right).$

Additional examples include the Poisson model, the exponential graphical model, and conditionally-specified mixed graphical models [40, 7].

When asymptotically normal estimates of $\beta^{\mathrm{I},*}_{j,k}$ and $\beta^{\mathrm{II},*}_{j,k}$ are available, one can perform a calibrated test of $H^{0}_{j,k}$ based on the difference between the estimates. In many cases, asymptotically normal estimates can be obtained using well-established methodology. For instance, when the log-partition function $U(\boldsymbol{\beta}^{g,*})$ is available in closed form and is tractable, one can obtain estimates via (penalized) maximum likelihood. This is a standard approach in the Gaussian setting, in which case the log-partition function is easy to compute. However, this is not the case for other exponential family models. Likelihood-based estimation strategies are thus generally difficult to implement. In this paper, we consider two alternative strategies that have been proposed to overcome these computational challenges and are more broadly applicable.

The first approach we discuss is neighborhood selection [7, 27, 40]. Consider a sub-class of exponential family graphical models for which the conditional density function for any node $X^{g}_{j}$ given the remaining nodes belongs to a univariate exponential family model. Because the log-partition function in univariate exponential family models is available in closed form, it is computationally feasible to estimate each conditional density function. By estimating the conditional density functions, one can identify the neighbors of nodes $j$, that is, the nodes upon which the conditional distribution depends. This approach was first proposed as an alternative to maximum likelihood estimation for estimating Gaussian graphical models [27]. To describe our approach, we focus on the Gaussian case, though this approach is more widely applicable and can be used for modeling dependencies among, e.g., Poisson, binomial, and exponential random variables as well [7, 40].

In Gaussian graphical models, the dependency of node $j$ on all other nodes can be determined based on the linear model

$\displaystyle\mathds{E}\left[X^{g}_{j}|X^{g}_{1},\ldots,X^{g}_{p}\right]=\beta^{g,*}_{j,0}+\sum_{k\neq j}\beta^{g,*}_{j,k}X^{g}_{k}.$

(4)

The regression coefficients $\beta^{g,*}_{j,k}$ measure the strength of linear association between nodes $j$ and $k$ conditional upon all other nodes and are zero if and only if nodes $j$ and $k$ are conditionally independent; $\beta^{g,*}_{j,0}$ is an intercept term and is zero if all nodes are mean-centered. (We acknowledge a slight abuse of notation here, as the regression coefficients in (4) are not equivalent to parameters in (2). However, either estimand fully characterizes conditional independence.) In the low-dimensional setting (i.e., $p\ll n^{g}$), statistically efficient and asymptotically normal estimates of the regression coefficients can be readily obtained via ordinary least squares. In high-dimensions (i.e., $p\geq n^{g}$), the ordinary least squares estimates are inconsistent, so to obtain consistent estimates we typically rely upon regularized estimators such as the LASSO and the elastic net [33, 51]. Regularized estimators are generally biased and have intractable sampling distributions, and as such, are unsuitable for performing formal statistical inference. However, several methods have recently emerged for obtaining asymptotically normal estimates by correcting the bias of regularized estimators [20, 12, 47].

The second computationally efficient approach we consider is to estimate the density function using the score matching framework of Hyvärinen [17, 18]. Hyvärinen derives a loss function for estimation of density functions for continuous random variables that is based on the gradient of the log-density with respect to the observations. As such, the score matching loss does not depend on the log-partition function in exponential family models. Moreover, when the joint distribution for $X^{g}$ belongs to an exponential family model, the loss is quadratic in the unknown parameters, allowing for efficient computation. In low dimensions, the minimizer of the score matching loss is consistent and asymptotically normal. In high dimensions, one can obtain asymptotically normal estimates by minimizing a regularized version of the score matching loss to obtain an initial estimate [23, 44] and subsequently correcting for the bias induced by regularization [43].

We now consider the setting in which the within-group networks depend on covariates. We denote by $W^{g}$ a $q$-dimensional random vector of covariate measurements for group $g$, and we define $\mathcal{W}$ as the sample space of $W^{g}$. Let $W^{g}_{i,r}$ refer to the value of covariate $r$ for subject $i$ in group $g$, and let $W^{g}_{i}=(W^{g}_{i,1},\ldots,W^{g}_{i,q})^{\top}$ be a $q$-dimensional vector containing all covariates for subject $i$ in group $g$. We assume the number of covariates is small relative to the sample size (i.e., $q\ll n^{g}$).

To study the dependence of the within-group networks on the covariates, we specify a model for the nodes $X^{g}$ given the covariates $W^{g}$ that allows for the inter-node dependencies to vary as a function of $W^{g}$. The model defines a function $\eta^{g}_{j,k}:\mathcal{W}\to\mathds{R}$ that takes as input a vector of covariates and returns a measure of association between nodes $j$ and $k$ for a subject in group $g$ with identical covariates. One can interpret $\eta^{g,*}_{j,k}$ as a conditional version of $\beta^{g,*}_{j,k}$, given the covariates.

We assume that $\eta^{g,*}_{j,k}$ can be written as a low-dimensional linear basis expansion in $W^{g}$ of dimension $d$ — that is,

$\displaystyle\eta^{g,*}_{j,k}\left(W^{g}\right)=\left\langle\phi\left(W^{g}\right),\alpha_{j,k}^{g,*}\right\rangle,$

(5)

where $\phi:\mathds{R}^{q}\to\mathds{R}^{d}$ is a map from a set of covariates to its expansion, $\alpha_{j,k}^{g,*}$ is a $d$-dimensional vector, and $\langle\cdot,\cdot\rangle$ denotes the vector inner product. Let $\phi_{c}(w)$ refer to the $c$-th element of $\phi(w)$. One can take the simple approach of specifying $\phi$ as a linear basis, $\phi(w)=\left(1,w_{1},\ldots,w_{q}\right)$ for $w\in\mathds{R}^{q}$, though more flexible choices such as polynomial or B-spline bases can also be considered. It may be preferable to specify $\phi$ so that $\eta^{g,*}_{j,k}$ is an additive function of the covariates. This allows one to easily assess the effect of any specific covariate on the network by estimating the sub-vector of $\alpha^{g,*}_{j,k}$ that is relevant to the covariate of interest.

When the association between nodes $j$ and $k$ does not depend on group membership, $\eta^{\mathrm{I,*}}_{j,k}(w)=\eta^{\mathrm{II,*}}_{j,k}(w)$ for all $w$, and $\alpha^{\mathrm{I,*}}_{j,k}=\alpha^{\mathrm{II,*}}_{j,k}$. In other words, if one subject from group $\mathrm{I}$ and another subject from group $\mathrm{II}$ have identically-valued covariates, the corresponding measure of association between nodes $j$ and $k$ is also the same. In the covariate-adjusted setting, we say that nodes $j$ and $k$ are differentially connected if there exists $w$ such that $\eta^{\mathrm{I,*}}_{j,k}(w)\neq\eta^{\mathrm{II,*}}_{j,k}(w)$, or equivalently, if $\alpha^{\mathrm{I,*}}_{j,k}\neq\alpha^{\mathrm{II,*}}_{j,k}$. We can thus assess differential connectivity between nodes $j$ and $k$ by testing the null hypothesis

$\displaystyle G^{0}_{j,k}:\alpha^{\mathrm{I,*}}_{j,k}=\alpha^{\mathrm{II,*}}_{j,k}.$

(6)

Similar to the unadjusted setting, when asymptotically normal estimates of $\alpha^{\mathrm{I},*}_{j,k}$ and $\alpha^{\mathrm{II},*}_{j,k}$ are available, a calibrated test can be constructed based on the difference between the estimates.

We now specify a form for the conditional distribution of $X^{g}$ given $W^{g}$ as a generalization of the exponential family pairwise interaction model (2). We assume the conditional density for $X^{g}$ given $W^{g}$ can be expressed as

$\displaystyle f^{g,*}(x|w)\propto\exp\left(\sum_{j=1}^{p}\mu_{j}(x_{j})+\sum_{j=1}^{p}\sum_{k=1}^{j}\eta^{g,*}_{j,k}(w)\psi_{j,k}(x_{j},x_{k})+\sum_{j=1}^{p}\sum_{c=1}^{d}\theta_{j,c}^{g,*}\zeta_{j,c}\left(x_{j},\phi_{c}(w)\right)\right),$

(7)

where $w=(w_{1},\ldots,w_{q})^{\top}$, and the proportionality is up to a normalizing constant that does not depend on $x$. Above, $\zeta_{j,c}$ is a fixed and known function, and the main effects of the covariates on $X^{g}$ are represented by the scalar parameters $\theta^{g,*}_{j,c}$. The conditional dependence between nodes $j$ and $k$, given all other nodes and given that $W^{g}=w$ is quantified by $\eta^{g,*}_{j,k}(w)$, and $\eta^{g,*}_{j,k}(w)=0$ if and only if nodes $j$ and $k$ are conditionally independent at $w$. One can thus view $\eta^{g,*}_{j,k}$ as a conditional version of $\beta^{g,*}_{j,k}$ in (7).

Either of the estimation strategies introduced in Section 2.1 can be used to perform covariate-adjusted inference. When the conditional distribution of each node given the remaining nodes and the covariates belongs to a univariate exponential family model, the covariate-dependent network can be estimated using neighborhood selection because the node conditional distributions can be estimated efficiently with likelihood-based methods. Alternatively, we can estimate the conditional density function (7) using score matching.

As a working example, we again consider estimation of covariate-dependent Gaussian networks using neighborhood selection. Suppose the conditional distribution of $X^{g}$ given $W^{g}$ takes the form

$\displaystyle f^{g,*}(x|w)\propto\exp\left(-\sum_{j=1}^{p}\sum_{k=1}^{j}\eta^{g,*}_{j,k}(w)x_{j}x_{k}-\sum_{j=1}^{p}\sum_{c=1}^{d}\theta_{j,c}^{g,*}x_{j}\phi_{c}(w)\right).$

(8)

Then the dependencies of node $j$ on all other nodes can be determined based on the following varying coefficient model [14]:

$\displaystyle\mathds{E}\left[X^{g}_{j}|X_{1}^{g},\ldots,X_{p}^{g},W^{g}\right]=\eta_{j,0}^{g,*}(W^{g})+\sum_{k\neq j}\eta^{g,*}_{j,k}\left(W^{g}\right)X_{k}^{g}.$

(9)

The varying coefficient model is a generalization of the linear model that treats the regression coefficients as functions of the covariates. In (9), $\eta^{g,*}_{j,k}(w)$ returns a regression coefficient that quantifies the linear relationship between nodes $j$ and $k$ for subjects in group $g$ with covariates equal to $w$. Then $X^{g}_{j}$ and $X^{g}_{k}$ are conditionally independent given all other nodes and given $W^{g}=w$ if and only if $\eta^{g,*}_{j,k}(w)=0$. The varying coefficients $\eta^{g,*}_{j,k}$ can thus be viewed as a conditional version of the regression coefficients in (4). (We have again abused the notation, as the varying coefficient functions in (9) are not equal to the parameters in (8), though both functions are zero for the same values of $w$). The intercept term $\eta^{g,*}_{j,0}$ accounts for the main effect of $W^{g}$ on $X^{g}_{j}$. We can remove this main effect term by first centering the nodes $X^{g}_{j}$ about their conditional mean given $W^{g}$ (which can be estimated by performing a linear regression of $X^{g}_{j}$ on $\phi(W^{g})$).

In Sections 3 and 4, we discuss construction of asymptotically normal estimators of $\alpha^{g,*}_{j,k}$ in the low- and high-dimensional settings using neighborhood selection and score matching. Before proceeding, we first examine the connection between the null hypotheses $H^{0}_{j,k}$ and $G^{0}_{j,k}$.

Hypotheses $H^{0}_{j,k}$ in (1) and $G^{0}_{j,k}$ in (6) are related but not equivalent. It is possible that $H^{0}_{j,k}$ holds while $G^{0}_{j,k}$ fails and vice versa. We provide an example below. Suppose we are using neighborhood selection to perform differential network analysis in the Gaussian setting, so we are making a comparison of linear regression coefficients between the two groups. Suppose further that the within-group networks depend on single scalar covariate $W^{g}$, and the nodes are centered about their conditional mean given $W^{g}$. One can show that the regression coefficients $\beta^{g,*}_{j,k}$ are equal to the average of their conditional versions $\eta^{g,*}_{j,k}(W^{g})$. That is, $\beta^{g,*}_{j,k}=\mathds{E}[\eta^{g,*}_{j,k}(W^{g})]$. Now, suppose $G^{0}_{j,k}$ holds. If $W^{\mathrm{I}}$ and $W^{\mathrm{II}}$ do not share the same distribution (e.g., the covariate tends to take higher values in group $\mathrm{I}$ than in group $\mathrm{II}$), the average conditional inter-node association may differ, and $H^{0}_{j,k}$ may not hold. Although the conditional association between nodes, given the covariate, does not differ by group, the average conditional association does differ, as illustrated in Figure 1a. In such a scenario, the difference in the average conditional association is induced by the dependence of the covariate on group membership and the dependence of the inter-node association on the covariate. Thus, inequality of $\beta^{\mathrm{I},*}_{j,k}$ and $\beta^{\mathrm{II},*}_{j,k}$ does not necessarily capture a meaningful association between the network and group membership. Similarly when $H^{0}_{j,k}$ holds, it is possible that $\eta^{\mathrm{I},*}_{j,k}\neq\eta^{\mathrm{II},*}_{j,k}$. For instance, suppose that the distribution of the covariate is the same in both groups, and $\mathds{E}[\eta^{g}(W^{g})]=0$ in both groups. If the between-node association depends more strongly upon the covariates in one group than the other, $G^{0}_{j,k}$ will be false. This example is depicted in Figure 1b. In this scenario, adjusting for covariates should provide improved power to detect differential connections. We note that for other distributions, it does not necessarily hold that $\beta^{g,*}_{j,k}=\mathds{E}[\eta^{g,*}_{j,k}(W^{g})]$, but regardless, there is generally no equivalence between hypotheses $H^{0}_{j,k}$ and $G^{0}_{j,k}$.

We first discuss testing the unadjusted null hypothesis $H^{0}_{j,k}$ in (1), where the $\beta^{g,*}_{j,k}$ are the regression coefficients in (4). Suppose, for now, that we are in the low-dimensional setting, so the number of nodes $p$ is smaller than the sample sizes $n^{g}$, $g\in\{\mathrm{I},\mathrm{II}\}$.

It is well-known that the regression coefficients can be characterized as the minimizers of the expected least squares loss — that is,

$\displaystyle\boldsymbol{\beta}_{j}^{g,*}=(\beta^{g,*}_{j,1},\ldots,\beta^{g,*}_{j,p})^{\top}=\underset{\beta_{1},\ldots,\beta_{p}\in\mathds{R}}{\text{arg min}}\,\mathds{E}\left[\left(X^{g}_{j}-\sum_{k\neq j}X^{g}_{k}\beta_{k}\right)^{2}\right].$

One can obtain an estimate $\hat{\boldsymbol{\beta}}_{j}^{g}=(\hat{\beta}^{g}_{j,1},\ldots,\hat{\beta}^{g}_{j,p})$ of $\boldsymbol{\beta}^{g,*}_{j}=(\beta^{g,*}_{j,1}\ldots.,\beta^{g,*}_{j,p})$ by minimizing the empirical average of the least squares, taking

$\displaystyle\hat{\boldsymbol{\beta}}^{g}_{j}=\underset{\beta_{1},\ldots,\beta_{p}\in\mathds{R}}{\text{arg min}}\,\frac{1}{n^{g}}\left\|\mathbf{X}^{g}_{j}-\sum_{k\neq j}\mathbf{X}^{g}_{k}\beta_{k}\right\|_{2}^{2},$

where $\|\cdot\|_{2}$ denotes the $\ell_{2}$ norm. The ordinary least squares estimate $\hat{\boldsymbol{\beta}}^{g}_{j}$ is available in closed form and is easy to compute. The estimates $\hat{\beta}^{g}_{j,k}$ are unbiased, and, under mild assumptions, are approximately normally distributed for sufficiently large $n^{g}$ — that is,

$\displaystyle\hat{\beta}^{g}_{j,k}\sim N\left(\beta^{g,*}_{j,k},\tau^{g}_{j,k}\right),$

with $\tau^{g}_{j,k}>0$ (though $\tau^{g}_{j,k}$ can be calculated in closed form, we omit the expression for brevity).

We construct a test of $H_{j,k}^{0}$ based on the difference between the estimates of the group-specific regression coefficients, $\hat{\beta}^{\mathrm{I}}_{j,k}-\hat{\beta}^{\mathrm{II}}_{j,k}$. When $H^{0}_{j,k}$ holds, $\hat{\beta}^{\mathrm{I}}_{j,k}-\hat{\beta}^{\mathrm{II}}_{j,k}$ is normally distributed with mean zero and variance $\tau^{\mathrm{I}}_{j,k}+\tau^{\mathrm{II}}_{j,k}$. Given a consistent estimate $\hat{\tau}^{g}_{j,k}$ of the variance, we can use the test statistic

$\displaystyle T_{j,k}=\frac{\left(\hat{\beta}_{j,k}^{\mathrm{I}}-\hat{\beta}^{\mathrm{II}}_{j,k}\right)^{2}}{\hat{\tau}^{\mathrm{I}}_{j,k}+\hat{\tau}^{\mathrm{II}}_{j,k}},$

which follows a chi-square distribution with one degree of freedom under the null for $n^{\mathrm{I}}$ and $n^{\mathrm{II}}$ sufficiently large. A p-value for $H^{0}_{j,k}$ can be calculated as

$\displaystyle\rho_{j,k}=P\left(\chi^{2}_{1}>T_{j,k}\right).$

In the low-dimensional setting, performing a covariate-adjusted test is similar to performing the unadjusted test. We can obtain an estimate $\hat{\boldsymbol{\alpha}}^{g}_{j}=\left((\hat{\alpha}^{g}_{j,1})^{\top},\ldots,(\hat{\alpha}^{g}_{j,p})^{\top}\right)^{\top}$ of $\boldsymbol{\alpha}^{g,*}_{j}=\left((\alpha^{g,*}_{j,1})^{\top},\ldots,(\alpha^{g,*}_{j,p})^{\top}\right)^{\top}$ by minimizing the empirical average of the least squares loss

$\displaystyle\hat{\boldsymbol{\alpha}}^{g}_{j}=\underset{\alpha_{j,1},\ldots,\alpha_{j,p}\in\mathds{R}^{d}}{\text{arg min}}\,\frac{1}{n^{g}}\sum_{i=1}^{n^{g}}\left(X^{g}_{i,j}-\sum_{k\neq j}\left\langle\phi\left(W_{i}^{g}\right),\alpha_{j,k}\right\rangle X_{i,k}^{g}\right)^{2}.$

(10)

To simplify the presentation, we introduce additional notation that allows us to rewrite (10) in a condensed form. Let $\mathcal{V}^{g}_{k}$ be the $n^{g}\times d$ matrix

$\displaystyle\mathcal{V}^{g}_{k}=\begin{pmatrix}X_{1,k}^{g}\times\phi\left(W_{1}^{g}\right)\\ \vdots\\ X_{n^{g},k}^{g}\times\phi\left(W^{g}_{n^{g}}\right)\end{pmatrix}.$

(11)

We can now equivalently express (10) as

$\displaystyle\hat{\boldsymbol{\alpha}}^{g}_{j}=\underset{\alpha_{j,1},\ldots,\alpha_{j,p}\in\mathds{R}^{d}}{\text{arg min}}\,\frac{1}{n^{g}}\left\|\mathbf{X}^{g}_{j}-\sum_{k\neq j}\mathcal{V}^{g}_{k}\alpha_{j,k}\right\|_{2}^{2}.$

(12)

Again, $\hat{\alpha}^{g}_{j,k}$ is an unbiased and approximately normal for sufficiently large $n^{g}$, satisfying

$$\hat{\alpha}^{g}_{j,k}\sim N\left(\alpha^{g,*}_{j,k},\Omega^{g}_{j,k}\right),$$

where $\Omega^{g}_{j,k}$ is a positive definite matrix of dimension $d\times d$ (though a closed form expression is available, we omit it here for brevity).

We construct a test of $G^{0}_{j,k}$ based on $\hat{\alpha}^{\mathrm{I}}_{j,k}-\hat{\alpha}^{\mathrm{II}}_{j,k}$. Under the null hypothesis, $\hat{\alpha}^{\mathrm{I}}_{j,k}-\hat{\alpha}^{\mathrm{II}}_{j,k}$ follows a normal distribution with mean zero and variance $\Omega^{\mathrm{I}}_{j,k}+\Omega^{\mathrm{II}}_{j,k}$. Given a consistent estimate $\hat{\Omega}^{g}_{j,k}$ of $\Omega^{g}_{j,k}$, we can test $G^{0}_{j,k}$ using the test statistic

$\displaystyle S_{j,k}=\left(\hat{\alpha}_{j,k}^{\mathrm{I}}-\hat{\alpha}_{j,k}^{\mathrm{II}}\right)^{\top}\left(\hat{\Omega}^{\mathrm{I}}_{j,k}+\hat{\Omega}^{\mathrm{II}}_{j,k}\right)^{-1}\left(\hat{\alpha}_{j,k}^{\mathrm{I}}-\hat{\alpha}_{j,k}^{\mathrm{II}}\right).$

Under the null, the test statistic follows a chi-squared distribution with $d$ degrees of freedom, and a p-value can therefore be calculated as

$\displaystyle P\left(\chi^{2}_{d}>S_{j,k}\right).$

The methods described in Section 3.1 are only appropriate when the number of nodes $p$ is small relative to the sample size. Model (9) has $(p-1)d$ parameters, so the least squares estimator of Section 3.1 provides stable estimates as long as $n^{\mathrm{I}}$ and $n^{\mathrm{II}}$ are larger than $(p-1)d$. However, in the high-dimensional setting, where the the number of parameters exceeds the sample size, the ordinary least squares estimates behave poorly.

To fit the varying coefficient model (9) in the high-dimensional setting, we use a regularized estimator that relies upon an assumption of sparsity in the networks. The sparsity assumption requires that within each group only a small number of nodes are partially correlated, meaning that in (9), only a few of the vectors $\alpha^{g,*}_{j,k}$ are nonzero. To leverage the sparsity assumption, we propose to use the group LASSO estimator [46]:

$\displaystyle\tilde{\boldsymbol{\alpha}}_{j}^{g}=\underset{\alpha_{j,1},\ldots,\alpha_{j,p}\in\mathds{R}^{d}}{\text{arg min}}\,\frac{1}{n^{g}}\left\|\mathbf{X}^{g}_{j}-\sum_{k\neq j}\mathcal{V}^{g}_{k}\alpha_{j,k}\right\|_{2}^{2}+\lambda\sum_{k\neq j}\left\|\alpha_{j,k}\right\|_{2},$

(13)

where $\lambda>0$ is a tuning parameter. The group LASSO provides a sparse estimate and sets some $\tilde{\alpha}_{j,k}$ to be exactly zero, resulting in networks with few edges. The level of sparsity of $\tilde{\boldsymbol{\alpha}}^{g}_{j}$ is determined by $\lambda$, with higher $\lambda$ values forcing more $\tilde{\alpha}_{j,k}$ to zero. We discuss selection of the tuning parameter in Section 5.1.

Though the group LASSO provides a consistent estimate of $\boldsymbol{\alpha}^{g,*}_{j}$, the estimate is not approximately normally distributed. The group LASSO estimate of ${\alpha}^{g,*}_{j,k}$ retains a bias that diminishes at the same rate as the standard error. As a result, the group LASSO estimator has a non-standard sampling distribution that cannot be derived analytically and is therefore unsuitable for hypothesis testing.

We can obtain approximately normal estimates of $\alpha^{g,*}_{j,k}$ by correcting the bias of $\tilde{\alpha}^{g}_{j,k}$, as was first proposed to obtain normal estimates for the classical $\ell_{1}$-penalized version of the LASSO [12, 47]. These “de-biased” or “de-sparsified” estimators can been shown to be approximately normal with moderately large samples even in the high-dimensional setting; they are therefore suitable for hypothesis testing. Our approach is to use a de-biased version of the group LASSO. Bias correction in group LASSO problems is well studied [11, 16, 28], so we are able to perform covariate-adjusted inference by applying previously-developed methods.

The bias of the group LASSO estimate can be written as

$\displaystyle\delta^{g}_{j,k}=\mathds{E}\left[\tilde{\alpha}^{g}_{j,k}\right]-\alpha^{g,*}_{j,k},$

(14)

where $\delta^{g}_{j,k}$ is a nonzero $d$-dimensional vector (recall $d$ is the dimension of $\alpha^{g,*}_{j,k}$). Our approach is to obtain an estimate of the bias $\tilde{\delta}_{j,k}$ and to use a de-biased estimator, defined as

$\displaystyle\check{\alpha}^{g}_{j,k}=\tilde{\alpha}^{g}_{j,k}-\tilde{\delta}^{g}_{j,k}.$

(15)

For a suitable choice of $\tilde{\delta}_{j,k}$, the bias-corrected estimator is approximately normal for a sufficiently large sample size $n^{g}$ under mild conditions, i.e.,

$\displaystyle\check{\alpha}^{g}_{j,k}\sim N\left(\alpha^{g,*}_{j,k},\Omega^{g}_{j,k}\right),$

(16)

where the variance $\Omega^{g}_{j,k}$ is a positive definite matrix, for which we obtain an estimate $\check{\Omega}^{g}_{j,k}$. We provide a derivation for the bias-correction and the form of our variance estimate in Appendix A.

Similar to Section 2.1, we test the null hypothesis $G^{0}_{j,k}$ in (6) using the test statistic

$\displaystyle S_{j,k}=\left(\check{\alpha}_{j,k}^{\mathrm{I}}-\check{\alpha}_{j,k}^{\mathrm{II}}\right)^{\top}\left(\check{\Omega}_{j,k}^{\mathrm{I}}+\check{\Omega}_{j,k}^{\mathrm{II}}\right)^{-1}\left(\check{\alpha}_{j,k}^{\mathrm{I}}-\check{\alpha}_{j,k}^{\mathrm{II}}\right).$

(17)

The test statistic asymptotically follows a chi-squared distribution with $d$ degrees of freedom under the null hypothesis.

We begin by providing a brief summary of the score matching framework [17, 18]. Let $Z\in\mathcal{Z}\subseteq\mathds{R}^{p}$ be a random vector generated from a distribution with density function $h^{*}$. For any candidate density $h$, we denote the gradient and Laplician of the log-density by

$\displaystyle\nabla\log{h}(z)=\left\{\frac{\partial}{\partial z_{j}}\log h(x)\right\}\in\mathds{R}^{p};\quad\quad\Delta\log h(z)=\sum_{j=1}^{p}\frac{\partial^{2}}{\partial z_{j}^{2}}\log h(z_{j}).$

The score matching loss $L$ is defined as a measure of divergence between a candidate density function $h$ and the true density $h^{*}$:

$\displaystyle L(h)=\int\left\|\nabla\log h(z)-\nabla\log h^{*}(z)\right\|_{2}^{2}h^{*}(z)dz=\mathds{E}\left[\left\|\nabla\log h(Z)-\nabla\log h^{*}(Z)\right\|_{2}^{2}\right].$

(18)

It is apparent that the score matching loss is minimized when $h=h^{*}$. A natural approach to constructing an estimator for $h^{*}$ would then be to minimize the empirical score matching loss given observations $Z_{1},\ldots,Z_{n}$, defined as

$\displaystyle L_{n}(h)=\frac{1}{n}\sum_{i=1}^{n}\left\|\nabla\log h\left(Z_{i}\right)-\nabla\log h^{*}\left(Z_{i}\right)\right\|_{2}^{2}.$

Because the score matching loss function takes as input the gradient of the log density function, the loss does not depend on the normalizing constant. This makes score matching appealing when the normalizing constant is intractable.

The empirical loss seemingly depends on prior knowledge of $h^{*}$. However, if $h(z)$ and $\|h(z)\|_{2}$ both tend to zero as $z$ approaches the boundary of $\mathcal{Z}$, a partial integration argument can be used to show that the score matching loss can be expressed as

$\displaystyle L(h)=\int\left\{\Delta\log h(z)+\frac{1}{2}\left\|\nabla\log h(z)\right\|_{2}^{2}\right\}h^{*}(z)dz+\text{const.},$

(19)

where ‘const.’ is a term that does not depend on $h$. We can therefore estimate $h^{*}$ by minimizing an empirical version of the score matching loss that does not depend on $h^{*}$. We can express the empirical loss as

$\displaystyle L_{n}(h)=\frac{1}{n}\sum_{i=1}^{n}\Delta\log h(Z_{i})+\frac{1}{2}\left\|\nabla\log h(Z_{i})\right\|_{2}^{2}.$

The score matching loss is particularly appealing for exponential family distributions with continuous support, as it leads to a quadratic optimization function [23]. However, when $Z$ is non-negative, the arguments used to express (18) as (19) fail because $h(z)$ and $\|\nabla h(z)\|_{2}$ do not approach zero at the boundary. We can overcome this problem by instead considering the generalized score matching framework [44, 18] as an extension that is suitable for non-negative data. Let $v_{1},\ldots,v_{p}:\mathds{R}^{+}\to\mathds{R}^{+}$ be positive and differentiable functions, let $v(z)=\left(v_{1}(z_{1}),\ldots,v_{p}(z_{p})\right)^{\top}$, let $\dot{v}_{j}$ denote the derivative of $v_{j}$, and let $\circ$ denote the element-wise product operator. The generalized score matching loss is defined as

$\displaystyle L(h)=\int\left\|\left\{\nabla\log h(z)-\nabla\log h^{*}(z)\right\}\circ v^{1/2}(z)\right\|_{2}^{2}h^{*}(z)dz,$

(20)

and is also minimized when $h=h^{*}$. As for the original score matching loss (18), the generalized score matching loss seemingly depends on prior knowledge of $h^{*}$. However, under mild technical conditions on $h$ and $v$ (see Appendix B.1), the loss in (20) can be rewritten as

$\displaystyle L(h)=\int\bigg{[}\sum_{j=1}^{p}$

$\displaystyle\dot{v}_{j}(z_{j})\left\{\frac{\partial\log h(z_{j})}{\partial z_{j}}\right\}+v_{j}(z_{j})\left\{\frac{\partial^{2}\log h(z)}{\partial z^{2}_{j}}\right\}+\frac{1}{2}v_{j}(z_{j})\left\{\frac{\partial\log h(z)}{\partial z_{j}}\right\}^{2}\bigg{]}h^{*}(z)dz.$

(21)

The generalized score matching loss thus no longer depends on $h^{*}$, and an estimator can be constructed by minimizing the empirical version of (21) with respect to $h$. To this end, the original generalized score matching estimator considered $v_{j}(z_{j})=z_{j}^{2}$ [18]. In this case, it becomes necessary to estimate high moments of $h^{*}$, leading to poor performance of the estimator. It has been shown that by instead taking $v$ as a slowly increasing function, such as $v_{j}(z_{j})=\log(1+v_{j})$, one obtains improved theoretical results and better empirical performance [44].

In this sub-section, we discuss construction of asymptotically normal estimators for the parameters of the exponential family pairwise interaction model (7) using the generalized score matching framework. To simplify our presentation, we consider the setting in which we are only interested in studying the connectedness between one node $X^{g}_{j}$ and all other neighboring nodes in the network. To this end, it suffices to estimate the conditional density of $X^{g}_{j}$ given all other nodes and the covariates $W^{g}$. A similar approach to the one we describe below can also be used to estimate the entire joint density (7). For simplicity, we assume that in (7), there exist functions $\psi$ and $\zeta$ such that $\psi=\psi_{j,k}$ for all $(j,k)$ and $\zeta_{j,c}=\zeta$ for all $(j,c)$, and that $\mu_{j}=0$. For $x=(x_{1},\ldots,x_{p})^{\top}$ and $w=(w_{1},\ldots,w_{q})^{\top}$ the conditional density can thus be expressed as

$\displaystyle f_{j}^{g,*}(x_{j}|x_{1},\ldots,x_{p},w)\propto\exp\left(\sum_{j=1}^{p}\left\langle\alpha^{g,*}_{j,k},\phi\left(w\right)\right\rangle\psi(x_{j},x_{k})+\sum_{c=1}^{d}\theta_{j,c}^{g,*}\zeta\left(x_{j},\phi_{c}(w)\right)\right),$

(22)

where the density is up to a normalizing constant that does not depend on $x_{j}$.

We first explicitly define the score matching loss for the conditional density function (22). Let $\boldsymbol{\alpha}^{g,*}_{j}=\left((\alpha^{g,*}_{j,1})^{\top},\ldots,(\alpha^{g,*}_{j,p})^{\top}\right)^{\top}$, and similarly let $\boldsymbol{\theta}^{g,*}_{j}=(\theta^{g,*}_{j,1},\ldots,\theta^{g,*}_{j,p})^{\top}$. Let $\dot{\psi}$ and $\ddot{\psi}$ denote the first and second derivatives of $\psi$ with respect to $x_{j}$, and similarly, let $\dot{\zeta}$ and $\ddot{\zeta}$ denote the first and second derivatives of $\zeta$ with respect to $x_{j}$. We define a non-negative function $v_{j}:\mathds{R}_{+}\to\mathds{R}_{+}$, and let $\dot{v}_{j}$ denote the first derivative of $v_{j}$. Then for candidate parameters $\boldsymbol{\alpha}_{j}=\left(\alpha^{\top}_{j,1},\ldots,\alpha^{\top}_{j,p}\right)^{\top}$ and $\boldsymbol{\theta}_{j}=(\theta_{j,1},\ldots,\theta_{j,d})^{\top}$, the empirical generalized score matching loss for the conditional density of $X_{j}^{g}$ given all other nodes and the covariates can be expressed as

	$\displaystyle L^{g}_{n,j}(\boldsymbol{\alpha},\boldsymbol{\theta})=\frac{1}{2n^{g}}$	$\displaystyle\sum_{i=1}^{n^{g}}v_{j}\left(X_{i,j}^{g}\right)\bigg{\{}\sum_{k=1}^{p}\left\langle\alpha_{j,k},\phi\left(W^{g}_{i}\right)\right\rangle\dot{\psi}\left(X^{g}_{i,j},X^{g}_{i,k}\right)+\sum_{c=1}^{d}\theta_{j,c}\dot{\zeta}\left(X^{g}_{i,j},\phi_{c}\left(W^{g}_{i}\right)\right)\bigg{\}}^{2}+$
	$\displaystyle\frac{1}{n^{g}}$	$\displaystyle\sum_{i=1}^{n^{g}}v_{j}\left(X^{g}_{i,j}\right)\bigg{\{}\sum_{k=1}^{p}\left\langle\alpha_{j,k},\phi\left(W^{g}_{i}\right)\right\rangle\ddot{\psi}\left(X^{g}_{i,j},X^{g}_{i,k}\right)+\sum_{c=1}^{d}\theta_{j,c}\ddot{\zeta}\left(X^{g}_{i,j},\phi_{c}\left(W^{g}_{i}\right)\right)\bigg{\}}+$
	$\displaystyle\frac{1}{n^{g}}$	$\displaystyle\sum_{i=1}^{n^{g}}\dot{v}_{j}\left(X^{g}_{i,j}\right)\bigg{\{}\sum_{k=1}^{p}\left\langle\alpha_{j,k},\phi\left(W^{g}_{i}\right)\right\rangle\dot{\psi}\left(X^{g}_{i,j},X^{g}_{i,k}\right)+\sum_{c=1}^{d}\theta_{j,c}\dot{\zeta}\left(X^{g}_{i,j},\phi_{c}\left(W^{g}_{i}\right)\right)\bigg{\}}.$		(23)

The true parameters $\boldsymbol{\alpha}^{g,*}_{j}$ and $\boldsymbol{\theta}^{g,*}_{j}$ can characterized as the minimizers of the population score matching loss $\mathds{E}\left[L^{g}_{n,j}\left(\boldsymbol{\alpha}_{j},\boldsymbol{\theta}_{j}\right)\right]$, as discussed in Section 4.1.

The loss function in (23) is quadratic in parameters $\boldsymbol{\alpha}^{g,*}_{j}$ and $\boldsymbol{\theta}_{j}^{g,*}$ and can thus be solved efficiently. When the sample size $n^{g}$ is much larger than the number of unknown parameters $(p+1)d$, one can estimate $\boldsymbol{\alpha}^{g,*}_{j}$ and $\boldsymbol{\theta}_{j}^{g,*}$ by simply minimizing $L^{g}_{n,j}$ with respect to the unknown parameters. The empirical loss function is quadratic in $(\boldsymbol{\alpha}_{j},\boldsymbol{\theta}_{j})$, so the minimizer of the loss is available in closed form and can be computed efficiently. Moreover, we can readily establish asymptotic normality of the parameter estimates using results from classical M-estimation theory [34]. To avoid including cumbersome notation, we reserve the details for Appendix B.2.

When the sample size is smaller than the number of parameters, the minimizer of $L^{g}_{n,j}$ is no longer consistent. Similar to Section 3.2, we use regularization to obtain a consistent estimator in the high-dimensional setting. We define the $\ell_{2}$-regularized generalized score matching estimator as

$\displaystyle\left(\tilde{\boldsymbol{\alpha}}^{g}_{j},\tilde{\boldsymbol{\theta}}^{g}_{j}\right)=\underset{\boldsymbol{\alpha}_{j},\boldsymbol{\theta}_{j}}{\text{arg min}}\,\,L^{g}_{n,j}(\boldsymbol{\alpha}_{j},\boldsymbol{\theta}_{j})+\lambda\sum_{j=1}^{p}\left\|\alpha_{j,k}\right\|_{2},$

(24)

where $\lambda>0$ is a tuning parameter. Similar to the group LASSO estimator (13), the regularization term in (24) induces sparsity in the estimate $\tilde{\boldsymbol{\alpha}}_{j}^{g}$ and sets some $\tilde{\alpha}^{g}_{j,k}$ to be exactly zero. The tuning parameter controls the level of sparsity, where more vectors $\tilde{\alpha}^{g}_{j,k}$ are zero for higher $\lambda$. In Appendix B.3, we establish consistency of the regularized score matching estimator assuming sparsity of $\tilde{\boldsymbol{\alpha}}^{g}_{j}$ and some additional regularity conditions.

As is the case for the group LASSO estimator, the regularized score matching estimator has an intractable limiting distribution because its bias and standard error diminish at the same rate. We can obtain an asymptotically normal estimate by subtracting from the initial estimate an estimate of the bias. In Appendix B.4, we construct such a bias-corrected estimate $\check{\alpha}^{g}_{j,k}$ that, for sufficiently large $n^{g}$, satisfies

$\displaystyle\check{\alpha}^{g}_{j,k}\sim N\left(\alpha^{g,*}_{j,k},\Omega^{g}_{j,k}\right),$

for a positive definite matrix $\Omega^{g}_{j,k}$. Given bias-corrected estimates and a consistent estimate $\check{\Omega}^{g}_{j,k}$ of $\Omega^{g}_{j,k}$, we can test the null hypothesis (6) using the test statistic

$\displaystyle S_{j,k}=\left(\check{\alpha}^{\mathrm{I}}_{j,k}-\check{\alpha}^{\mathrm{II}}_{j,k}\right)^{\top}\left(\check{\Omega}^{\mathrm{I}}_{j,k}+\check{\Omega}^{\mathrm{II}}_{j,k}\right)^{-1}\left(\check{\alpha}^{\mathrm{I}}_{j,k}-\check{\alpha}^{\mathrm{II}}_{j,k}\right).$

Under the null hypothesis, the test statistic follows a chi-squared distribution with $d$ degrees of freedom.

We first discuss implementation of the neighborhood selection approach. The group LASSO estimate (13) does not exist in closed form, in contrast to the ordinary least squares estimate (12). To solve (13), we use the efficient algorithm implemented in the publicly available R package gglasso [42].

The group LASSO estimator requires selection of a tuning parameter $\lambda$, which controls the sparsity of the estimate. We select the tuning parameter by performing $K$-fold cross-validation, using $K=10$ folds. Since the selection of $\lambda$ is sensitive to the scale of the columns of $\mathcal{V}_{k}^{g}$ in (11), we scale the columns by their standard deviations prior to cross-validating. After fitting the group LASSO with the selected tuning parameter, we convert the estimates back to their original scale by dividing the estimates by the standard deviations of the columns of $\mathcal{V}_{k}^{g}$.

In what follows, we describe our simulation setting. In short, we generate data from the varying coefficient model (9), where we treat nodes $1$ through $(p-1)$ as predictors, and treat node $p$ as the response. We first randomly generate data for nodes $1$ through $(p-1)$ in groups $\mathrm{I}$ and $\mathrm{II}$ from the same multivariate normal distribution. We then construct $\eta^{g,*}_{j,k}$ and generate data for two covariates $W_{i}^{g}=(W_{i,1}^{g},W_{i,2}^{g})^{\top}$ so that one covariate acts as a confounding variable, and the other covariate should improve statistical power to detect differential associations after adjustment.

To simulate data for nodes $1$ through $(p-1)$, we first generate a random graph with $(p-1)$ nodes and an edge density of .05 from a power law distribution with power parameter 5 [30]. Denoting the edge set of the graph by $E$, we generate the $(p-1)\times(p-1)$ matrix $\Theta$ as

$$\Theta_{j,k}=\begin{cases}0&(j,k)\notin E\\ .5&(j,k)\in E\text{ with 50\% probability}\\ -.5&(j,k)\in E\text{ with 50\% probability}\end{cases},$$

with $\Theta_{j,k}=\Theta_{k,j}$. Defining by $a^{*}$ the smallest eigenvalue of $\Theta$, we set $\Sigma=(\Theta-(a^{*}-.1)I)^{-1}$, where $I$ is the identity matrix. We then draw $(X_{i,1}^{g},\ldots,X_{i,p-1}^{g})^{\top}$ from a multivariate normal distribution with mean zero and covariance $\Sigma$ for $i=1,\ldots,n^{g}$ for each group $g$.

We generate $W^{\mathrm{I}}_{i,1}$ from a $\text{Beta}(3/2,1)$ distribution and $W^{\mathrm{II}}_{i,1}$ from a $\text{Beta}(1,3/2)$ distribution. We center and scale both $W^{\mathrm{I}}_{i,1}$ and $W^{\mathrm{II}}_{i,1}$ to the $(-1,1)$ interval. We generate $W^{\mathrm{I}}_{i,2}$ and $W_{i,2}^{\mathrm{II}}$ each from a Uniform$(-1,1)$ distribution.

We consider two different choices for the varying coefficient functions $\eta^{g,*}_{j,k}$:

• •

  Linear Polynomial:

  	$\displaystyle\eta^{\mathrm{I,*}}_{p,1}(w_{1},w_{2})=.5+.5w_{1};$	$\displaystyle\eta^{\mathrm{II,*}}_{p,1}(w_{1},w_{2})=.5+.5w_{1}$
  	$\displaystyle\eta^{\mathrm{I,*}}_{p,2}(w_{1},w_{2})=.5+.25w_{2};$	$\displaystyle\eta^{\mathrm{II,*}}_{p,2}(w_{1},w_{2})=.5+.75w_{2}$
  	$\displaystyle\eta^{\mathrm{I,*}}_{p,3}(w_{1},w_{2})=0;$	$\displaystyle\eta^{\mathrm{II,*}}_{p,3}(w_{1},w_{2})=.5,$

  and $\eta^{g,*}_{p,k}=0$ for $k\geq 4$.

• •

  Cubic Polynomial:

  	$\displaystyle\eta^{\mathrm{I,*}}_{p,1}(w_{1},w_{2})=.5+.5\left(w_{1}+w_{1}^{2}+w_{1}^{3}\right);$	$\displaystyle\eta^{\mathrm{II,*}}_{p,1}(w_{1},w_{2})=.5+.5\left(w_{1}+w_{1}^{2}+w_{1}^{3}\right)$
  	$\displaystyle\eta^{\mathrm{I,*}}_{p,2}(w_{1},w_{2})=.5+.25\left(w_{2}+w_{2}^{3}\right);$	$\displaystyle\eta^{\mathrm{II,*}}_{p,2}(w_{1},w_{2})=.5+.75\left(w_{2}+w_{2}^{3}\right)$
  	$\displaystyle\eta^{\mathrm{I,*}}_{p,3}(w_{1},w_{2})=0;$	$\displaystyle\eta^{\mathrm{II,*}}_{p,3}(w_{1},w_{2})=.5,$

  and $\eta^{g,*}_{p,k}=0$ for $k\geq 4$.

The first covariate $W_{i,1}^{g}$ confounds the association between nodes $p$ and 1. The distribution of $W_{i,1}^{g}$ depends on group membership, and $W_{i,1}^{g}$ affects the association between genes $p$ and $1$. However, $\eta^{\mathrm{I},*}_{p,1}(w)=\eta^{\mathrm{II},*}_{p,1}(w)$ for all $w$. Thus, $G^{0}_{p,1}$ in (6) holds while $H^{0}_{p,1}$ in (1) fails, as depicted in Figure 1a. Failing to adjust for $W_{1}^{g}$ should therefore result in an inflated type-I error rate for the hypothesis $G^{0}_{p,1}$. Adjusting for the second covariate $W_{i,2}^{g}$ should improve the power to detect the differential connection between nodes $p$ and $2$. We have constructed $\eta^{g,*}_{p,2}$ so that $\mathds{E}\left[\eta^{\mathrm{I},*}\left(W^{\mathrm{I}}\right)\right]=\mathds{E}\left[\eta^{\mathrm{II},*}\left(W^{\mathrm{II}}\right)\right]$, though the association between nodes $p$ and 2 depends more strongly on $W^{g}$ in group $\mathrm{II}$ than in group $\mathrm{I}$. Thus, $H^{0}_{p,2}$ holds while $G^{0}_{p,2}$ fails, as depicted in Figure 1b. The association between nodes $p$ and $3$ does not depend on either covariate, though the association differs by group. Thus, one should be able to identify a differential connection using either the adjusted or unadjusted test. Node $p$ is conditionally independent of all other nodes in both groups.