Pseudolikelihood for Multilayer Graphical Models

Why Use a Product of Conditional Likelihoods?

A Gaussian graphical model represents conditional dependencies among measured features. In a multilayer omics network, an edge may connect two features in the same layer or connect features from different layers. Estimating the full joint model can be computationally demanding, particularly when the number of features is large relative to the number of samples. Pseudolikelihood offers an alternative: fit the conditional distribution of each feature given all the others, then combine those conditional likelihood contributions into one objective.

This construction connects naturally to nodewise regression, but its interpretation needs care. A product of conditional densities is generally not itself a normalized joint density. It is an estimation objective that uses local conditional information; it does not make the conditional contributions independent observations.

Definition: Let $X=(X_1,\ldots,X_p)$ be a random vector and let $x_i=(x_{i1},\ldots,x_{ip})$ denote sample $i$, for $i=1,\ldots,n$. The node-conditional pseudolikelihood is $$ PL(\theta)=\prod_{i=1}^{n}\prod_{j=1}^{p} f_{\theta,j}(x_{ij}\mid x_{i,-j}), $$ where $x_{i,-j}$ is the vector of all features except $j$, and $f_{\theta,j}$ is the model for feature $j$ conditional on the others. Its log is the sum of the node-conditional log-likelihood contributions. The parameter $\theta$ denotes the model parameters, which may include edge parameters shared across conditionals.

Gaussian Conditionals and Edge Parameters

For continuous features, suppose the rows are independent observations from a multivariate normal distribution with mean $\mu$ and positive-definite covariance matrix $\Sigma$. Write the precision matrix as $\Omega=\Sigma^{-1}$. For each feature $j$, the conditional distribution is

X_j\mid X_{-j}=x_{-j} \sim N\left(\mu_j-\frac{\Omega_{j,-j}}{\Omega_{jj}}(x_{-j}-\mu_{-j}),\ \frac{1}{\Omega_{jj}}\right).

Here the first argument of $N$ is the conditional mean and the second is the conditional variance. If the features are centered, the conditional mean can be written as $\sum_{k\ne j}\beta_{jk}x_k$, where $\beta_{jk}=-\Omega_{jk}/\Omega_{jj}$. For Gaussian graphical models, a conditional-independence edge between $j$ and $k$ is present exactly when $\Omega_{jk}\ne 0$. Thus a nonzero conditional regression coefficient indicates an edge in that node's conditional model. The coefficient itself is not the precision-matrix entry: the scale factor $\Omega_{jj}$ matters.

For a multilayer model, partition the feature indices into layer sets $V_1,\ldots,V_L$. Each possible pair $j,k$ belongs to a within-layer or cross-layer edge class according to which sets contain its endpoints. Pseudolikelihood supplies an estimation objective for these edge parameters; the layer labels do not change the conditional-density formula.

Model assumptions: The Gaussian formulas require a multivariate normal distribution with a positive-definite covariance matrix. The usual likelihood argument also assumes independent sample rows and an adequately specified conditional model. For a shared precision-matrix interpretation, the conditional distributions must be compatible with one common joint distribution. Penalization can make estimation feasible when $p$ is large, but does not remove the need to consider these assumptions.

How the Objective Is Built

For a centered Gaussian model, the contribution from node $j$ and sample $i$ is the log density of a normal regression with conditional mean $\sum_{k\ne j}\beta_{jk}x_{ik}$ and residual variance $\tau_j^2=1/\Omega_{jj}$. Ignoring constants that do not depend on the parameters, the negative log pseudolikelihood is

-\ell_{\mathrm{PL}}(\theta) = \sum_{j=1}^{p} \left[ \frac{n}{2}\log(\tau_j^2) + \frac{1}{2\tau_j^2} \sum_{i=1}^{n} \left(x_{ij}-\sum_{k\ne j}\beta_{jk}x_{ik}\right)^2 \right].

Without constraints, each node's regression contribution can be optimized as a regression problem. A network estimator often adds a penalty to encourage sparsity; for example, a penalty proportional to the absolute values of edge parameters can shrink weak associations to zero. A multilayer analysis can also use penalties that distinguish edge classes, but the resulting estimates still require a clear rule for enforcing or summarizing symmetry.

For a compatible Gaussian model, the two conditional coefficients for one pair satisfy

\beta_{jk}=-\frac{\Omega_{jk}}{\Omega_{jj}}, \qquad \beta_{kj}=-\frac{\Omega_{jk}}{\Omega_{kk}}.

