Integrative Factorization of Bidimensionally Linked Matrices

TL;DR

BIDIFAC is a new method for integrating and analyzing bidimensionally linked biomedical data matrices, capturing shared and unique variability across multiple cohorts and platforms.

Contribution

It introduces a novel bidimensional factorization approach that extends nuclear norm penalization for complex multi-cohort, multi-platform data integration.

Findings

Successfully applied to breast cancer genomics data from TCGA

Effectively separates shared and unique data patterns

Provides R code for implementation and simulation

Abstract

Advances in molecular "omics'" technologies have motivated new methodology for the integration of multiple sources of high-content biomedical data. However, most statistical methods for integrating multiple data matrices only consider data shared vertically (one cohort on multiple platforms) or horizontally (different cohorts on a single platform). This is limiting for data that take the form of bidimensionally linked matrices (e.g., multiple cohorts measured on multiple platforms), which are increasingly common in large-scale biomedical studies. In this paper, we propose BIDIFAC (Bidimensional Integrative Factorization) for integrative dimension reduction and signal approximation of bidimensionally linked data matrices. Our method factorizes the data into (i) globally shared, (ii) row-shared, (iii) column-shared, and (iv) single-matrix structural components, facilitating the investigation of shared and unique patterns of variability. For estimation we use a penalized objective function that extends the nuclear norm penalization for a single matrix. As an alternative to the complicated rank selection problem, we use results from random matrix theory to choose tuning parameters. We apply our method to integrate two genomics platforms (mRNA and miRNA expression) across two sample cohorts (tumor samples and normal tissue samples) using the breast cancer data from TCGA. We provide R code for fitting BIDIFAC, imputing missing values, and generating simulated data.

Tables

Table 1. Table 1: Summary of the simulation studies for Design 1, where prediction and imputation errors are computed using Equations ( 12 ) and ( 13 ) respectively.

		Prediction							Imputation
SNR	Model	$G$	$R$	$C$	$I$	$G+C$	$R+I$	$S$	Cell	Column	Row
0.5	BIDIFAC	$-$	$-$	0.81	0.87	0.80	0.87	0.81	0.86	0.95	1.00
	UNIFAC	$-$	$-$	0.80	0.87	0.80	0.87	0.81	0.86	0.94	1.00
	JIVE(T)	$-$	$-$	1.01	1.40	1.01	1.40	1.05	1.49	1.13	1.02
	JIVE(P)	$-$	$-$	1.09	2.67	1.09	2.67	1.33	2.05	1.13	1.02
	AJIVE(T)	$-$	$-$	4.58	1.00	4.58	1.00	1.48	$-$	$-$	$-$
	SLIDE	$-$	$-$	0.97	1.02	0.97	1.02	0.84	$-$	$-$	$-$
	SVD(T)	$-$	$-$	1.00	3.61	1.00	3.61	1.20	1.87	1.02	1.02
	SVD(soft)	$-$	$-$	1.00	0.89	1.00	0.89	0.85	0.90	1.00	1.00
1	BIDIFAC	$-$	$-$	0.35	0.41	0.35	0.41	0.36	0.43	0.80	1.00
	UNIFAC	$-$	$-$	0.35	0.41	0.35	0.41	0.36	0.43	0.79	1.00
	JIVE(T)	$-$	$-$	0.18	0.23	0.18	0.23	0.19	0.25	0.63	1.01
	JIVE(P)	$-$	$-$	0.47	0.54	0.47	0.54	0.23	0.40	0.78	1.01
	AJIVE(T)	$-$	$-$	1.98	1.00	1.98	1.00	0.30	$-$	$-$	$-$
	SLIDE	$-$	$-$	0.19	0.22	0.19	0.22	0.19	$-$	$-$	$-$
	SVD(T)	$-$	$-$	1.00	1.73	1.00	1.73	0.22	0.30	1.01	1.01
	SVD(soft)	$-$	$-$	1.00	0.66	1.00	0.66	0.40	0.49	1.00	1.00
2	BIDIFAC	$-$	$-$	0.11	0.13	0.11	0.13	0.11	0.16	0.65	1.00
	UNIFAC	$-$	$-$	0.11	0.13	0.11	0.13	0.11	0.16	0.65	1.00
	JIVE(T)	$-$	$-$	0.04	0.05	0.04	0.05	0.04	0.05	0.54	1.01
	JIVE(P)	$-$	$-$	0.35	0.32	0.35	0.32	0.05	0.14	0.81	1.01
	AJIVE(T)	$-$	$-$	0.07	0.08	0.07	0.08	0.04	$-$	$-$	$-$
	SLIDE	$-$	$-$	0.05	0.06	0.05	0.06	0.06	$-$	$-$	$-$
	SVD(T)	$-$	$-$	1.00	1.37	1.00	1.37	0.05	0.06	1.01	1.01
	SVD(soft)	$-$	$-$	1.00	0.74	1.00	0.74	0.13	0.18	1.00	1.00
Mixed	BIDIFAC	$-$	$-$	0.37	0.41	0.37	0.41	0.31	0.38	0.80	1.00
	UNIFAC	$-$	$-$	0.37	0.41	0.37	0.41	0.31	0.36	0.77	1.00
	JIVE(T)	$-$	$-$	0.20	0.27	0.20	0.27	0.22	0.29	0.64	1.01
	JIVE(P)	$-$	$-$	0.55	0.73	0.55	0.73	0.29	0.47	0.83	1.01
	AJIVE(T)	$-$	$-$	2.28	0.92	2.28	0.92	0.64	$-$	$-$	$-$
	SLIDE	$-$	$-$	0.51	0.53	0.51	0.53	0.38	$-$	$-$	$-$
	SVD(T)	$-$	$-$	1.00	1.78	1.00	1.78	0.74	0.35	1.01	1.01
	SVD(soft)	$-$	$-$	1.00	0.70	1.00	0.70	0.59	0.43	1.00	1.00

Table 2. Table 2: Summary of the simulation studies for Design 2.

		Prediction							Imputation
SNR	Model	$G$	$R$	$C$	$I$	$G+C$	$R+I$	$S$	Cell	Column	Row
0.5	BIDIFAC	0.67	0.60	0.81	0.89	0.73	0.83	0.76	0.81	0.91	0.91
	UNIFAC	1.00	1.00	0.85	0.91	0.80	0.87	0.82	0.86	0.95	1.00
	JIVE(T)	1.00	1.00	3.05	3.86	1.00	1.38	1.04	1.47	1.14	1.02
	JIVE(P)	1.00	1.00	2.60	6.58	1.12	2.51	1.36	2.13	1.17	1.02
	AJIVE(T)	1.00	1.00	10.74	1.00	4.05	1.00	1.47	$-$	$-$	$-$
	SLIDE	1.00	1.00	1.00	1.27	0.96	0.97	0.85	$-$	$-$	$-$
	SVD(T)	1.00	1.00	1.00	8.51	1.00	3.32	1.19	1.86	1.02	1.02
	SVD(soft)	1.00	1.00	1.00	0.95	1.00	0.89	0.86	0.90	1.00	1.00
1	BIDIFAC	0.26	0.26	0.34	0.42	0.30	0.37	0.32	0.39	0.77	0.77
	UNIFAC	1.00	1.00	0.65	0.68	0.34	0.41	0.36	0.44	0.79	1.00
	JIVE(T)	1.00	1.00	1.62	1.76	0.18	0.23	0.19	0.25	0.63	1.01
	JIVE(P)	1.00	1.00	1.72	2.18	0.45	0.48	0.23	0.40	0.77	1.01
	AJIVE(T)	1.00	1.00	5.77	1.00	1.76	1.00	0.30	$-$	$-$	$-$
	SLIDE	1.00	1.00	0.99	0.93	0.19	0.22	0.19	$-$	$-$	$-$
	SVD(T)	1.00	1.00	1.00	4.82	1.00	1.56	0.22	0.29	1.01	1.01
	SVD(soft)	1.00	1.00	1.00	1.19	1.00	0.63	0.40	0.49	1.00	1.00
2	BIDIFAC	0.09	0.09	0.11	0.14	0.10	0.12	0.10	0.14	0.63	0.64
	UNIFAC	1.00	1.00	0.76	0.76	0.11	0.13	0.11	0.16	0.65	1.00
	JIVE(T)	1.00	1.00	1.33	1.39	0.04	0.05	0.04	0.06	0.54	1.01
	JIVE(P)	1.00	1.00	1.42	1.70	0.29	0.25	0.05	0.14	0.78	1.01
	AJIVE(T)	1.00	1.00	1.39	1.57	0.05	0.06	0.04	$-$	$-$	$-$
	SLIDE	1.00	1.00	1.02	1.19	0.05	0.06	0.06	$-$	$-$	$-$
	SVD(T)	1.00	1.00	1.00	4.07	1.00	1.22	0.05	0.06	1.01	1.01
	SVD(soft)	1.00	1.00	1.00	2.00	1.00	0.67	0.13	0.18	1.00	1.00
Mixed	BIDIFAC	0.30	0.28	0.37	0.48	0.31	0.35	0.27	0.33	0.77	0.78
	UNIFAC	1.00	1.00	0.69	0.79	0.36	0.40	0.31	0.37	0.77	1.00
	JIVE(T)	1.00	1.00	1.65	1.83	0.20	0.27	0.22	0.29	0.64	1.01
	JIVE(P)	1.00	1.00	1.85	2.53	0.53	0.65	0.30	0.46	0.87	1.01
	AJIVE(T)	1.00	1.00	6.35	1.11	2.07	0.90	0.60	$-$	$-$	$-$
	SLIDE	1.00	1.00	1.01	1.37	0.49	0.49	0.34	$-$	$-$	$-$
	SVD(T)	1.00	1.00	1.00	4.94	1.00	1.61	0.25	0.36	1.01	1.01
	SVD(soft)	1.00	1.00	1.00	1.39	1.00	0.66	0.36	0.43	1.00	1.00

Table 3. Table 3: Proportion of variance explained (Equation ( 10 )) for each component, estimated by BIDIFAC using the breast cancer data. The parentheses denote the rank of the estimated components.

	Global	Global+Row	Global+Col	Global+Row+Col	Signal
Tumor mRNA	0.14 (34)	0.32 (68)	0.45 (93)	0.58 (127)	0.67 (173)
NAT mRNA	0.23 (34)	0.50 (68)	0.44 (41)	0.66 (75)	0.78 (83)
Tumor miRNA	0.09 (34)	0.46 (67)	0.30 (93)	0.63 (126)	0.76 (175)
NAT miRNA	0.13 (34)	0.66 (67)	0.24 (41)	0.75 (74)	0.76 (79)

Table 4. Table 5: Mean relative recovery error for joint and individual signals, with standard error in parentheses, across different scenarios.

	$d=100,n=100$	$d=500,n=100$	$d=100,n=500$	$d=500,n=500$
