Worked Example Fitting Two Matched Omics Blocks

Fitting a Common Factor Across Matched Blocks

Suppose each subject has measurements from two omics blocks, such as expression and protein abundance. A factor model can summarize a pattern that varies across subjects and appears in both blocks. The practical task is to fit the pattern, quantify each subject’s position along it, and check how well it reconstructs the observed measurements. Matching is essential: each row in both blocks must refer to the same subject in the same order.

For this worked model, measurements are centered across subjects and placed on comparable scales before fitting. Write \(X^{(1)}\) and \(X^{(2)}\) for the resulting \(n\)-by-\(p_1\) and \(n\)-by-\(p_2\) matrices. Concatenate them by columns to form \(X=[X^{(1)}\;X^{(2)}]\), an \(n\)-by-\(p\) matrix where \(p=p_1+p_2\). A row is a subject; a column is a feature.

The Gaussian assumption gives a likelihood-based reason to minimize squared residuals: for independent Gaussian errors with common variance, maximizing the conditional likelihood over scores and loadings is equivalent to minimizing \(\sum_{i,j}(X_{ij}-z_i\lambda_j)^2\). For a single-factor fit, use the normalization \(\sum_i z_i=0\) and \(z'z=n\). Without a scale constraint, multiplying all scores by a constant and dividing all loadings by the same constant leaves the fitted values unchanged. A sign change in both scores and loadings also leaves them unchanged.

Least-Squares Fit and Its Justification

For any fixed, centered score vector \(z\) with \(z'z=n\), the loading that minimizes the squared residuals is \(\hat{\lambda}=X'z/n\). The fitted matrix is then \(\hat{X}=z\hat{\lambda}'\). Searching over all allowed score vectors gives the best rank-one least-squares reconstruction of \(X\). This is the leading-component solution from the singular value decomposition (SVD), with the score vector rescaled to satisfy the stated normalization.

To see why the SVD gives the least-squares solution, take a unit-length vector \(u\) and write \(z=\sqrt{n}\,u\). For fixed \(u\), the fitted value in feature \(j\) is the projection of column \(X_{\cdot j}\) onto \(u\), namely \(u(u'X_{\cdot j})\). The total squared residual is therefore

The first term does not depend on \(u\). Thus minimizing residual sum of squares is equivalent to maximizing \(u'XX'u\) subject to \(u'u=1\). The maximum is attained by a leading eigenvector of \(XX'\), equivalently the leading left singular vector of \(X\). This establishes the rank-one least-squares result; it does not establish that one factor is scientifically adequate or that the Gaussian assumptions are appropriate.

Scores, Reconstruction, and Residuals

After fitting the loadings, reconstruct a subject’s measurements by multiplying their score by each loading. A positive score predicts values in the directions of positive loadings and opposite directions for negative loadings; a negative score reverses those directions. A loading’s magnitude indicates how strongly that feature changes along the fitted pattern, in the units used for the input data.

There are different ways to estimate scores once loadings are available. With scores treated as unknown fixed quantities, the least-squares score for a subject is \(\hat{z}_i=X_{i\cdot}\hat{\lambda}/(\hat{\lambda}'\hat{\lambda})\). If instead the model assigns a standard normal prior \(z_i\sim N(0,1)\), the posterior mean score is \(\hat{z}_i=\hat{\lambda}'X_{i\cdot}/(\hat{\lambda}'\hat{\lambda}+\sigma^2)\). The latter shrinks scores toward zero to reflect residual noise. These estimators answer slightly different questions; state which is used.

Assumptions, Diagnostics, and Alternatives

A useful interpretation depends on more than a successful numerical fit. Check that the rows in both blocks refer to the same subjects, that the features were centered consistently, and that scaling has not allowed a high-variance block to dominate the joint fit. The spherical Gaussian model assumes independent residuals with a common variance across features and subjects. Correlated residuals, strong outliers, or markedly different feature-specific noise levels can make its scores and loadings misleading.

Inspect residuals \(X-\hat{X}\) by subject and feature, and compare residual variability between blocks. Look for structure remaining in the residuals, such as a feature group that consistently moves together or a subject with unusually large residuals. Check whether the score and loading patterns are stable to reasonable preprocessing choices. These checks help detect misspecification; a small SSE alone does not show that the factor is stable or that a one-factor model is adequate. The earlier material on residual diagnostics and overfitting in latent factor models provides broader context for these checks.

The rank-one least-squares fit is closely related to principal component analysis: both select a direction that explains as much total squared variation as possible. A probabilistic factor model adds an explicit error distribution and can support uncertainty calculations, but the simple spherical version still imposes restrictive residual assumptions. Alternatives may allow feature-specific residual variances, correlated residuals, or robust error distributions. Other integration methods can target cross-block association rather than total reconstruction error. They need not return the same scores: the result depends on the model and objective, not just on the data.

A common mistake is to treat a positive score as evidence that a biological process is “activated,” or to read a large loading as proof of a causal effect. The correction is to describe the fitted pattern in terms of the observed, processed features: for example, “subjects with higher scores tend to have higher centered values for these positively loaded features and lower values for these negatively loaded features.” A shared pattern across two blocks is an association in these measurements. It does not by itself establish a mechanism, direction of influence, or clinical meaning.