Deriving EM for a Gaussian Multi-Omics Model

Why Use EM for Whole Missing Modalities?

In a multi-omics study, a subject may have measurements for one modality but none for another. Estimating a separate covariance from each observed feature pair can produce estimates that do not form a valid covariance matrix. The earlier material on the observed-data likelihood showed how each subject contributes a marginal Gaussian likelihood for the coordinates that are available. The expectation-maximization (EM) algorithm offers another route: it repeatedly calculates the expected complete-data moments under the current model, then updates the model using those moments.

Here the target is a single multivariate Gaussian model for the joint feature vector. The method does not require filling in one fixed set of guessed values for missing measurements. Instead, the E-step accounts for uncertainty in missing coordinates through their conditional variances and covariances.

Definition: Let \(Y_i\) be the \(d\)-dimensional vector of all features for subject \(i\), with \(Y_i\sim N_d(\mu,\Sigma)\), independently across \(i=1,\ldots,n\). Here \(\mu\) is the mean vector and \(\Sigma\) is a positive-definite covariance matrix. Write \(O_i\) for the indices observed for subject \(i\), and \(M_i\) for the missing indices. The whole-modality pattern means that these sets contain or omit modality blocks, rather than arbitrary individual feature values. The parameters to estimate are \(\mu\) and \(\Sigma\).

The observed-data likelihood depends on \(\mu\) and on the appropriate marginal submatrices of \(\Sigma\). EM treats the unavailable entries of \(Y_i\) as unobserved complete data, not as known values. The resulting updates rely on the conditional distribution under the current parameter values.

Model and Assumptions

The Gaussian model assumes that subjects are independent draws from one multivariate normal population with a common mean and covariance. Thus, the model allows features from different modalities to be correlated, but it does not by itself account for batches, covariates, or distinct subpopulations. Those would require an expanded model.

For likelihood-based inference from the observed values to be valid without explicitly modelling the availability mechanism, missingness must be ignorable: the mechanism is missing completely at random, or missing at random given observed information included in the analysis, and its parameters are distinct from \(\mu\) and \(\Sigma\). Structural availability can be ignorable in this sense when its design is independent of the unobserved measurements. If availability depends on unobserved values even after conditioning on observed information, the Gaussian EM updates alone do not address that dependence.

A further requirement is that the observed patterns contain enough information to identify the parameters being estimated. In particular, estimating a cross-modality covariance requires information linking the modalities, such as subjects with both modalities observed, or other valid model restrictions. The EM formulas can still be written when a parameter is not identifiable, but the data will not determine a unique estimate of it.

The E-Step: Expected Complete-Data Moments