UNIFAC	0.058 (0.002)	0.033 (0.001)	0.027 (0.001)	0.010 (0.001)
JIVE	0.106 (0.003)	0.110 (0.003)	0.023 (0.001)	0.022 (0.001)
AJIVE	0.096 (0.003)	0.100 (0.003)	0.021 (0.001)	0.021 (0.001)
SLIDE	0.132 (0.005)	0.130 (0.004)	0.026 (0.001)	0.024 (0.001)

Equations

$$\displaystyle{\bf A}_{00}=\left[\begin{array}[]{ccc}{\bf A}_{11}&\dots&{\bf A}_{1q}\\ \vdots&\ddots&\vdots\\ {\bf A}_{p1}&\dots&{\bf A}_{pq}\end{array}\right],$$

$${\mathbf{S}}_{ij}={\mathbf{G}}_{ij}+{\mathbf{R}}_{ij}+{\mathbf{C}}_{ij}+{\mathbf{I}}_{ij},$$

$${\mathbf{G}}_{00}={\mathbf{U}}_{00}^{(G)}{\mathbf{V}}_{00}^{(G)T},{\mathbf{R}}_{i0}={\mathbf{U}}_{i0}^{(R)}{\mathbf{V}}_{i0}^{(R)T},{\mathbf{C}}_{0j}={\mathbf{U}}_{0j}^{(C)}{\mathbf{V}}_{0j}^{(C)T},{\mathbf{I}}_{ij}={\mathbf{U}}_{ij}^{(I)}{\mathbf{V}}_{ij}^{(I)T}.$$

$$f_1(\{{\mathbf{U}}_{ij}^{(\cdot )},{\mathbf{V}}_{ij}^{(\cdot )}|i=0,\dots ,p,j=0,\dots ,q,(i,j)\ne (0,0)\})$$

$$=$$

$$+$$

$$\underset{\mathbf{U}:n\times r,\mathbf{V}:m\times r}{\min }||\mathbf{X}-\mathbf{U}{\mathbf{V}}^T|{|}_F^2+\lambda (||\mathbf{U}|{|}_F^2+||\mathbf{V}|{|}_F^2)=\underset{\mathbf{Z}:r(\mathbf{Z})\le r}{\min }||\mathbf{X}-\mathbf{Z}|{|}_F^2+2\lambda ||\mathbf{Z}|{|}_{\ast }$$

$$\displaystyle f_{2}(\mbox{\boldmath$\Theta$})=\frac{1}{2}\sum_{i=1}^{p}\sum_{j=1}^{q}||{\bf X}_{ij}-{\bf G}_{ij}-{\bf R}_{ij}-{\bf C}_{ij}-{\bf I}_{ij}||_{F}^{2}$$

$$+$$

$$\widehat{\mbox{\boldmath$\Theta$}}_{1}=\{\widehat{{\bf U}}_{i0}^{(G)}\widehat{{\bf V}}_{0j}^{(G)T},\widehat{{\bf U}}_{i0}^{(R)}\widehat{{\bf V}}_{ij}^{(R)T},\widehat{{\bf U}}_{ij}^{(C)}\widehat{{\bf V}}_{0j}^{({\bf C})T},\widehat{{\bf U}}_{ij}^{(I)}\widehat{{\bf V}}_{ij}^{(I)T}\mid i=1\ldots p,j=1,\ldots q\}$$

$$\widehat{\mbox{\boldmath$\Theta$}}_{2}=\{\widehat{{\bf G}}_{ij},\widehat{{\bf R}}_{ij},\widehat{{\bf C}}_{ij},\widehat{{\bf I}}_{ij}\mid i=1\ldots p,j=1,\ldots q\}$$

$$R_{{\mathbf{X}}_{ij}}^2({\mathbf{G}}_{ij})=1-\frac{||{\mathbf{X}}_{ij}-{\mathbf{G}}_{ij}|{|}_F^2}{||{\mathbf{X}}_{ij}|{|}_F^2}$$

$${\lambda }_{ij}=\sqrt{m_i}+\sqrt{n_j},\text{mbox}wherei=0,\dots ,p,\text{mbox}andj=0,\dots ,q.$$

$$\text{mbox}PredErr(\hat{\mathbf{G}})=\frac{{\sum }_{i=1,j=1}^{p,q}||{\mathbf{G}}_{ij}-{\hat{\mathbf{G}}}_{ij}|{|}_F^2}{{\sum }_{i=1,j=1}^{p,q}||{\mathbf{G}}_{ij}|{|}_F^2}.$$

$$\text{mbox}ImputeErr=\frac{{\sum }_{i=1,j=1}^{p,q}{\sum }_{(r,s)\in {\mathcal{I}}_{ij}}|{\mathbf{S}}_{ij}[r,s]-{\hat{\mathbf{S}}}_{ij}[r,s]{|}^2}{{\sum }_{i=1,j=1}^{p,q}{\sum }_{(r,s)\in {\mathcal{I}}_{ij}}|{\mathbf{S}}_{ij}[r,s]{|}^2}.$$

$$\widehat{\mbox{\boldmath$\Theta$}}_{1}=\{\widehat{{\bf U}}_{i0}^{(G)}\widehat{{\bf V}}_{0j}^{(G)T},\widehat{{\bf U}}_{i0}^{(R)}\widehat{{\bf V}}_{ij}^{(R)T},\widehat{{\bf U}}_{ij}^{(C)}\widehat{{\bf V}}_{0j}^{({\bf C})T},\widehat{{\bf U}}_{ij}^{(I)}\widehat{{\bf V}}_{ij}^{(I)T}\mid i=1\ldots p,j=1,\ldots q\}$$

$$\widehat{\mbox{\boldmath$\Theta$}}_{2}=\{\widehat{{\bf G}}_{ij},\widehat{{\bf R}}_{ij},\widehat{{\bf C}}_{ij},\widehat{{\bf I}}_{ij}\mid i=1\ldots p,j=1,\ldots q\}$$

$$||{\hat{\mathbf{U}}}_{ij}^{(I)}|{|}_F^2+||{\hat{\mathbf{V}}}_{ij}^{(I)}|{|}_F^2=\underset{\begin{array}{c}\mathbf{U}{\mathbf{V}}^T={\hat{\mathbf{I}}}_{ij}\\ \mathbf{U}:m_1\times \min (m_1,n_1)\\ \mathbf{V}:n_1\times \min (m_1,n_1)\end{array}}{\min }(||\mathbf{U}|{|}_F^2+||\mathbf{V}|{|}_F^2).$$

$$||{\hat{\mathbf{U}}}_{ij}^{(I)}|{|}_F^2+||{\hat{\mathbf{V}}}_{ij}^{(I)}|{|}_F^2=2||{\hat{\mathbf{I}}}_{ij}|{|}_{\ast }\text{ for all }i>0,j>0.$$

$$||{\hat{\mathbf{U}}}_{i0}^{(R)}|{|}_F^2+||{\hat{\mathbf{V}}}_{i0}^{(R)}|{|}_F^2$$

$$||{\hat{\mathbf{U}}}_{0j}^{(C)}|{|}_F^2+||{\hat{\mathbf{V}}}_{0j}^{(C)}|{|}_F^2$$

$$||{\hat{\mathbf{U}}}_{00}^{(G)}|{|}_F^2+||{\hat{\mathbf{V}}}_{00}^{(G)}|{|}_F^2$$

$$\widetilde{\mbox{\boldmath$\Theta$}}_{1}=\{\widetilde{{\bf U}}_{i0}^{(G)}\widetilde{{\bf V}}_{0j}^{(G)T},\widetilde{{\bf U}}_{i0}^{(R)}\widetilde{{\bf V}}_{ij}^{(R)T},\widetilde{{\bf U}}_{ij}^{(C)}\widetilde{{\bf V}}_{0j}^{(C)T},\widetilde{{\bf U}}_{ij}^{(I)}\widetilde{{\bf V}}_{ij}^{(I)T}\mid i=1\ldots p,j=1,\ldots q\}$$

$$2f_{2}(\widetilde{\mbox{\boldmath$\Theta$}}_{2})=f_{1}(\widetilde{\mbox{\boldmath$\Theta$}}_{1})\geq f_{1}(\widehat{\mbox{\boldmath$\Theta$}}_{1})=2f_{2}(\widehat{\mbox{\boldmath$\Theta$}}_{2}),$$

$$\displaystyle f_{2}\left(\alpha\widetilde{\mbox{\boldmath$\Theta$}}^{(1)}+(1-\alpha)\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}\right)\leq\alpha f_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(1)})+(1-\alpha)f_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}),$$

$$\{\alpha {\mathbf{G}}_{ij}^{(1)}+(1-\alpha ){\mathbf{G}}_{ij}^{(2)},\alpha {\mathbf{R}}_{ij}^{(1)}+(1-\alpha ){\mathbf{R}}_{ij}^{(2)},\alpha {\mathbf{C}}_{ij}^{(1)}+(1-\alpha )$$

$$\mid i=1,\dots ,p,j=1,\dots q\}.$$

$$\displaystyle f_{2}^{\text{LS}}(\mbox{\boldmath$\Theta$})$$

$$\displaystyle f_{2}^{\text{PEN}}(\mbox{\boldmath$\Theta$})$$

$$||{\mathbf{X}}_{ij}-(\alpha {\mathbf{A}}_{ij}^{(1)}+(1-\alpha ){\mathbf{A}}_{ij}^{(2)})|{|}_F^2\le \alpha ||{\mathbf{X}}_{ij}-{\mathbf{A}}_{ij}^{(1)}|{|}_F^2+(1-\alpha )||{\mathbf{X}}_{ij}-{\mathbf{A}}_{ij}^{(2)}|{|}_F^2$$

$$\displaystyle f_{2}^{\text{LS}}\left(\alpha\widetilde{\mbox{\boldmath$\Theta$}}^{(1)}+(1-\alpha)\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}\right)\leq\alpha f^{\text{LS}}_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(1)})+(1-\alpha)f^{\text{LS}}_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}).$$

$$||\alpha {\mathbf{A}}^{(1)}+(1-\alpha ){\mathbf{A}}^{(2)}|{|}_{\ast }\le \alpha ||{\mathbf{A}}^{(1)}|{|}_{\ast }+(1-\alpha )||{\mathbf{A}}^{(2)}|{|}_{\ast }$$

$$\displaystyle f_{2}^{\text{PEN}}\left(\alpha\widetilde{\mbox{\boldmath$\Theta$}}^{(1)}+(1-\alpha)\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}\right)\leq\alpha f^{\text{PEN}}_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(1)})+(1-\alpha)f^{\text{PEN}}_{2}(\widetilde{\mbox{\boldmath$\Theta$}}^{(2)}).$$

Full text

Keywords: BIDIFAC; Bidimensional data; Data integration; Dimension reduction; Principal component analysis

1 Introduction

1.1 Overview