They have the same zero pattern, but need not have the same numerical value when the conditional variances differ. If nodewise regressions are fitted separately, finite-sample estimates can disagree about whether a pair is selected. One must specify a rule, such as retaining an edge if either conditional estimate is nonzero (an OR rule) or only if both are nonzero (an AND rule), or instead optimize a parameterization that imposes shared symmetric edge parameters. These choices can yield different networks.

Worked Example: Conditional Likelihood Is Not the Joint Likelihood

Consider a hypothetical centered bivariate Gaussian model with covariance matrix $\Sigma=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}$. Its determinant is $1-0.25=0.75$, so

\Omega=\Sigma^{-1} =\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix} =\begin{pmatrix}4/3&-2/3\\-2/3&4/3\end{pmatrix}.

Each conditional variance is $1/(4/3)=0.75$, and each conditional mean is $-(-2/3)/(4/3)$ times the other feature, or $0.5$ times that feature. Evaluate the two conditional densities at the hypothetical observation $x=(1,0)$. The residuals are $1-0.5(0)=1$ for the first node and $0-0.5(1)=-0.5$ for the second. Their product has log

\log PL(x) =-\log(2\pi\cdot0.75) -\frac{1^2+(-0.5)^2}{2(0.75)} =-\log(4.7124)-0.8333 \approx -2.384.

The joint log density at the same observation is

\log f(x) =-\log(2\pi)-\frac{1}{2}\log(0.75) -\frac{1}{2}\left(\frac{4}{3}\right) \approx -2.361.

The two values differ, even though the conditionals come from a valid common Gaussian model. The product of conditionals is therefore not the joint density in this example. Its role is to provide a tractable objective for estimating model parameters.

Worked Estimation and Interpretation

Worked Example: Estimating a Cross-Layer Edge in Two Features

Suppose three hypothetical samples have centered values for one RNA feature $X_1$ and one protein feature $X_2$: $(1,1)$, $(-1,0)$, and $(0,-1)$. Assume independent observations from a centered, nonsingular bivariate Gaussian model. The sample means are both zero. Using the Gaussian maximum-likelihood covariance estimate, which divides the centered cross-product matrix by $n=3$, gives

\widehat{\Sigma} =\frac{1}{3} \begin{pmatrix} 1^2+(-1)^2+0^2 & 1(1)+(-1)(0)+0(-1)\\ 1(1)+(-1)(0)+0(-1) & 1^2+0^2+(-1)^2 \end{pmatrix} = \begin{pmatrix}2/3&1/3\\1/3&2/3\end{pmatrix}.

Its determinant is $(2/3)(2/3)-(1/3)^2=1/3$. Inverting gives

\widehat{\Omega} = \begin{pmatrix}2&-1\\-1&2\end{pmatrix}.

The nodewise regression slope for $X_1$ on $X_2$ is the sample covariance divided by the variance of $X_2$: $(1/3)/(2/3)=0.5$. The reverse slope is also $0.5$. The conditional variance for either node is $2/3-(1/3)^2/(2/3)=2/3-1/6=0.5$. These calculations agree with the precision-matrix relations: $-\widehat{\Omega}_{12}/\widehat{\Omega}_{11}=1/2$, and $1/\widehat{\Omega}_{11}=1/2$.

Because one feature is RNA and the other is protein, the nonzero off-diagonal precision entry represents a fitted cross-layer conditional-dependence edge. This tiny example illustrates the parameter mapping, not reliable edge inference: three samples provide little information about uncertainty, and no claim of biological causation follows from the fitted association.

Penalization and a Conditional Regression Calculation

When there are many features, unpenalized regressions can produce dense estimates. An absolute-value penalty can shrink conditional coefficients and set some to zero. The following calculation isolates one conditional regression contribution to show how the penalty acts; it is not, by itself, a complete compatible network fit.

Worked Example: A Lasso-Penalized Conditional Slope

Take four hypothetical centered predictor values $x=(1,-1,1,-1)$, and let $z=(1,1,-1,-1)$. Define the response values by $y=0.4x+\sqrt{0.84}\,z$, so $x^\mathsf{T}x/4=1$ and $x^\mathsf{T}y/4=0.4$, since $x^\mathsf{T}z=0$. For a conditional regression with residual variance fixed at 1, minimize the negative log-likelihood per observation plus an absolute-value penalty:

Q(\beta)=\frac{1}{8}\sum_{i=1}^{4}(y_i-\beta x_i)^2+0.15|\beta|.