For subject \(i\), let \(m_i=E(Y_i\mid Y_{i,O_i};\mu^{(t)},\Sigma^{(t)})\) and \(V_i=\operatorname{Var}(Y_i\mid Y_{i,O_i};\mu^{(t)},\Sigma^{(t)})\), where the superscript \(t\) denotes the current iteration. The conditional expectation of the outer product is \(E(Y_iY_i'\mid Y_{i,O_i})=m_im_i'+V_i\). This identity follows by writing \(Y_i=m_i+(Y_i-m_i)\): the conditional mean of the residual is zero, so the cross-product terms have conditional expectation zero, while the residual outer product has expectation \(V_i\).

For a subject with observed coordinates \(O_i\) and missing coordinates \(M_i\), partition the current mean and covariance accordingly. The observed entries of \(m_i\) equal their recorded values, and the observed rows and columns of \(V_i\) are zero. The missing entries and their conditional covariance are

m_{i,M_i}=\mu^{(t)}_{M_i}+ \Sigma^{(t)}_{M_iO_i} \left(\Sigma^{(t)}_{O_iO_i}\right)^{-1} \left(Y_{i,O_i}-\mu^{(t)}_{O_i}\right), \qquad V_{i,M_iM_i}=\Sigma^{(t)}_{M_iM_i}- \Sigma^{(t)}_{M_iO_i} \left(\Sigma^{(t)}_{O_iO_i}\right)^{-1} \Sigma^{(t)}_{O_iM_i}.

These are the conditional mean and covariance of a partitioned multivariate normal distribution. They require the observed covariance block to be invertible. With a positive-definite \(\Sigma^{(t)}\), every nonempty principal submatrix is positive definite and therefore invertible. If no coordinates are missing, \(m_i=Y_i\) and \(V_i=0\); if none are observed, the conditional moments are simply the current \(\mu^{(t)}\) and \(\Sigma^{(t)}\), though such a subject contributes no observed-data likelihood information.

Worked Example: Conditional Moments for One Missing Feature

Suppose the current mean and covariance for two features, A and B, are \(\mu^{(t)}=(1,0)'\) and \(\Sigma^{(t)}=\begin{pmatrix}4&2\\2&3\end{pmatrix}\). For a subject with A observed as 3 and B missing, calculate the conditional moments needed by the E-step.

The conditional mean of B is \(0+(2/4)(3-1)=1\). Its conditional variance is \(3-(2/4)2=2\). Thus \(m_i=(3,1)'\), and \(V_i\) has zero entries for the observed A coordinate and variance 2 for B: \(V_i=\begin{pmatrix}0&0\\0&2\end{pmatrix}\).

The expected outer product is \(m_im_i'+V_i=\begin{pmatrix}9&3\\3&1\end{pmatrix}+\begin{pmatrix}0&0\\0&2\end{pmatrix}=\begin{pmatrix}9&3\\3&3\end{pmatrix}\). The value 3 in the B,B entry includes both the squared conditional mean, \(1^2\), and the conditional variance, 2. Omitting that variance would understate the expected second moment.

The M-Step and Its Derivation

For complete data, the Gaussian log-likelihood, apart from a constant, is \(-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\sum_i(Y_i-\mu)'\Sigma^{-1}(Y_i-\mu)\). In EM, the M-step maximizes its conditional expectation given the observed data and the current parameters. The expected complete-data log-likelihood is therefore

Q(\mu,\Sigma\mid\mu^{(t)},\Sigma^{(t)}) =C-\frac{n}{2}\log|\Sigma| -\frac{1}{2}\sum_{i=1}^{n} E\!\left[(Y_i-\mu)'\Sigma^{-1}(Y_i-\mu)\mid Y_{i,O_i}\right],

where \(C\) does not depend on the proposed \(\mu\) or \(\Sigma\). Define \(S_1=\sum_i m_i\) and \(S_2=\sum_i(m_im_i'+V_i)\), using E-step moments computed under the current parameters. Expanding the quadratic term shows that the part depending on \(\mu\) is minimized at \(\mu=S_1/n\). Substituting this mean, the covariance-dependent part is maximized by matching \(\Sigma\) to the average expected centered outer product. The updates are

\mu^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}m_i, \qquad \Sigma^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n} \left\{V_i+(m_i-\mu^{(t+1)})(m_i-\mu^{(t+1)})'\right\}.

Equivalently, the covariance update is \(S_2/n-\mu^{(t+1)}\mu^{(t+1)'}\). The denominator is \(n\), as in the maximum-likelihood estimate for a Gaussian model, not \(n-1\), which is used for the usual unbiased sample covariance when complete observations are available.

Worked Example: One EM Iteration With Two Modalities

Consider three hypothetical subjects, each with one measurement for modality A and one for modality B. The observations are: subject 1 has A=1 and B=2; subject 2 has A=0 and B missing; subject 3 has A missing and B=1. Start EM at \(\mu^{(t)}=(0,0)'\) and \(\Sigma^{(t)}=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\). This covariance is positive definite because its determinant is \(1-0.5^2=0.75>0\).

For subject 1, both values are observed, so \(m_1=(1,2)'\), \(V_1=0\), and \(m_1m_1'+V_1=\begin{pmatrix}1&2\\2&4\end{pmatrix}\).

For subject 2, A=0 is observed. The conditional mean of B is \(0+0.5(0-0)=0\), and its conditional variance is \(1-0.5^2=0.75\). Thus \(m_2=(0,0)'\), and \(m_2m_2'+V_2=\begin{pmatrix}0&0\\0&0.75\end{pmatrix}\).

For subject 3, B=1 is observed. The conditional mean of A is \(0+0.5(1-0)=0.5\), and its conditional variance is again 0.75. Thus \(m_3=(0.5,1)'\), and \(m_3m_3'+V_3=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\).

The updated mean is \(\mu^{(t+1)}=((1+0+0.5)/3,(2+0+1)/3)'=(0.5,1)'\). The sum of expected outer products is \(\begin{pmatrix}2&2.5\\2.5&5.75\end{pmatrix}\), so \(S_2/3=\begin{pmatrix}2/3&5/6\\5/6&23/12\end{pmatrix}\). Subtracting \(\mu^{(t+1)}\mu^{(t+1)'}=\begin{pmatrix}1/4&1/2\\1/2&1\end{pmatrix}\) gives \(\Sigma^{(t+1)}=\begin{pmatrix}5/12&1/3\\1/3&11/12\end{pmatrix}\), approximately \(\begin{pmatrix}0.4167&0.3333\\0.3333&0.9167\end{pmatrix}\).

The determinant from the unrounded entries is \((5/12)(11/12)-(1/3)^2=55/144-16/144=39/144=13/48\approx0.2708>0\). The diagonal entries are positive, so this two-dimensional covariance update is positive definite. The iteration has updated both modalities jointly while retaining the conditional variance contributed by the missing measurements.

Why the Likelihood Does Not Decrease

The EM update has a useful monotonicity property under the stated model and regularity conditions. Let \(q^{(t)}\) be the conditional distribution of the complete data given the observed data under the current parameter values. For any candidate parameters \(\theta=(\mu,\Sigma)\), the observed-data log-likelihood satisfies the identity

\log p_{\theta}(Y_{\mathrm{obs}}) = E_{q^{(t)}}[\log p_{\theta}(Y_{\mathrm{complete}})] +H(q^{(t)}) +\operatorname{KL}\!\left(q^{(t)}\,\middle\|\,p_{\theta}(Y_{\mathrm{mis}}\mid Y_{\mathrm{obs}})\right),

where \(H(q^{(t)})\) is the entropy of the current conditional distribution and \(\operatorname{KL}\) is the Kullback–Leibler divergence, which is nonnegative. At \(\theta=\theta^{(t)}\), the current conditional distribution equals the model's conditional distribution, so the divergence is zero. The M-step increases or leaves unchanged the expected complete-data log-likelihood, while \(H(q^{(t)})\) is fixed during that step. Since the divergence at the new parameters cannot be negative, the observed-data log-likelihood at the new parameters is at least as large as at the old parameters. This guarantees nondecrease, not a unique or global maximum in every problem.

Worked Example: An Availability Pattern That Cannot Identify Covariance

Suppose a hypothetical study has subjects with A observed only and other subjects with B observed only; no subject has both modalities measured. To isolate the issue, take a bivariate Gaussian with means zero and variances one, and covariance \(\rho\), where \(-1<\rho<1\). An A-only observation \(a\) has marginal density \(N(0,1)\), and a B-only observation \(b\) also has marginal density \(N(0,1)\). Neither marginal density contains \(\rho\).

For example, if the observed A values are 0 and 1 and the observed B values are -1 and 1, the observed-data likelihood is the product of four standard normal densities at those values. It is exactly the same for every \(\rho\) in \((-1,1)\). Thus, the available data cannot distinguish covariance 0 from covariance 0.5, or from any other admissible value. An EM run may return a covariance depending on its starting value or numerical details, but that does not make the covariance data-identified. The missing issue is information linking A and B, not a failure of the arithmetic in the update.

Diagnostics, Alternatives, and Limitations

Before interpreting an EM estimate, inspect the availability patterns and the numerical behavior of the iterations. Useful checks include:

Pattern support: determine which modality pairs are jointly observed and whether these overlaps connect the blocks whose covariance is of interest. A disconnected or absent overlap can leave cross-block dependence unidentified.
Numerical validity: check that \(\Sigma^{(t)}\) remains symmetric and positive definite, and monitor the observed-data log-likelihood. The likelihood should not decrease beyond small numerical tolerance under a correctly implemented exact EM step.
Convergence and initialization: inspect changes in the likelihood and parameter estimates, and try multiple valid starting values when the likelihood may have multiple stationary points or a weakly identified direction.
Model fit and missingness: assess whether a single Gaussian distribution is plausible for the observed data, and consider whether the availability mechanism could depend on unobserved values. EM does not repair a misspecified distribution or a nonignorable missingness mechanism.

Complete-case analysis uses only subjects with every modality observed. It is simple, but can discard substantial information and can be biased if the retained subjects are not representative. Pairwise covariance estimation uses different subsets for different feature pairs; it can be convenient descriptively, but the resulting entries need not form a positive-semidefinite matrix, so they may not define any joint Gaussian model. Gaussian EM instead uses all available observations through a coherent joint model, but its validity depends on the model and identifiability assumptions.

The covariance update is positive semidefinite because it is an average of expected centered outer products. It is positive definite only when the accumulated information spans all feature directions; sparse or poorly connected patterns can yield a singular update. In high dimensions, the number of covariance parameters can also exceed what the observed patterns support. Regularization or a lower-dimensional model may be needed, but those choices change the estimation problem and are not supplied automatically by EM. Finally, a converged estimate is not evidence that the missingness assumptions or Gaussian model are correct.

Common pitfall: Replacing each missing value by its conditional mean and then calculating an ordinary covariance treats the imputed value as certain. The correction is to include \(V_i\) in each expected second moment: \(E(Y_iY_i'\mid Y_{i,O_i})=m_im_i'+V_i\). This term preserves uncertainty about the missing measurements.

Check Your Understanding

Answer each question, then compare with the solutions below.

For one subject, why is \(E(Y_iY_i'\mid Y_{i,O_i})\) generally not just \(m_im_i'\)?
In the two-modality iteration, what would the conditional variance of a missing feature be if the current correlation were zero and both variances were one?
Why can no subjects jointly observed for A and B prevent estimation of their covariance when only A-only and B-only patterns occur?
State one reason the EM observed-data likelihood might fail to be nondecreasing in a numerical implementation.

Solutions. (1) The conditional second moment includes both the outer product of the conditional mean and the conditional covariance: \(E(Y_iY_i'\mid Y_{i,O_i})=m_im_i'+V_i\). The variance term is positive when missing values remain uncertain. (2) With zero correlation, the conditional mean does not use the observed feature and the conditional variance remains its marginal variance, 1. (3) The observed likelihood is then a product of separate marginal likelihoods; neither marginal depends on the cross-modality covariance, so different covariance values yield the same observed-data likelihood. (4) Possible causes include an incorrect conditional-moment calculation, an invalid or poorly conditioned covariance update, or numerical error in evaluating matrix inverses or determinants. Exact EM under its assumptions should not decrease the observed-data likelihood.