Several recent methodological developments have been motivated by the integration of multiple sources of genetic, genomic, epigenomic, transcriptomic, proteomic, and other omics data. Successful integration of these disparate but related sources is essential for a complete understanding of the molecular underpinnings of human diseases, by providing essential tools for novel hypotheses (Hawkins et al., 2010) and improving statistical power. For example, TWAS and MetaXcan/PrediXcan have shown improved power for gene-based association testing, by integrating expression quantitative trait loci (eQTL) and genome-wide association (GWAS) data (Gamazon et al., 2015; Gusev et al., 2016). Similarly, the iBAG approach (Wang et al., 2012) has shown improved power for gene-based association by integrating messenger RNA (mRNA) levels with epigenomic data (e.g., DNA methylation). Moreover, recent studies showed that protein-level variations explain additional individual phenotypic differences not explained by the mRNA levels (Wu et al., 2013).

In addition to integrating multiple sources of high-dimensional data, integrating high-dimensional data across multiple patient cohorts can also improve interpretation and statistical power. For example, the integration of genome-wide data from multiple types of cancers can improve classification of oncogenes or tumor suppressors (Kumar et al., 2015) and may improve clinical prognoses (Liu et al., 2018).

Most statistical methods for the integration of high-dimensional matrices apply to data that are linked vertically (e.g., one cohort measured with more than one platforms, such as mRNA and miRNA) or horizontally (e.g., mRNA expression measured for multiple cohorts) (Tseng, Ghosh and Feingold, 2012). However, linked structures in molecular biomedical data are often more complex. In particular, the integration of bidimensionally linked data (e.g., more than one heterogeneous groups of subjects measured by more than one platform) is largely unaddressed. In this paper, we propose a new statistical method for the low-rank structural factorization of large bidimensionally linked datasets. This can be used to accomplish three important tasks: (i) missing value imputation (ii) dimension reduction, and (iii) the interpretation of lower-dimensional patterns that are shared across matrices or unique to particular matrices.

1.2 Motivating Example

The Cancer Genome Atlas (TCGA) is the most comprehensive and well-curated study of the cancer genome, with data for 6 different omics data sources from 11,000 patients representing 33 different cancer tumor types as well as rich clinical phenotypes. We consider integrating a cohort of breast cancer (BRCA) tumor samples, and a cohort of normal adjacent tissue (NAT) samples, from TCGA. NAT samples are often used for differential analyses, e.g., to identify genes with mean differential expression between cancer and normal tissue. However, such analyses do not address the molecular heterogeneity or trans-omic interactions that characterize cancer cells. Noticeably, Aran et al. (2017) conducted a comprehensive study on NAT across different cancers using TCGA and the Genotype-Tissue Expression (GTEx) program data and showed that the expression levels of NAT from breast, colon, liver, lung, and uterine tumors yield different clustering from their respective tumor tissues. In addition, Huang et al. (2016), using TCGA data, suggested that NATs not only serve as controls to tumor tissues but also provide useful information on patients’ survival that tumor samples do not. More detailed investigations of the molecular heterogeneity between tumor and NAT tissue are limited by available statistical methods, especially for multi-omics data. Our premise is that comprehensive analysis of multiple omics data sources across both tumor tissues and NAT would distinguish the joint signals that are shared across different omics profiles (e.g., mRNA and miRNA) and those that are only attributed to the tumors.

1.3 Existing Methods on Joint Matrix Factorization

Principal components analysis (PCA) and related techniques such as the singular value decomposition (SVD) are popular for the dimension reduction of a single data matrix ${\bf X}:m\times n$, resulting in the low-rank approximation ${\bf X}\approx{\bf U}{\bf V}^{T}$. Here, ${\bf U}$ are row loadings and ${\bf V}$ are column scores that together explain variation in ${\bf X}$. There is also a growing literature on the simultaneous dimension reduction of multiple data matrices ${\bf X}_{i}$ with size $m_{i}\times n$, which estimate low-rank signals that are jointly shared across data matrices. To capture joint variation, concatenated PCA assumes ${\bf X}_{i}={\bf U}_{i}{\bf V}^{T}$ for each matrix ${\bf X}_{i}$, i.e., the scores are shared across matrices. The iCluster (Shen, Olshen, and Ladanyi, 2009) and irPCA (Liu et al., 2016) approaches make this assumption for the integration of multi-source biomedical data. Alternatively, more flexible approaches allow for structured variation that may be shared across matrices or specific to individual matrices. The Joint and Individual Variations Explained (JIVE) method (Lock et al., 2013) decomposes joint and individual low-rank signals across matrices via the decomposition ${\bf X}_{i}={\bf U}_{i}{\bf V}^{T}+{\bf W}_{i}{\bf V}_{i}^{T}+{\bf E}_{i}$. In the context of vertical integration, the joint and individual scores ${\bf V}$ and ${\bf V}_{j}$ have been applied to risk prediction (Kaplan and Lock, 2017) and clustering (Hellton and Thoresen, 2016) for high-dimensional data. Several related techniques such as AJIVE (Feng et al., 2018) and SLIDE (Gaynanova and Li, 2017) have been proposed (Zhou et al., 2016), as well as extensions that allow the adjustment of covariates (Li and Jung, 2017) or accommodate heterogeneity in the distributional assumptions for different sources (Li and Gaynanova, 2018; Zhu et al., 2018).

The aforementioned methods focus exclusively on data that share a single dimension (i.e., either horizontally or vertically), and extension to matrices that are linked both vertically and horizontally is not straightforward. O’Connell and Lock (2019) decompose shared and individual low-rank structure for three interlinked matrices ${\bf X},{\bf Y},{\bf Z}$ where ${\bf X}$ and ${\bf Y}$ are shared vertically and ${\bf X}$ and ${\bf Z}$ are shared horizontally. However, their approach is not directly applicable to more general forms of bidimensionally linked data, and it suffers from potential convergence to a local minimum of the objective during estimation.

1.4 Our Contribution

We propose the first unified framework to decompose bidimensionally linked matrices into globally shared, horizontally shared (i.e., row-shared), vertically shared (i.e., column-shared), or individual structural components. Our specific aims are to (i) separate shared and individual structures, (ii) separate the shared components into one of globally-shared, column-shared, or row-shared structures, and (iii) maintain the low-rank structures for the signals. Our approach extends soft singular value thresholding (SSVT), i.e., nuclear norm penalization, for a single matrix. It requires optimizing a single convex objective function, which is relatively computationally efficient. It also facilitates a simple and intuitive approach based on random matrix theory for model specification, rather than complex and computationally expensive procedures to select tuning parameters or model ranks. Although our primary focus is bidimensional integration, our approach includes a novel method for vertical-only or horizontal-only integration as a special case. We show in simulation studies that our method outperforms existing methods, including JIVE.

The rest of this article is organized as follows. Section 2 describes the proposed method, denoted by BIDIFAC (bidimensional integrative factorization) for bidimensionally linked matrices, and addresses estimation, tuning parameter selection and imputation algorithms. Section 3 constructs simulated data under various scenarios and compares our method to existing methods in terms of structural reconstruction error and imputation performance. In Section 4, we apply BIDIFAC to the breast cancer and NAT data obtained from TCGA and illustrate the utility of the model. We conclude with some points of discussions in Section 5.

2 Methods

2.1 Notations and Definitions

Consider a set of $pq$ matrices $\{{\bf A}_{ij}:m_{i}\times n_{j}\mid i=1,\dots,p,\;j=1,\dots,q\}$, which may be concatenated to form the matrix

[TABLE]

where ${\bf A}_{00}:m_{0}\times n_{0}$ with with $m_{0}=\sum_{i=1}^{p}m_{i}$ and $n_{0}=\sum_{j=1}^{q}n_{j}$. Analogously, we define the column-concatenated matrices ${\bf A}_{i0}=[{\bf A}_{i1},\dots,{\bf A}_{iq}]$ for $i=1,\ldots,p$ and the row-concatenated matrices ${\bf A}_{0j}=[{\bf A}_{1j}^{T},\dots,{\bf A}_{pj}^{T}]^{T}$ for $j=1,\ldots,q$. We first define terms to characterize the relationships among these data matrices.

Definition 1**.**

In the arrangement (4), a set of data matrices with the structure of $\{{\bf X}_{ij}|i=1,\cdots,p,j=1,\cdots,q\}$ follows a $p\times q$ bidimensionally linked structure. The elements of $\{{\bf X}_{ij}|j=1,\cdots,q\}$ are row-shared and the elements of $\{{\bf X}_{ij}|i=1,\cdots,p\}$ are column-shared. The elements of $\{{\bf X}_{ij}|i=1,\cdots,p,j=1\cdots,q\}$ are globally-shared if every ${\bf X}_{ij}$ in the set is column-shared with ${\bf X}_{i^{\prime}j}$ and row-shared with ${\bf X}_{ij^{\prime}}$ for $i^{\prime}=1,\cdots,p$ and $j^{\prime}=1,\cdots,q$.

2.2 Model Specification

We assume that each matrix is decomposed by ${\bf X}_{ij}={\bf S}_{ij}+{\bf E}_{ij}$, where ${\bf S}_{ij}$ is a low-rank signal matrix and ${\bf E}_{ij}$ is a full-rank white noise. We further assume that ${\bf S}_{ij}$ can be decomposed as

[TABLE]

where ${\bf G}_{ij}$ is globally shared structure, ${\bf R}_{ij}$ is row-shared structure, ${\bf C}_{ij}$ is column-shared structure, and ${\bf I}_{ij}$ is individual structure for matrix ${\bf X}_{ij}$. The shared nature of the terms are apparent from their factorized forms. Defining the parameter set as \mbox{\boldmath\Theta}=\{{\bf G}_{00},{\bf R}_{i0},{\bf C}_{0j},{\bf I}_{ij}|i=1,\cdots,p,j=1,\cdots,q\}, we write each term as a product of row loadings ${\bf U}$ and column scores ${\bf V}$:

[TABLE]

The loadings and scores are shared across matrices for the global components ${\bf G}$, i.e., row-shared matrices have common global loadings ${\bf U}_{i0}^{(G)}$, and column-shared matrices have common global scores ${\bf V}_{0j}^{(G)}$. The row-shared structures ${\bf R}$ have common loadings, and the column-shared structures ${\bf C}$ have common scores. The dimensions of ${\bf U}$ and ${\bf V}$ depend on the global shared rank $r_{00}$, column-shared ranks $r_{0j}$, row-shared ranks $r_{i0}$, and individual ranks $r_{ij}$: ${\bf U}_{ij}^{(\cdot)}~{}:~{}m_{i}\times r_{ij}$ and ${\bf V}_{ij}^{(\cdot)}~{}:~{}n_{j}\times r_{ij}$ for $i=0,\cdots,p$, and $j=0,\cdots q$.

A diagram with the proposed notation for $2\times 2$ linked structure is shown in Figure 1.