The derivative of the squared-error term is $\beta-x^\mathsf{T}y/4=\beta-0.4$. For a positive coefficient, setting the derivative of $Q$ to zero gives $\beta-0.4+0.15=0$, hence $\widehat{\beta}=0.25$. For a zero coefficient, the absolute-value penalty can offset a score only when $|0.4|\leq0.15$, which is false here; the positive solution is therefore the minimizer. The penalty reduces the unpenalized slope $0.4$ to $0.25$.

This is a calculation for one node-conditional regression with a specified residual variance. In a full Gaussian network, conditional variances and the consistency of parameters across nodes also matter. A fitted nonzero conditional slope should not automatically be reported as a symmetric precision-matrix entry.

Justification, Diagnostics, and Alternatives

The computational appeal follows from the form of the objective: each term uses one response feature and the remaining features as predictors. Rather than repeatedly evaluating a full $p$-dimensional joint density, an algorithm can work with a collection of conditional regression problems. This can be useful when $p$ is large or the joint likelihood optimization is difficult. It does not make the data contributions independent across nodes: for the same sample, every conditional term uses the same feature vector.

Under a correctly specified, identifiable model with independent observations and suitable regularity conditions, maximizing a compatible pseudolikelihood can yield consistent estimates as the sample size grows. In high dimensions, consistency of a penalized estimator additionally depends on conditions such as adequate signal strength, sparsity, and suitable design behavior. The penalty trades bias for a sparser, more stable estimate; it does not guarantee that every selected edge is real.

Useful checks include examining conditional residuals for gross departures from the assumed Gaussian form, checking whether estimates are sensitive to feature scaling and penalty choice, and comparing the two nodewise estimates for each undirected pair. Large directional disagreements can signal sampling variability, incompatible separately fitted conditionals, or a scaling problem. A positive-definite covariance model is also required for the stated Gaussian joint interpretation. These checks can reveal problems but do not prove that the model is correct.

Approach	What it optimizes	Main consideration
Full Gaussian likelihood	A joint density parameterized by a covariance or precision matrix	Maintains the joint model directly, but may be computationally demanding.
Nodewise regression	Separate conditional regressions, often with sparsity penalties	Computationally convenient; estimates for the two ends of an edge may disagree.
Pseudolikelihood	The product, or sum of log contributions, of node conditionals	Uses all nodewise conditional information, but is not generally a joint density and needs careful edge and uncertainty interpretation.

Pseudolikelihood is not a universal replacement for joint likelihood. It can be less statistically efficient than a correctly specified full likelihood, and naive uncertainty calculations may be wrong because the conditional contributions from a sample are dependent. In particular, treating the product as an ordinary joint density can lead to unjustified likelihood-ratio tests or standard errors. The interpretation should remain conditional association, not causation.

Common Mistakes

Calling the product a joint density. The node conditionals overlap in the same observations and their product need not integrate to one. Treat pseudolikelihood as an estimation objective unless a separate result establishes a joint-density identity.
Equating a regression slope with a precision entry. In the Gaussian model, $\beta_{jk}=-\Omega_{jk}/\Omega_{jj}$. The slope and precision entry have different scales and opposite signs.
Assuming two separately fitted slopes must match. Even in a compatible population model their magnitudes can differ because the conditional variances differ. Finite-sample and penalized estimates can also have different zero patterns.
Interpreting an edge causally. A Gaussian graphical edge encodes conditional dependence under the model; it does not establish a causal effect.

Key takeaway: Gaussian pseudolikelihood combines node-conditional likelihood contributions to make estimation more tractable. Interpret edges through the conditional-independence structure, account for the relation between regression coefficients and precision parameters, and do not mistake the objective for a joint density.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a Gaussian graphical model, what does $\Omega_{jk}=0$ imply about $X_j$ and $X_k$ conditional on the other features?
Why is the product of node-conditional densities generally not the joint density?
If $\Omega_{jk}=-0.6$ and $\Omega_{jj}=2$, what is the conditional regression coefficient $\beta_{jk}$?
What problem can arise when two nodewise regressions give different edge-selection results?
Does a selected pseudolikelihood edge establish that one omics feature causes the other? Why or why not?

Solutions. (1) The two features are conditionally independent given all other features. (2) The conditional factors reuse the same feature vector and their product is not generally normalized as a joint density. (3) $\beta_{jk}=-(-0.6)/2=0.3$. (4) The undirected edge may be present in one conditional estimate but absent in the other, so the analyst needs a stated symmetry or aggregation rule. (5) No; it represents conditional association under the model and does not establish causation.