To understand the factorized forms in (6), it is instructive to consider their interpretation under the motivating example of Section 1.2. Say ${\bf X}_{11}$ gives gene (mRNA) expression for tumor samples, ${\bf X}_{12}$ gives gene expression for NAT samples, ${\bf X}_{21}$ gives miRNA expression for tumor samples, and ${\bf X}_{22}$ gives miRNA expression for NAT samples. Then, the loadings of the global structures ${\bf G}_{00}$ give trans-omic signatures ${\bf U}_{00}^{(G)}$ that explain substantial variability across both tumor and NAT samples with associated scores ${\bf V}_{00}^{(G)}$. The loadings of the column-shared structures ${\bf C}_{01}$ include trans-omic signatures ${\bf U}_{01}^{(C)}$ that explain substantial variability in the tumor samples with associated scores ${\bf V}_{01}^{(C)}$ but not the NAT samples. The loadings of the row-shared structures ${\bf R}_{10}$ include gene signatures ${\bf U}_{10}^{(R)}$ that explain substantial variability across both tumor and NAT samples, but are unrelated to miRNA. The loadings of the individual structures ${\bf I}$ include gene signatures unrelated to miRNA that explain variability in only the tumor samples, ${\bf U}_{11}^{(I)}$.

Marginally, ${\bf G}_{i0}+{\bf R}_{i0}$ denotes the components shared by ${\bf X}_{i0}$ and ${\bf G}_{0j}+{\bf C}_{0j}$ denotes the components shared by ${\bf X}_{0j}$. In practice, we assume the errors ${\bf E}_{ij}$ are Gaussian white noise with mean zero and variance $\sigma_{ij}^{2}$. For vertical integration, $q=1$ and the ${\bf G}_{ij}$ and ${\bf R}_{ij}$ terms are redundant with ${\bf C}_{ij}$ and ${\bf I}_{ij}$, respectively. Thus, we suppress these terms and the decomposition of the model reduces to the “joint and individual” structures as in one-dimensional factorization methods, including JIVE. Note, however, that we do not require pairwise orthogonality constraints on ${\bf G}_{ij}$, ${\bf R}_{ij}$, ${\bf C}_{ij}$, and ${\bf I}_{ij}$, which will be discussed in the remaining sections.

2.3 Estimation

Without any penalization, estimation of $\Theta$ would have an identifiability issue. We first consider minimizing the sum of squared errors over all matrices with different levels of matrix $L_{2}$ penalties on the $\{{\bf U}_{ij}^{(\cdot)},{\bf V}_{ij}^{(\cdot)}\}$. Our objective function is

[TABLE]

where $||\cdot||_{F}$ denotes the Frobenious norm and each $\lambda_{ij}$ is a non-negative penalty factor. We upper bound the ranks by setting $r_{ij}=\min(m_{i},n_{j})$ for $i=0,\ldots,p$ and $j=0,\ldots,q$; the actual ranks of the solution may be lower, as discussed below. The objective function (7) is a convex function of each of ${\bf U}_{ij}^{(\cdot)}$ and ${\bf V}_{ij}^{(\cdot)}$, given all the others fixed. One may use an alternating least squares (ALS) with a matrix $L_{2}$ penalty to iteratively update each of $\{{\bf U}_{00}^{(G)},{\bf V}_{00}^{(G)}\}$, $\{{\bf U}_{i0}^{(R)},{\bf V}_{i0}^{(R)}\}$, $\{{\bf U}_{0j}^{(C)},{\bf V}_{0j}^{(C)}\}$, $\{{\bf U}_{ij}^{(I)},{\bf V}_{ij}^{(I)}\}$ until convergence.

Alternatively, we reformulate (7) and motivate our model using nuclear norm penalties. A matrix ${\bf A}:m\times n$ with ordered singular values $\delta_{1},\delta_{2},\ldots$ has nuclear norm $||{\bf A}||_{*}=\sum_{i=1}^{\text{min}\{m,n\}}\delta_{i}.$ We first present a well-known result on the equivalence of nuclear norm penalization and matrix factorization for a single matrix in Proposition 1.

Proposition 1**.**

(Mazumder, Hastie, and Tibshirani, 2010)* For a matrix ${\bf X}:m\times n$,*

[TABLE]

where $r=min(m,n)$. Moreover, $\widehat{{\bf Z}}=\widehat{{\bf U}}\widehat{{\bf V}}^{T}$, where $\widehat{{\bf Z}}$ solves the right-hand side of (8) and $\{\widehat{{\bf U}},\widehat{{\bf V}}\}$ solves the left-hand side of (8).

Proposition 1 depends on Lemma 1, which is also shown in Mazumder, Hastie, and Tibshirani (2010).

Lemma 1**.**

For a matrix ${\bf Z}:m\times n$, $\underset{\begin{subarray}{c}{\bf U}{\bf V}^{T}={\bf Z}\\ {\bf U}:m\times\min(m,n)\\ {\bf V}:n\times\min(m,n)\end{subarray}}{\min}(||{\bf U}||_{F}^{2}+||{\bf V}||_{F}^{2})=2||{\bf Z}||_{*}.$

We extend the nuclear norm objective in (8) to our context as follows:

[TABLE]

Theorem 1 establishes the equivalence of (7) and (9), with proof in Appendix B.

Theorem 1**.**

Let

[TABLE]

minimize (7). Then,

[TABLE]

minimizes (9), where $\widehat{{\bf G}}_{ij}=\widehat{{\bf U}}_{i0}^{(G)}\widehat{{\bf V}}_{0j}^{(G)T},\widehat{{\bf R}}_{ij}=\widehat{{\bf U}}_{i0}^{(R)}\widehat{{\bf V}}_{ij}^{(R)T},\widehat{{\bf C}}_{ij}=\widehat{{\bf U}}_{ij}^{(C)}\widehat{{\bf V}}_{0j}^{(C)T},$ and $\widehat{{\bf I}}_{ij}=\widehat{{\bf U}}_{ij}^{(I)}\widehat{{\bf V}}_{ij}^{(I)T}$.

The objective (9) has several advantages. First, the function is convex, which we state in Theorem 2 and prove in Appendix B.

Theorem 2**.**

The objective $f_{2}(\cdot)$ in (9) is convex over its domain.

Fortunately, minimizing one term (${\bf G}_{00},{\bf R}_{i0},{\bf C}_{0j}$, or ${\bf I}_{ij}$) with the others fixed is straightforward via soft singular value thresholding (SSVT). We state the well-known equivalence between nuclear norm penalization and SSVT in Proposition 2; for a proof see Mazumder, Hastie, and Tibshirani (2010).

Proposition 2**.**

If the SVD of ${\bf X}$ is ${\bf U}_{X}{\bf D}_{X}{\bf V}_{X}^{T}$ and ${\bf D}_{X}$ has diagonal entries $\delta_{1}\geq\dots\geq\delta_{r}\geq 0$, the solution for ${\bf Z}$ in (8) is equal to $\widehat{{\bf Z}}={\bf U}_{X}{\bf D}_{X}(\lambda){\bf V}_{X}^{T}$, where ${\bf D}_{X}(\lambda)$ is a diagonal matrix with $\delta_{1},\dots,\delta_{r}$ replaced by $\max(\delta_{1}-\lambda,0),\ldots,\max(\delta_{r}-\lambda,0)$, respectively.

We use an iterative soft singular value thresholding (ISSVT) algorithm to solve (9), applying Proposition 2 to the appropriate residual matrices for ${\bf G}_{00},{\bf R}_{i0},{\bf C}_{0j}$, or ${\bf I}_{ij}$. For example, $\widehat{{\bf I}}_{ij}$ and $\widehat{{\bf R}}_{i0}$ are obtained by soft-thresholding the singular values of ${\bf X}_{ij}-\widehat{{\bf G}}_{ij}-\widehat{{\bf R}}_{ij}-\widehat{{\bf C}}_{ij}$ and ${\bf X}_{i0}-\widehat{{\bf G}}_{i0}-\widehat{{\bf I}}_{i0}$ towards 0, respectively. This iterative algorithm is guaranteed to converge to a coordinatewise-minimum, and convexity implies that it will be a global minimum if it is a local minimum. In practice, we find that iterative algorithms to solve either (7) or (9) converge to the same solution and are robust to their initial values. The detailed algorithm of ISSVT is provided in Appendix A. The resulting global, column-shared, row-shared and individual terms of the decomposition will have reduced rank depending on the penalty factors $\lambda_{ij}$. Moreover, this relation between the penalty factors and the singular values motivates a straightforward choice of tuning parameters using random matrix theory described in Section 2.6.

ISSVT can also be represented as a blockwise coordinate descent algorithm for the ALS objective (7) (see Appendix A), and it converges faster than ALS, so we use it as our default algorithm. However, the ALS approach may be extended to certain related contexts. For example, it can incorporate sparsity in the loadings via an additional penalty (e.g., $L_{1}$) on ${\bf U}_{ij}^{(\cdot)}$. Also, with only three linked matrices ${\bf X}_{11},{\bf X}_{12}$ and ${\bf X}_{21}$, as in O’Connell and Lock (2019), (9) cannot properly construct the globally shared components. The formulation of (7) with an ALS algorithm can handle such contexts.

2.4 Data Pre-Processing

In practice, the data matrices ${\bf X}_{ij}$ may have very different levels of variability or be measured on different scales. Thus, a straightforward application of the objective (7) or (9) is not appropriate without further processing. By default we center the matrices to have mean [math], so that each matrix has the same baseline. To resolve issues of scale, we propose dividing each data matrix ${\bf X}_{ij},i,j>0$ by an estimate of the square root of its noise variance $\hat{\sigma}_{ij}$, denoted by ${\bf X}_{ij,scale}$. We discuss estimating the Gaussian noise variance in Section 2.6. After scaling, each matrix has homogeneous unit noise variance, which motivates the proposed penalties. After all components, denoted by ${\bf G}_{ij,scale},{\bf R}_{ij,scale},{\bf C}_{ij,scale},$ and ${\bf I}_{ij,scale}$, are estimated, we transform the results back to the original scale by multiplying each matrix by $\hat{\sigma}_{ij}$.

We comment on estimating the noise variance $\sigma_{ij}^{2}$. Without any signal or with a very weak signal, the standard deviation of $\mbox{vec}({\bf X}_{ij})$, denoted as $\hat{\sigma}_{ij}^{SD}$, provides a nearly unbiased estimate of $\sigma_{ij}$. However, this estimate is biased and overly conservative with a high signal-to-noise ratio. An alternative is to use random matrix theory, and estimate $\sigma_{ij}$ by minimizing the Kolmogorov-Smirnov distance between the theoretical and empirical distribution functions of singular values, as in Shabalin and Nobel (2013). Their estimate $\hat{\sigma}_{ij}^{KS}$ is based on grid-search on a candidate set of $\sigma$. Recently, Gavish and Donoho (2017) proposed another estimator $\hat{\sigma}_{ij}^{MAD}$, also based on random matrix theory, which is defined as the median of the singular values of ${\bf X}_{ij}$ divided by the square root of the median of the Marcenko-Pastur distribution. Our simulations, not shown here, revealed that both $\hat{\sigma}_{ij}^{KS}$ and $\hat{\sigma}_{ij}^{MAD}$ well approximate the standard deviation of a true noise matrix when the data matrix consists of low rank signal, and we use $\hat{\sigma}_{ij}^{MAD}$ as a default throughout this paper for its simplicity. From here, we assume that BIDIFAC is applied to the data matrices with $\sigma_{ij}^{2}=1$.

2.5 Summarizing Results

Given that the mean of each matrix is 0, we propose proportion of variance explained $(R_{{\bf X}_{ij}}^{2}(\cdot))$ as a summary statistic. For example,

[TABLE]

provides a measure of the proportion of variability explained by the globally shared component. However, because orthogonality is not explicitly enforced in BIDIFAC, this equality does not hold in general and $R^{2}(\cdot)$ is not necessarily additive across terms (e.g., $R_{{\bf X}_{ij}}^{2}({\bf G}_{ij}+{\bf C}_{ij})\neq R_{{\bf X}_{ij}}^{2}({\bf G}_{ij})+R_{{\bf X}_{ij}}^{2}({\bf C}_{ij})$).

2.6 Selecting Tuning Parameters

The performance of the proposed method depends heavily on the choice of tuning parameters. In the literature, there are several approaches to select ranks in the context of vertical integration, including permutation testing (Lock et al., 2013), BIC (O’Connell and Lock, 2016), and cross-validation (Li and Jung, 2017). In our context, the issue of rank selection is analogous to selecting the tuning parameters $\lambda_{ij}$. Although cross-validation is a natural way of selecting tuning parameters in penalized regression, our objective involves too many parameters ($1+p+q+pq$) to be computationally feasible. Moreover, despite the rich literature on cross-validating SVD for a single matrix, it is not clear how to define the training and test set (e.g., randomly select cells, rows, columns, a whole matrix, etc). A general description of the difficulties in cross-validating matrices is provided by Owen and Perry (2009).

Admitting that cross validation in BIDIFAC is not straightforward and inefficient, we provide an alternative approach to select the tuning parameters based on random matrix theory. We first construct necessary conditions for each element of \widehat{\mbox{\boldmath\Theta}} to be nonzero.

Proposition 3**.**

The following conditions are necessary to allow for non-zero $\widehat{{\bf G}}_{00}$, $\widehat{{\bf R}}_{i0}$, $\widehat{{\bf C}}_{0j}$, and $\widehat{{\bf I}}_{ij}$:

$\mbox{max}_{j}\lambda_{ij}<\lambda_{i0}<\sum_{j}\lambda_{ij}$* for $i=1,\ldots,p$ and $\mbox{max}_{i}\lambda_{ij}<\lambda_{0j}<\sum_{i}\lambda_{ij}$ for $j=1,\ldots,q$* 2. 2.

$\mbox{max}_{i}\lambda_{i0}<\lambda_{00}<\sum_{i}\lambda_{i0}$* and $\mbox{max}_{j}\lambda_{0j}<\lambda_{00}<\sum_{j}\lambda_{0j}$.*

We provide a proof of Proposition 3 in Appendix B. Without loss of generality, suppose that each cell of ${\bf X}_{ij}$ has Gaussian noise with unit variance ($\sigma_{ij}^{2}=1$), which is independent within each matrix and across matrices. Based on random matrix theory, we propose using the following penalty factors:

[TABLE]

It is straightforward to show that our choice of tuning parameters meets the necessary requirement. Also, under the aforementioned assumptions, $\sqrt{m_{i}}+\sqrt{n_{j}}$ provides a tight upper bound for the largest singular value of ${\bf E}_{ij}$ (Rudelson and Vershynin, 2010). Thus, without any shared structure, the motivation for (11) is apparent by considering the penalty as a soft-thresholding operator on the singular values in Proposition 2. The penalty (11) is also used in the single matrix reconstruction method in Shabalin and Nobel (2013). The penalties for the shared components are decided analogously because the stack of column-/row- shared matrices are also Gaussian random matrices with unit noise variance. For example, the penalty for ${\bf R}_{i0}$ is determined given by an estimate of the largest singular value of its concatenated noise matrix ${\bf E}_{i0}$ with unit variance: $\sqrt{m_{i}}+\sqrt{n_{0}}$.

2.7 Imputation

A convenient feature of BIDIFAC is its potential for missing value imputation. For single block data, PCA or related low-rank factorizations can be used to impute missing values by iteratively updating missing entries with their low-rank approximation, and this approach has proven to be very accurate in many applications (Kurucz, Benczúr and Csalogány, 2007). An analogous algorithm has been used for imputation with joint matrix factorization (O’Connell and Lock, 2019), and this approach readily extends to BIDIFAC. Importantly, this allows for the imputation of data that are missing an entire row or column within a block, via an expectation-maximization (EM) approach. Our imputation algorithm is presented below.

Let $\mathcal{I}_{ij}=\{(r,s)~{}|~{}{\bf X}_{ij}[r,s]\mbox{ is missing}\}$. For each matrix, initialize the missing values by the column- and/or row-wise mean and denote the initial matrix by $\widehat{{\bf X}}_{ij}^{(old)}$. 2. 2.

1. (a)

Maximization: Apply BIDIFAC to $\{\widehat{{\bf X}}_{ij}^{(old)}\}_{i,j=1}^{p,q}$. 2. (b)

Expectation: For $(r,s)\in\mathcal{I}_{ij}$, replace $\widehat{{\bf X}}_{ij}^{(old)}[r,s]$ with $\widehat{{\bf S}}_{ij}[r,s]$ and denote the imputed matrix by $\widehat{{\bf X}}_{ij}^{(new)}$. 3. 3.

If $\sum_{i,j=1}^{p,q}\sum_{(r,s)\in\mathcal{I}_{ij}}|\widehat{{\bf X}}_{ij}^{(old)}[r,s]-\widehat{{\bf X}}_{ij}^{(new)}[r,s]|^{2}<\epsilon$, the algorithm converges. If not, re-apply (2) after replacing $\widehat{{\bf X}}_{ij}^{(old)}$ by $\widehat{{\bf X}}_{ij}^{(new)}$.

To improve computational efficiency, each maximization step uses the \widehat{\mbox{\boldmath\Theta}} from the previous maximization step as starting values. The imputation scheme can be used to impute entries or entire rows or columns of the constituent data matrices; however, in all cases the missing entities (entries, rows or columns) must be missing at random.

Our algorithm can be considered a regularized EM algorithm, using a model-based motivation for the objective function (7). Because each component is estimated as a product of two matrices, it is naturally translated as a probabilistic matrix factorization. The unpenalized objective with $\lambda_{ij}=0$ for all $i,j$ maximizes a Gaussian likelihood model if the noise variances are the same across the matrices (i.e., $\sigma^{2}=\sigma_{11}^{2}=\dots=\sigma_{pq}^{2}$). With penalization, the approach is analogous to maximizing the posterior distribution in a Bayesian context with Gaussian prior on the terms of ${\bf U}_{ij}^{(\cdot)}$ and ${\bf V}_{ij}^{(\cdot)}$. Specifically, if the entries of ${\bf U}_{ij}^{(\cdot)}$ and ${\bf V}_{ij}^{(\cdot)}$ are independent Gaussian with mean [math] and variance $\sigma^{2}/\lambda_{ij}$, then minimizing (7) is analogous to finding the posterior mode (Mnih and Salakhutdinov, 2008). This provides a theoretical foundation for the iterative imputation algorithm based on expectation-maximization (EM) algorithm. It is closely related to softImpute, proposed by Mazumder, Hastie, and Tibshirani (2010), a computationally efficient imputation algorithm for a single matrix based on Proposition 1. Prediction intervals for the imputed values may be obtained via a resampling approach, as described in Web Appendix E of O’Connell and Lock (2019).

3 Simulation Studies

3.1 Simulation Setup

In this section, we compare our model to existing factorization methods using simulated data. Because competing approaches apply only to uni-dimensionally linked matrices (vertical or horizontal), we constructed two simulation designs with $p=q=2$ for proper comparison. Design 1 does not include any row-shared or globally shared structure and may be considered as two separate sets of vertically-linked matrices. Design 2 includes all linked structures that are represented in our model: global, column-shared, row-shared and individual.

For each simulation design, we used existing methods to compare performance. First, we fit 2 separate JIVE models to $\{{\bf X}_{11},{\bf X}_{21}\}$ and $\{{\bf X}_{12},{\bf X}_{22}\}$, using both (i) true marginal ranks for joint and individual components and (ii) rank selection based on permutation testing, denoted by JIVE(T) and JIVE(P) respectively. The “true marginal rank” in this context means $r({\bf G}_{ij}+{\bf C}_{ij})$ for the joint component and $r({\bf R}_{ij}+{\bf I}_{ij})$ for the individual components given that $j$ is fixed. We similarly apply the AJIVE and SLIDE methods, where we used the rank of ${\bf S}_{ij}$ as the initial rank for AJIVE. We also consider an approach that is analogous to BIDIFAC but reduced to vertical integration only, i.e., with $\widehat{{\bf G}}_{ij}$ and $\widehat{{\bf R}}_{ij}$ set to ${\bf 0}_{m_{i}\times n_{j}}$, denoted as UNIFAC in this paper. We fit UNIFAC to $\{{\bf X}_{11},{\bf X}_{21}\}$ and $\{{\bf X}_{12},{\bf X}_{22}\}$ separately. For notational simplicity, we denote the joint and individual components estimated by JIVE(P), JIVE(T), AJIVE, SLIDE and UNIFAC by $\widehat{{\bf C}}_{ij}$ and $\widehat{{\bf I}}_{ij}$. We also applied soft singular value thresholding for each single matrix with the corresponding imputation algorithm, softImpute, with tuning parameters decided as in Section 2.6 $(\lambda_{ij}=\hat{\sigma}_{ij}\cdot(\sqrt{m_{i}}+\sqrt{n_{j}}))$. We denote this approach by SVD(soft). We also consider the performance of the hard-thresholding low-rank approximation of each matrix, denoted by SVD(T), using the true marginal rank for a single matrix ($r({\bf S}_{ij})$). Similar to the above, SVD(soft) and SVD(T) estimate ${\bf I}_{ij}$ (or, equivalently, ${\bf S}_{ij}$) components only and assume $\widehat{{\bf G}}_{ij}=\widehat{{\bf R}}_{ij}=\widehat{{\bf C}}_{ij}={\bf 0}_{m_{i}\times n_{j}}$. BIDIFAC, UNIFAC, and SVD(soft) are soft-thresholding methods, while JIVE(T), JIVE(P), and SVD(T) are based on hard-thresholding.

In our simulation studies, the number of rows and columns for each matrix was set to $100$: $m_{i}=n_{j}=100$. The rank of the total signal in each matrix, ${\bf S}_{ij}={\bf G}_{ij}+{\bf R}_{ij}+{\bf C}_{ij}+{\bf I}_{ij}$, was $10$. This total rank was distributed across each of the $4$ terms (or the two terms ${\bf C}_{ij}$ and ${\bf I}_{ij}$ for Design 1) via a multinomial distribution with equal probabilities. For clarity of the simulation studies and to allow for comparison with other methods, we enforced orthogonality among the shared structures, both within each matrix (i.e., ${\bf G}_{11}$ and ${\bf C}_{11}$ are orthogonal) and across matrices (i.e., ${\bf G}_{11}$ and ${\bf G}_{12}$ are orthogonal). Note, however, that our model does not enforce orthogonality when estimating parameters.

Each signal matrix, ${\bf S}_{ij}$, was generated by applying SVD to a Gaussian random matrix with mean 0 and unit variance, denoted as ${\bf Y}={\bf U}_{Y}{\bf D}_{Y}{\bf V}_{Y}^{T}$. Then ${\bf U}_{Y}$ and ${\bf V}_{Y}$ were rearranged accordingly to guarantee orthogonality and a $2\times 2$ linked structure. The singular values were randomly permuted within each shared component to allow for heterogeneity in the size of the joint signal across matrices, e.g., ${\bf C}_{11}$ and ${\bf C}_{21}$ have the same loadings and scores but with different order. Each signal matrix ${\bf S}_{ij}$ was standardized to have $||{\bf S}_{ij}||_{F}=1$. Finally, independent Gaussian noise was added to each signal matrix, where the noise variance was decided by the signal-to-noise ratio (SNR) defined by $1/(\sigma_{ij}\cdot\sqrt{m_{i}\cdot n_{j}})$. We first considered three SNRs in our simulation studies: 0.5, 1, and 2. For the true structure ${\bf S}_{ij}$ the expected value of $R_{{\bf X}_{ij}}^{2}({\bf S}_{ij})$ is $\mbox{SNR}^{2}/(1+\mbox{SNR}^{2})$, when SNRs are the same for all matrices. We also considered a scenario where the SNR is randomly selected separately for the different matrices, uniformly between 0.5 and 2, to accommodate heterogeneous noise variances.

We compared the performance of our method and the competing methods from two perspectives: prediction error and imputation performance. We computed prediction error for each term in the decomposition (${\bf G},{\bf R},{\bf C}$ or ${\bf I}$) as the relative reconstruction error:

[TABLE]

For a fair comparison between our method and the existing methods for uni-dimensionally shared matrices, we also report $\mbox{PredErr}(\widehat{{\bf G}}+\widehat{{\bf C}})$, $\mbox{PredErr}(\widehat{{\bf R}}+\widehat{{\bf I}})$. Finally, we also report $\mbox{PredErr}(\widehat{{\bf S}})$ to evaluate the overall signal reconstruction performances.

For imputation, we considered three scenarios where in each matrix (i) 200 randomly selected cells are missing, (ii) 2 columns are missing, and (iii) 2 rows are missing. To foster borrowing information from the shared structures in (ii) and (iii), there was no overlapping row/column that are missing simultaneously in shared matrices. To reduce computation cost, imputation using JIVE(P) used fixed ranks determined by the complete data, resulting in slightly inflated imputation performances. We evaluated the imputation performance using the scaled reconstruction error for missing cells, defined by

[TABLE]

3.2 Results

We repeated each simulation 200 times and averaged the performance. The results are shown in Tables 1 and 2. We summarize the results below from a few perspectives.

On performance of BIDIFAC: BIDIFAC was competitive against all compared models in separating signals into multiple shared structures. In Design 1 where $||{\bf G}_{ij}||_{F}^{2}=||{\bf R}_{ij}||_{F}^{2}=0$ , the loss from the model misspecification was negligible when comparing the prediction errors for column-shared matrices (${\bf C}$) and the sum of global and column-shared matrices (${\bf G}+{\bf C}$), as well as individual matrices (${\bf I}$) and the sum of row and individual matrices (${\bf R}+{\bf I}$). The negligible effect of misspecification can also be seen when compared to UNIFAC, where the global or row-shared structures are ignored. In Design 2, where BIDIFAC is the only method estimating global (${\bf G}$) and row-shared (${\bf R}$) components, the prediction errors were comparable to the other components regardless of SNR. Even when error variances differ by each matrix, BIDIFAC was not affected severely.

On soft and hard thresholding: When SNR is low, matrix completion via soft thresholding generally outperformed hard-thesholding approaches (JIVE, SLIDE and SVD(T)). It is easily seen by comparing SVD(soft) to SVD(T) or comparing JIVE(T)/JIVE(P)/SLIDE to UNIFAC. Especially, when SNR$=0.5$ in Design 1, JIVE(T) overfitted severely even when true ranks are given. JIVE(P) performed even worse, due to erroneous rank selection. SLIDE suffered less from overfitting, but it was because it selected 0 rank in most simulated data. As SNR increases, SLIDE, JIVE(P) and JIVE(T) outperformed UNIFAC and BIDIFAC in estimating ${\bf G}+{\bf C}$ and ${\bf R}+{\bf I}$ in both designs. This result is intuitive: soft-thresholding prevents over-fitting when the SNR is low, but over-penalizes the estimated signal when SNR is high. We found that AJIVE did not perform well unless SNR is 2, where it provides provides similar results to JIVE with the true ranks, JIVE(T).

On overall signal recovery: In Design 1, JIVE(T) can be considered as the gold standard for overall signal recovery as it reflects the true joint and individual rank structures, which are unknown in practice. Except when SNR is 0.5 (which is explained by the difference between soft and hard thresholding), JIVE(T), JIVE(P) and SVD(soft) performed better than BIDIFAC in recovering the overall signals. The overall signal recovery of BIDIFAC was the same as UNIFAC and better than SVD(soft), revealing that BIDIFAC and UNIFAC obtained additional power from the column-shared matrices and the effect of model misspecification of BIDIFAC is negligible. In Design 2, BIDIFAC performed better than UNIFAC and SVD(soft) as it closely matched the data generating process.

On imputation performance: JIVE(T) was the winner in both designs in imputing missing cells and columns except when SNR is 0.5. However, recall that JIVE(T) uses the true ranks which are unknown in practice. BIDIFAC performed close to or even better than JIVE(P) in both designs, revealing that signal detection does not necessarily guarantee imputation performance as it is also affected by appropriate rank selection and detecting shared structures. In Design 1, it is not surprising that all models suffered when an entire row is missing, as there was no row-shared structure. In Design 2, BIDIFAC was the only method that successfully imputes missing rows.

To summarize, the performance of BIDIFAC was promising in simulation studies even when it was misspecified. It appropriately separated signals into linked structures and did not overfit for low SNR. Among the compared models, BIDIFAC is the only model that well accommodates the cases where a whole data matrix is missing or both a whole column and a whole row is missing. Acknowledging that it is not possible to obtain the true ranks of shared and individual structures, BIDIFAC performed the best across different scenarios.

4 Data Analysis

4.1 TCGA Breast Cancer Data

In this section, we apply our method to breast cancer data from TCGA(Cancer Genome Atlas Network, 2012). We integrate mRNA and miRNA profiles of the tumor samples and NATs, which, to our knowledge, has not been previously investigated in a fully unified framework. The data used here are freely available from TCGA. Specifically, we obtained the level III raw count of mRNA (RNASeq-V2) and miRNA (miRNA-Seq) data using the R package TCGA2STAT (Wan et al., 2015). In our analysis we first removed the tumor data of those samples with matched NAT, so that tumor and NAT data correspond to two independent cohorts. We also filtered mRNAs and miRNAs with more than half zero counts for all individuals. Then we took $\log(1+\mbox{count})$ and centered each of mRNA and miRNA profiles to the mean of tumors and NATs, so that the mean of each row of each matrix is 0. Lastly, we selected the 500 mRNAs and miRNAs with maximum variability. The resulting data had $500$ mRNA and miRNA profiles for $660$ tumor samples and $86$ NATs.

Breast cancer tumor samples are classified into 5 intrinsic subtypes, based on expression levels of 50 pre-defined genes: Luminal A (LumA), Luminal B (LumB), HER2-enriched (HER2), Basal-enriched (Basal), and Normal-like tumors (Ciriello et al., 2015). In our TCGA data, 419 out of 660 tumor samples had the labeled subtypes.

4.2 Results

We first applied BIDIFAC to the processed data until convergence. We summarized the proportion of explained variance, as well as the estimated ranks of the components, in Table 3. The difference between $R_{{\bf X}_{ij}}^{2}({\bf S}_{ij})-R_{{\bf X}_{ij}}^{2}({\bf G}_{ij}+{\bf R}_{ij}+{\bf C}_{ij})$ among tumors and NATs suggests that most of the variability of the miRNA and mRNA profiles were attributed to shared structures. In particular, more than a half of the variability of both the mRNA and miRNA profiles from the NAT is attributed to the row-shared components (Global+Row); this makes sense, as tumor cells are derived from normal tissue and thus we expect much of the same patterns of variability that are present in normal tissue to also be present in tumors. Between 30-45% of the variability is explained by the shared structure between miRNA and mRNA across the four matrices. There is a large difference of sample sizes between tumor and normal tissues, which may have also affected the estimate of proportion of variance explained. Even though the rank of the estimated signals in NAT was close to the rank of the data matrix (86) linear independence among the terms of the decomposition was preserved in our decomposition.

We compared BIDIFAC to existing methods on one-dimensionally linked matrices, including UNIFAC, JIVE(P) and SLIDE. We used two criteria to compare models: subtype classification and survival analysis. Specifically, we hypothesize that accounting for multiple –omics profiles and removing the shared variations between tumors and normal tissues would increase the biological interpretation of cancer subtypes and patient’s survival. Thus, we focused on column-shared and individual components estimated by BIDIFAC and compared it to other methods using tumor mRNA and miRNA only.

We first compared how the estimated components are well distinguished by breast cancer subtypes. To summarize the subtype distinctions, we used the SWISS (Standardized WithIn Class Sum of Squares) score (Cabanski et al., 2010). Interpreted similar to ANOVA, a lower SWISS score implies more clear subtype distinction. We restricted the attention to the tumor samples with labeled subtypes. As summarized in Table 4, the SWISS score of ${\bf C}+{\bf I}$ of UNIFAC was superior to JIVE and SLIDE, suggesting that soft-thresholding better uncovered subtype heterogeneity when integrating mRNA and miRNA. Compared to UNIFAC, BIDIFAC slightly obtained better SWISS score with lower rank of the ${\bf C}+{\bf I}$. The scatterplot of the principal components from the column-shared structures of BIDFAC, shown in Figure 2, reveals that the subtypes are well distinguished. The column-shared structure from JIVE(P) had the lowest SWISS score overall, which may be because of its extremely low rank.

Using overall survival data for patients with tumors, we applied the Cox proportional hazards model (PHM) using scores from the estimated components as predictors and used the score test to assess its significance. The results are also shown in Table 4. At $\alpha=0.05$, all models suggested that scores from ${\bf C}+{\bf I}$ were associated. Narrowing down the scope to each structural component, we also found that the individual components of UNIFAC and BIDIFAC, even though not indicative of subtype distinction in SWISS scores, could provide additional information on patients’ survival. BIDIFAC, accounting for normal tissues, provided reduced rank for individual components compared to UNIFAC, which would contribute to the increased power of the test.

5 Conclusion

In this paper we propose BIDIFAC, the first unified framework to handle multiple matrices that are shared both vertically and horizontally. In contrast to existing methods on joint matrix factorization, we provide the estimator based on soft-thresholding and nuclear norm penalization, where the complicated problem of selecting tuning parameters is alleviated by using the well-known result from random matrix theory. We conducted extensive simulation studies to show the efficiency and flexibility of BIDIFAC compared to existing methods. We applied our method to TCGA breast cancer data, where mRNA and miRNA profiles are obtained separately from both tumor tissues and normal tissues adjacent to tumors (NATs). From this application we conclude that (i) patterns of variability in normal tissue are largely also present in tumor tissue for both mRNA and miRNA, (ii) patterns that are associated with survival and clinical subtypes across both mRNA and miRNA are largely not present in normal tissue, and (iii) patterns that distinguish the clinical subtypes are shared by mRNA and miRNA. Existing vertical integration methods establish (iii), but not (i) or (ii). In addition to the integration of tumor and NAT data, this methodology may be applied to a growing number of applications with bidimensionally linked matrices. A particularly intriguing potential application is the integration of “pan-omics pan-cancer” data; that is, the integration of multi-omic data for samples from multiple types of cancers (as defined by their tissue-of-origin).

We describe several limitations of our method. These limitation primarily relate to use of the Frobenious norm in the objective and applying random matrix theory to select tuning parameters. Importantly, our model assumes normality of each matrix to select appropriate tuning parameters, which may be violated in practice. For example, application to SNPs data would be limited as each cell of a matrix takes discrete values (e.g., $0,1,2$). Similarly, we do not consider the case where one or more matrices are binary, which is common for biomedical data. Also, our approach is sensitive to outliers, due to the use of the Frobenious norm in the objective function. Even when the Gaussian assumption holds, our model may be extended by adopting variable sparsity or overcoming deficiencies of soft-thresholding, as neither soft-thresholding nor hard-thresholding provides optimal signal recovery (Shabalin and Nobel, 2013). Although our convex objective guarantees convergence to a stationary point that is empirically consistent for different starting values, the theoretical identifiability properties of the resulting decomposition deserves further study. Lastly, our approach does not inherently model uncertainty in the underlying structural decomposition and missing value imputations; the result gives the mode of a Bayesian posterior, and extensions to fully Bayesian approaches are worth considering.

Here we have focused on capturing four types of low-rank signals in bidimensional data: global, column-shared, row-shared, and individual. Other joint signals are possible. For example, when $p=3,q=1$, it is possible that ${\bf X}_{11}$ and ${\bf X}_{21}$ share signal that is not present in ${\bf X}_{31}$. Similarly, our method would suffer when ${\bf X}_{11},{\bf X}_{12},{\bf X}_{21}$ share signal that is not present in ${\bf X}_{22}$, though this is perhaps less likely for most data applications. It is straightforward to extend our framework to accommodate these or other structures, with appropriate additions to the objective (7).

An interesting extension of our work is the factorization of higher-order arrays (i.e., tensors). For example, integrating multiple omics profiles for multiple cohorts that are measured at multiple time points would give additional insights on the linked structure in a longitudinal manner that cross-sectional designs cannot provide.

6 Availability

We provide R functions for fitting BIDIFAC, imputing missing values, and generating simulation data as in Section 3, which are available at (https://github.com/lockEF/bidifac).

Acknowledgement

We would like to thank co-editor, associate editor, and two anonymous reviewers for their constructive comments. This research was supported by the National Institute of Health (NIH) under the grant R21CA231214 and by the Minnesota Supercomputing Institute (MSI).

Appendix A Proposed Algorithm (ISSVT)

The ISSVT algorithm, proceeds as follows proceeds as follows:

Initialize \widehat{\mbox{\boldmath\Theta}}=\{\widehat{{\bf G}}_{00},\widehat{{\bf R}}_{i0},\widehat{{\bf C}}_{0j},\widehat{{\bf I}}_{ij}~{}|~{}i=1,\dots,p,j=1,\dots,q\}. 2. 2.

For $i=1,\dots,p$ and $j=1,\dots,q$, apply Proposition 1 to obtain a closed form solution of the following:

1. (a)

$\widehat{{\bf G}}_{00}^{(new)}=\underset{{\bf G}_{00}:r({\bf G}_{00})\leq r_{00}}{\arg\min}\dfrac{1}{2}||{\bf X}_{00}-\widehat{{\bf R}}_{00}-\widehat{{\bf C}}_{00}-\widehat{{\bf I}}_{00}-{\bf G}_{00}||_{F}^{2}+\lambda_{00}||{\bf G}_{00}||_{*}$ 2. (b)

$\widehat{{\bf R}}_{i0}^{(new)}=\underset{{\bf R}_{i0}:r({\bf R}_{i0})\leq r_{i0}}{\arg\min}\dfrac{1}{2}||{\bf X}_{i0}-\widehat{{\bf G}}_{i0}^{(new)}-\widehat{{\bf C}}_{i0}-\widehat{{\bf I}}_{i0}-{\bf R}_{i0}||_{F}^{2}+\lambda_{i0}||{\bf R}_{i0}||_{*}$ 3. (c)

$\widehat{{\bf C}}_{0j}^{(new)}=\underset{{\bf C}_{0j}:r({\bf C}_{0j})\leq r_{0j}}{\arg\min}\dfrac{1}{2}||{\bf X}_{0j}-\widehat{{\bf G}}_{0j}^{(new)}-\widehat{{\bf R}}_{0j}^{(new)}-\widehat{{\bf I}}_{0j}-{\bf C}_{0j}||_{F}^{2}+\lambda_{0j}||{\bf C}_{0j}||_{*}$ 4. (d)

$\widehat{{\bf I}}_{ij}^{(new)}=\underset{{\bf I}_{ij}:r({\bf I}_{ij})\leq r_{ij}}{\arg\min}\dfrac{1}{2}||{\bf X}_{ij}-\widehat{{\bf G}}_{ij}^{(new)}-\widehat{{\bf R}}_{ij}^{(new)}-\widehat{{\bf C}}_{ij}^{(new)}-{\bf I}_{ij}||_{F}^{2}+\lambda_{ij}||{\bf I}_{ij}||_{*}$ 3. 3.

The algorithm converges if f_{2}(\widehat{\mbox{\boldmath\Theta}})-f_{2}(\widehat{\mbox{\boldmath\Theta}}^{(new)})<\epsilon. If it does not converge, replace \widehat{\mbox{\boldmath\Theta}}^{(new)} with \widehat{\mbox{\boldmath\Theta}} and repeat Step 2.

The algorithm iteratively minimizes the objective $f_{2}$ (6) over blocks ${\bf G}_{00}$, $\{{\bf R}_{i0}\mid i=1,\ldots,p\}$, $\{{\bf C}_{0j}\mid j=1,\ldots,q\}$, and $\{{\bf I}_{ij}\mid i=1,\dots,p,j=1,\dots,q\}$. By Proposition 1, this is equivalent to a blockwise coordinate descent algorithm for $f_{1}$ (4), with corresponding update blocks $\{{\bf U}_{i0}^{(G)},{\bf V}_{0j}^{(G)}\mid i=1,\dots,p,j=1,\dots,q\}$, $\{{\bf U}_{ij}^{(C)},{\bf V}_{0j}^{(C)}\mid i=1,\dots,p,j=1,\dots,q\}$, $\{{\bf U}_{i0}^{(R)},{\bf V}_{ij}^{(R)}\mid i=1,\dots,p,j=1,\dots,q\}$, $\{{\bf U}_{ij}^{(R)},{\bf V}_{ij}^{(I)}\mid i=1,\dots,p,j=1,\dots,q\}$.

Appendix B Proofs

Theorem 1.

Let

[TABLE]

minimize (4). Then,

[TABLE]

minimizes (6), where $\widehat{{\bf G}}_{ij}=\widehat{{\bf U}}_{i0}^{(G)}\widehat{{\bf V}}_{0j}^{(G)T},\widehat{{\bf R}}_{ij}=\widehat{{\bf U}}_{i0}^{(R)}\widehat{{\bf V}}_{ij}^{(R)T},\widehat{{\bf C}}_{ij}=\widehat{{\bf U}}_{ij}^{(C)}\widehat{{\bf V}}_{0j}^{(C)T},$ and $\widehat{{\bf I}}_{ij}=\widehat{{\bf U}}_{ij}^{(I)}\widehat{{\bf V}}_{ij}^{(I)T}$.

Proof.

Because \widehat{\mbox{\boldmath\Theta}}_{1} minimizes $f_{1}$,

[TABLE]

It follows from Lemma 1 that

[TABLE]

Analogous arguments show that

[TABLE]

Thus, f_{1}(\widehat{\mbox{\boldmath\Theta}}_{1})=2f_{2}(\widehat{\mbox{\boldmath\Theta}}_{2}).

Consider an alternative estimate \widetilde{\mbox{\boldmath\Theta}}_{2}=\{\widetilde{{\bf G}}_{ij},\widetilde{{\bf R}}_{ij},\widetilde{{\bf C}}_{ij},\widetilde{{\bf I}}_{ij}\mid i=1\ldots p,j=1,\ldots q\}. By Lemma 1, 2f_{2}(\widetilde{\mbox{\boldmath\Theta}}_{2})=f_{1}(\widetilde{\mbox{\boldmath\Theta}}_{1}) for some

[TABLE]

where $\widetilde{{\bf G}}_{ij}=\widetilde{{\bf U}}_{i0}^{(G)}\widetilde{{\bf V}}_{0j}^{(G)T},\widetilde{{\bf R}}_{ij}=\widetilde{{\bf U}}_{i0}^{(R)}\widetilde{{\bf V}}_{ij}^{(R)T},\widetilde{{\bf C}}_{ij}=\widetilde{{\bf U}}_{ij}^{(C)}\widetilde{{\bf V}}_{0j}^{(C)T},$ and $\widetilde{{\bf I}}_{ij}=\widetilde{{\bf U}}_{ij}^{(I)}\widetilde{{\bf V}}_{ij}^{(I)T}$. Thus, because \widehat{\mbox{\boldmath\Theta}}_{1} minimizes $f_{1}$,

[TABLE]

and we conclude that \widehat{\mbox{\boldmath\Theta}}_{2} minimizes $f_{2}$.

∎

Theorem 2.

The objective $f_{2}(\cdot)$ in (6) is convex over its domain.

Proof.

Consider \widetilde{\mbox{\boldmath\Theta}}^{(m)}=\{{\bf G}_{ij}^{(m)},{\bf R}_{ij}^{(m)},{\bf C}_{ij}^{(m)},{\bf I}_{ij}^{(m)}\mid i=1,\ldots,p,j=1,\ldots q\} for $m=1,2$ and $\alpha\in[0,1]$. It suffices to show

[TABLE]

where \alpha\widetilde{\mbox{\boldmath\Theta}}^{(1)}+(1-\alpha)\widetilde{\mbox{\boldmath\Theta}}^{(2)} is given by

[TABLE]

Decompose f_{2}(\mbox{\boldmath\Theta})=f_{2}^{\text{LS}}(\mbox{\boldmath\Theta})+f_{2}^{\text{PEN}}(\mbox{\boldmath\Theta}), where

[TABLE]

By convexity of the least squares objective

[TABLE]

for any ${\bf A}_{ij}^{(1)}$ and ${\bf A}_{ij}^{(2)}$. Applying (15) for each $(i,j)$, where ${\bf A}^{(m)}_{ij}={\bf G}^{(m)}_{ij}+{\bf R}^{(m)}_{ij}+{\bf C}^{(m)}_{ij}+{\bf I}^{(m)}_{ij}$ for $m=1,2$, gives

[TABLE]

By convexity of the nuclear norm operator,

[TABLE]

for any ${\bf A}^{(1)}$ and ${\bf A}^{(2)}$. Applying (17) to each additive term in $f_{2}^{\text{PEN}}$ gives

[TABLE]

Thus, (16) and (18) imply (14). ∎

Proposition 3.

The following conditions are necessary to allow for non-zero $\widehat{{\bf G}}_{00}$, $\widehat{{\bf R}}_{i0}$, $\widehat{{\bf C}}_{0j}$, and $\widehat{{\bf I}}_{ij}$:

$\mbox{max}_{j}\lambda_{ij}<\lambda_{i0}<\sum_{j}\lambda_{ij}$* for $i=1,\ldots,p$ and $\mbox{max}_{i}\lambda_{ij}<\lambda_{0j}<\sum_{i}\lambda_{ij}$ for $j=1,\ldots,q$* 2. 2.

$\mbox{max}_{i}\lambda_{i0}<\lambda_{00}<\sum_{i}\lambda_{i0}$* and $\mbox{max}_{j}\lambda_{0j}<\lambda_{00}<\sum_{j}\lambda_{0j}$.*

Proof.

Consider a violation of the left-hand inequality in condition 1.: $\lambda_{ij}\geq\lambda_{i0}$. Define \widehat{\mbox{\boldmath\Theta}} to be identical to \widetilde{\mbox{\boldmath\Theta}} but with $\widehat{{\bf I}}_{ij}=\mathbf{0}$ and $\widehat{{\bf R}}_{ij}=\widetilde{{\bf R}}_{ij}+\widetilde{{\bf I}}_{ij}$. By convexity of the nuclear norm, $||\widehat{{\bf R}}_{i0}||_{*}\leq||\widetilde{{\bf R}}_{i0}||_{*}+||\widetilde{{\bf I}}_{ij}||_{*}$, and it follows that

[TABLE]

Thus, regardless of the data ${\bf X}_{00}$, the objective f_{2}(\mbox{\boldmath\Theta}) is minimized with $\widehat{{\bf I}}_{ij}={\bf 0}_{m_{i}\times n_{j}}$. An analogous argument show that a violation of $\lambda_{ij}<\lambda_{0j}$ implies $\widehat{{\bf I}}_{ij}={\bf 0}_{m_{i}\times n_{j}}$ for some $i,j$. Moreover, analogous arguments show that a violation of $\max_{i}\lambda_{i0}<\lambda_{00}$ implies $\widehat{{\bf R}}_{i0}={\bf 0}_{m_{i}\times n_{0}}$ for some $i$, and that a violation of $\max_{j}\lambda_{0j}<\lambda_{00}$ implies $\widehat{{\bf C}}_{0j}={\bf 0}_{m_{0}\times n_{j}}$ for some $j$.

Now, consider a violation of the right-hand inequality of condition 1.: $\lambda_{i0}\geq\sum_{j}\lambda_{ij}$. Define \widehat{\mbox{\boldmath\Theta}} to be identical to \widetilde{\mbox{\boldmath\Theta}} but with $\widehat{{\bf R}}_{i0}={\bf 0}_{m_{i}\times n_{0}}$ and $\widehat{{\bf I}}_{ij}=\widetilde{{\bf R}}_{ij}+\widetilde{{\bf I}}_{ij}$ for $j=1,\ldots,q$. Then,

[TABLE]

Thus, regardless of the data ${\bf X}_{00}$, the objective f_{2}(\mbox{\boldmath\Theta}) is minimized with $\widehat{{\bf R}}_{i0}={\bf 0}_{m_{i}\times n_{0}}$. An analogous argument show that a violation of $\sum_{0j}<\sum_{i}\lambda_{ij}$ implies $\widehat{{\bf C}}_{0j}={\bf 0}_{m_{0}\times n_{j}}$ for some $j$. Moreover, analogous arguments show that a violation of $\lambda_{00}<\sum_{i}\lambda_{i0}$ or $\lambda_{00}<\sum_{j}\lambda_{0j}$ imply $\widehat{{\bf G}}_{00}={\bf 0}_{m_{0}\times n_{0}}$. $\square$ ∎

Appendix C Noiseless simulation

Here we generate joint and individual signals under the UNIFAC (vertically linked) model, and assess their recovery without residual error. We contrast the decomposition provided by the penalized objective $f_{2}(\cdot)$, with that obtained by alternative methods that enforce orthogonality of the estimated components.

We generate $\mathbf{X}_{1}:d_{1}\times n$ and $\mathbf{X}_{2}:d_{2}\times n$ as a sum of low-rank column-shared and individual structures that are independent but not necessarily orthogonal. That is,

[TABLE]

where

[TABLE]

and the entries of $\mathbf{U}_{1}^{(C)}:d_{1}\times r,\mathbf{U}_{2}^{(C)}:d_{1}\times r,\mathbf{V}:n\times r,\mathbf{U}_{1}^{(I)}:d_{1}\times r_{1},\mathbf{V}_{1}:n\times r_{2},\mathbf{U}_{2}^{(I)}:d_{2}\times r_{2},\mathbf{V}_{2}:n\times r_{2}$ are generated independently from a $\mathcal{N}(0,1)$ distribution.

We fix $r=r_{1}=r_{2}=10$ and generate $10$ datasets under each of four scenarios with different row and column dimensions: (1) $d_{1}=d_{2}=n=100$, (2) $d_{1}=d_{2}=100,n=500$, (3) $d_{1}=d_{2}=500,n=100$, (4) $d_{1}=d_{2}=n=500$. For each generated dataset we apply (i) UNIFAC, (ii) JIVE, (iii) AJIVE, and (iv) SLIDE, where the correctly specified ranks are used for methods (ii–iv). For UNIFAC, the noise variance is set to a small value ($\sigma=0.0001$) and the default penalties are used. In each case we compute the mean relative error in recovering underlying joint and individual structures:

[TABLE]

The results are summarized in Table 5. All methods decompose the underlying joint and individual signals with negligible error as the dimension of the sample size ($n$) and dimensions ($d=d_{2}=d_{2}$) increase. UNIFAC recovers the joint and individual signals with comparable or substantially improved accuracy across the four scenarios, despite the use of the correct ranks for the other three methods. The error in recovery for methods (ii-iv) is due to the orthogonality constraints, which are necessary for identifiability of the decomposition without additional penalization. The independent joint and individual structures are not exactly orthogonal, but will approach orthogonality as $n\rightarrow\infty$, and thus performance improves for $n=500$ vs $n=100$. JIVE and AJIVE both assume orthogonality of the rows of the joint and individual structure, ${\bf C}$ and ${\bf I}$; thus, their performance is comparable across scenarios and differences are solely due to imprecision of the computational algorithm for JIVE. SLIDE additionally assumes orthogonality of the individual structures ${\bf I}_{1}$ and ${\bf I}_{2}$, which results in slightly less accurate recovery for lower $n$. Even under independence the underlying joint and individual signals are not precisely orthogonal. The performance of UNIFAC for higher dimension $d$ demonstrates the potential to recover the true joint and individual signals more accurately by relaxing orthogonality constraints.

We also compare the recovery of joint and individual signals as their ranks increase relative to the data dimensions. We consider the scenario with $d_{1}=d_{2}=500,n=100$, and generate $10$ datasets for each of $r_{1}=r_{2}=r=\{1,2,\ldots,50\}$. For each dataset we estimate the decomposition via UNIFAC, AJIVE, or SLIDE; AJIVE is used instead of JIVE because in the noiseless scenario with given ranks they give the same underlying decomposition and AJIVE is more computationally efficient. The resulting mean relative errors (19) are shown in Figure 3. For ranks greater than $r=r_{1}=r_{2}=33$, the sum of the ranks of column-shared and individual structures ($r+r_{1}+r_{2}$) is greater than the rank of the observed signal: rank$([\mathbf{X}_{1}^{T},\mathbf{X}_{2}^{T}])\leq 100$. Thus, for these cases a SLIDE decomposition with the given ranks does not exist because the condition of orthogonality between $\hat{J}$, $\hat{I}_{1}$ and $\hat{I}_{2}$ cannot be satisfied. The recovery errors for the AJIVE decomposition also increase sharply at this point, while the trend remains stable for UNIFAC. In general, UNIFAC provides better recovery of the generated joint and individual signal as the ranks increase. However, the recovery error does increase steadily as the ranks get larger. Moreover, if the sum of the ranks is greater than the rank of the overall signal, this implies linear dependence among the underlying components (here, linear dependence among $\mathbf{V}$, $\mathbf{V}_{1}$, and $\mathbf{V}_{2}$), which complicates their interpretation.

Appendix D Data Analysis: Residual Diagnostics

Here we consider the distribution of residuals after the application of BIDIFAC for the TCGA application described in Section 4 of the main manuscript. That is, we consider the residual matrices

[TABLE]

BIDIFAC is best motivated when the error terms ${\bf E}_{ij}$ are approximately Gaussian, as the objective identifies the mode of a Bayesian model with a Gaussian likelihood (Section 2.7) and the selection of the tuning parameters is based on the assumption of Gaussian error (Section 2.6). Figure 4 shows the distribution of residuals for each of the four datasets considered (tumor mRNA, tumor miRNA, normal mRNA, normal miRNA) overlayed with two Gaussian densities: one giving the theoretical distribution of residuals implied by the pre-hoc estimate of the noise variance for each dataset ($\hat{\sigma}_{ij}^{MAD}$), and another giving the empirical Gaussian fit resulting from the sample mean and standard deviation of the observed residuals. None of the residual histograms show strong departures from Gaussianity, and the variance estimate used to tune the model ($\hat{\sigma}_{ij}^{MAD}$) fits the observed residuals reasonably well in each case.