Data-fusion using factor analysis and low-rank matrix completion

A common resolution to the identifiability issue in the file-matching problem is to assume that $\bm{Y}$ and $\bm{Z}$ are conditionally independent given $\bm{X}$. Under this assumption, ${\bm{\Sigma}}_{YZ}={\bm{\Sigma}}_{YX}{\bm{\Sigma}}_{XX}^{-1}{\bm{\Sigma}}_{XZ}$. The conditional independence assumption is commonly invoked to justify statistical matching (D’Orazio et al., 2006a; Rässler, 2002). However, it is untestable and can have a large influence on downstream results (Barry, 1988; Rodgers, 1984).

To avoid making assumptions, the generative model can be treated as being partially identified (Gustafson, 2015). There is a growing literature focused on developing uncertainty bounds for the matching problem that reflect the limited information in the sample (D’Orazio et al., 2006b; Conti et al., 2012, 2016). For the multivariate normal model, the goal is to estimate the a feasible set of values for $\bm{\bm{\Sigma}}_{YZ}$, rather than to deliver a point estimate (Kadane, 2001; Moriarity and Scheuren, 2001; Ahfock et al., 2016). The restriction that the covariance matrix be positive definite creates an identified set

$\displaystyle\Theta_{\text{I}}$

$\displaystyle=\left\{\bm{\bm{\Sigma}}_{YZ}\in\mathbb{R}^{p_{Y}\times p_{Z}}:\begin{pmatrix}\bm{\bm{\Sigma}}_{XX}&\bm{\bm{\Sigma}}_{XY}&\bm{\bm{\Sigma}}_{XZ}\\ \bm{\bm{\Sigma}}_{YX}&\bm{\bm{\Sigma}}_{YY}&\bm{\bm{\Sigma}}_{YZ}\\ \bm{\bm{\Sigma}}_{ZX}&\bm{\bm{\Sigma}}_{ZY}&\bm{\bm{\Sigma}}_{ZZ}\end{pmatrix}\text{is positive definite}\right\}.$

(3)

A disadvantage of the partial identification approach is that the level of uncertainty on $\bm{\Sigma}$ can be very high. Bounds on the parameters can be too wide for any practical interpretation.

Low-rank matrix completion may be applied to fill in the partial estimate of the covariance matrix or to impute the missing observations in the partially observed data matrix. To complete the estimate of the covariance matrix $\widehat{\bm{\Sigma}}$, one may consider the regularised least-squares objective

	$\displaystyle f(\bm{M},\bm{\Psi})$	$\displaystyle=\dfrac{1}{2}\sum_{(i,j)\in\Omega}(\widehat{\bm{\Sigma}}_{ij}-\bm{M}_{ij}-\bm{\Psi}_{ij})^{2}+\lambda\lVert\bm{M}\rVert_{*},$
		$\displaystyle\text{subject to }\bm{M}^{\mathsf{T}}=\bm{M},$

where $\lVert\bm{M}\rVert_{*}$ is the nuclear norm of the matrix $\bm{M}$ (Candes and Plan, 2010; Mazumder et al., 2010) and $\Omega$ contains the indices of the known elements of $\widehat{\bm{\Sigma}}$. The main barrier to attempting matrix completion on $\widehat{\bm{\Sigma}}$ is that the presence of the diagonal matrix ${\bm{\Psi}}$ along with the symmetry constraint on $\bm{M}$ necessitates modifications to existing algorithms for low-rank matrix completion. The algorithm for symmetric low-rank matrix completion with block missingness developed in Bishop and Byron (2014) could be a useful starting point.

A more natural approach in the file-matching problem is to use a low-rank matrix completion algorithm to impute the missing data. Partitioning the complete-dataset $\mathcal{D}$ as

$\displaystyle\mathcal{D}$

$\displaystyle=\begin{pmatrix}\bm{X}_{A}&\bm{Y}_{A}&\bm{Z}_{A}\\ \bm{X}_{B}&\bm{Y}_{B}&\bm{Z}_{B}\end{pmatrix}$

where the first row corresponds to dataset A, and the second row corresponds to dataset B, the missing $\bm{Z}_{A}$ and $\bm{Y}_{B}$ observations can be estimated by applying matrix completion to $\mathcal{D}$. The modelling assumption is that complete-dataset can be expressed as $\mathcal{D}=\mathcal{D}_{0}+\mathcal{E}$ for some rank $q$ matrix $\mathcal{D}_{0}$ plus noise $\mathcal{E}$. To estimate $\mathcal{D}_{0}$, the following least squares objective can be used

$\displaystyle f(\bm{G},\bm{H})$

$\displaystyle=\dfrac{1}{2}\sum_{(i,j)\in\Omega}(\mathcal{D}_{ij}-[\bm{G}\bm{H}]_{ij})^{2},$

where $\bm{G}$ is a $n\times q$ matrix, $\bm{H}$ is a $q\times p$ matrix and $\Omega$ contains the indices of the observed elements in $\mathcal{D}$. With $(\widehat{\bm{G}},\widehat{\bm{H}})\in\operatorname*{argmin}_{(\bm{G},\bm{H})}f(\bm{G},\bm{H})$, the low-rank estimate of $\mathcal{D}_{0}$ is

$\displaystyle\widehat{\mathcal{D}}_{0}$

$\displaystyle=\widehat{\bm{G}}\widehat{\bm{H}}=\begin{pmatrix}\widehat{\bm{X}}_{A}&\widehat{\bm{Y}}_{A}&\widehat{\bm{Z}}_{A}\\ \widehat{\bm{X}}_{B}&\widehat{\bm{Y}}_{B}&\widehat{\bm{Z}}_{B}\end{pmatrix}.$

(4)

Using the low-rank reconstructions $\widehat{\bm{Z}}_{A}$ and $\widehat{\bm{Y}}_{B}$ as imputed values for the missing-data $\bm{Z}_{A}$ and $\bm{Y}_{B}$, an estimate of $\bm{\Sigma}_{YZ}$ is then given by

$\displaystyle\widehat{\bm{\Sigma}}_{YZ}$

$\displaystyle=\text{cov}\left(\begin{bmatrix}{\bm{Y}}_{A}\\ \widehat{\bm{Y}}_{B}\end{bmatrix},\begin{bmatrix}\widehat{\bm{Z}}_{A}\\ {\bm{Z}}_{B}\end{bmatrix}\right).$

The theoretical work in Koltchinskii et al. (2011) can be used to give finite-sample bounds on the reconstruction error $\lVert\mathcal{D}_{0}-\widehat{\mathcal{D}}_{0}\rVert_{2}^{2}$ under assumptions on the noise matrix $\mathcal{E}$. However, in the context of the file-matching problem, it is not immediately clear from existing theoretical work whether exact recovery of the true covariance $\bm{\Sigma}_{YZ}$ is possible in the large-sample limit.

A probabilistic factor analysis model with $q$ factors has the latent variable representation

$\displaystyle\begin{pmatrix}\bm{X}\\ \bm{Y}\\ \bm{Z}\end{pmatrix}$

$\displaystyle=\begin{pmatrix}\bm{\Lambda}_{X}\\ \bm{\Lambda}_{Y}\\ \bm{\Lambda}_{Z}\end{pmatrix}\bm{F}+\bm{\epsilon},$

(5)

where $\bm{F}\sim\mathcal{N}(\bm{0},\bm{I}_{q})$, $\bm{\epsilon}\sim\mathcal{N}(\bm{0},\bm{\Psi})$, and

$\displaystyle\bm{\Psi}=\begin{pmatrix}\bm{\Psi}_{X}&\bm{0}&\bm{0}\\ \bm{0}&\bm{\Psi}_{Y}&\bm{0}\\ \bm{0}&\bm{0}&\bm{\Psi}_{Z}\\ \end{pmatrix}.$

Kamakura and Wedel (2000) discuss the application of factor analysis models to the statistical file-matching problem, but do not consider the identifiability of the model, or develop an algorithm for parameter estimation. We address both of these important points. We prove the identifiability of the factor analysis model in the file-matching problem under the following assumptions.

Assumption 2 is based on the sufficient conditions for the identifiability of the factor analysis model given in Anderson and Rubin (1956). For the file-matching-problem, Assumption 2 enforces the requirement that $q<(p_{X}+p_{Y})/2$ and $q<(p_{X}+p_{Z})/2$. Assumption 1 requires that $q\leq p_{X}$.

The key to recovery of $\bm{\Sigma}_{YZ}$ is that the matrix $\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}$ is of rank $q\leq p_{X}$. Lemma 1 shows that it is possible to complete the matrix $\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}$ given the missing-data pattern in Table 1.

The proof is given in the Appendix. Bishop and Byron (2014) establish more general results for the low-rank completion of symmetric matrices with structured missingness. Although Lemma 1 is sufficient for low-rank matrix completion problems, identifiability of the factor analysis model is a more complex issue. Specifically, for two parameter sets $\{\bm{\Lambda},\bm{\Psi}\}$ and $\{\bm{\Lambda}^{*},\bm{\Psi}^{*}\}$, the equality $\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi}=\bm{\Lambda}^{*}\bm{\Lambda}^{*\mathsf{T}}+\bm{\Psi}^{*}$ does not imply that $\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}=\bm{\Lambda}^{*}\bm{\Lambda}^{*\mathsf{T}}$ (Anderson and Rubin, 1956; Shapiro, 1985). Assumption 2 allows us to make a statement about the identifiability of the factor analysis model in the file-matching problem.

The proof is given in the Appendix. In terms of the generative model in Section 3.1, Theorem 1 states that if parameter sets $\bm{\Theta}_{1}$ and $\bm{\Theta}_{2}$ satisfy $f(\bm{x},\bm{y};\bm{\Theta}_{1})=f(\bm{x},\bm{y};\bm{\Theta}_{2})$ and $f(\bm{x},\bm{z};\bm{\Theta}_{1})=f(\bm{x},\bm{z};\bm{\Theta}_{2})$ then the factor loadings in $\bm{\Theta}_{1}$ and $\bm{\Theta}_{2}$ are equal up to post-multiplication on the right by an orthogonal matrix. The conclusion is that recovery of $\bm{\Sigma}_{YZ}$ in the file-matching scenario is possible given the observable sub-blocks of the population covariance matrix $\bm{\Sigma}$. Assuming that the sample covariance matrices converge in probability to the population quantities as $n_{A},n_{B}\to\infty$, consistent estimation of $\bm{\Sigma}_{YZ}$ is possible.

Under mild regularity conditions on the factor loadings, $q<p/2$ is sufficient for identifiability of the factor analysis model with complete-cases (Anderson and Rubin, 1956; Shapiro, 1985). In the file-matching scenario, Assumption 2 enforces the stronger condition that $q<(p_{X}+p_{Y})/2$ and $q<(p_{X}+p_{Z})/2$. If the number of factors $q$ falls within the intervals $(p_{X}+p_{Y})/2<q<p/2$ and $(p_{X}+p_{Z})/2<q<p/2$, Assumption 2 is violated. However, one may hope that the model is still identifiable, given that consistent estimation is possible with complete-cases.

Some insight can be gained by comparing the number of equations and the number of unknowns (Ledermann, 1937). There are $pq$ unknown parameters in $\bm{\Lambda}$ and $p$ unknown parameters in $\bm{\Psi}$. Given the symmetry in $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi}$, there are $p(p+1)/2$ equations, and the adoption of a rotation constraint on $\bm{\Lambda}$ yields an additional $q(q-1)/2$ equations (Anderson and Rubin, 1956). The number of equations minus the number of unknowns is $C=[(p-q)^{2}-p-q]/2$. A useful rule of thumb is to expect the the factor analysis model to be identifiable if $C\geq 0$, and nonidentifiable if $C<0$ (Ledermann, 1937; Anderson and Rubin, 1956). Bekker and ten Berge (1997) show that $C>0$ is sufficient for the identifiability of the factor model in most circumstances. In the file-matching problem, the equality constraints due to $\bm{\Sigma}_{YZ}$ are not present and the degrees of freedom is now reduced to $C_{M}=[(p-q)^{2}-p-q]/2-p_{Y}p_{Z}$. If $C_{M}$ is negative and $C$ is positive it may signal a situation where the missing-data has caused the factor model to be nonidentifiable.

Table 2 reports the maximum number of allowable factors $q$ using the degrees of freedom criteria $C$ and $C_{M}$, and the restrictions of Assumption 2 in a file-matching problem with $p_{X}=p_{Y}=p_{Z}=p/3$. Assumption 2 places stronger requirements on the number of factors $q$ than the degrees of freedom requirement $C_{M}$. The number of allowable factors using $C_{M}$ is smaller than $C$, demonstrating the information loss due to the missing-data.

As mentioned in Section 2.3, procedures for low-rank matrix completion focus on the estimation of $\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}$ rather than the joint estimation of $\bm{\Lambda}$ and $\bm{\Psi}$. A maximum likelihood approach to the factor analysis model facilitates joint estimation of $\bm{\Lambda}$ and $\bm{\Psi}$. With complete-cases, the EM updates for factor analysis can be written in closed form (Rubin and Thayer, 1982). Let $\bm{S}$ denote the sample scatter matrix

$\displaystyle\bm{S}$

$\displaystyle=\sum_{i=1}^{n}\begin{pmatrix}\bm{x}_{i}\\ \bm{y}_{i}\\ \bm{z}_{i}\end{pmatrix}\begin{pmatrix}\bm{x}_{i}\\ \bm{y}_{i}\\ \bm{z}_{i}\end{pmatrix}^{\mathsf{T}}.$

(8)

Given current parameter estimates of $\bm{\bm{\Lambda}}$ and $\bm{\bm{\Psi}}$, define $\bm{\beta}\equiv\bm{\Lambda}^{\mathsf{T}}(\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi})^{-1}$. The M-step is

	$\displaystyle\bm{\Lambda}^{\text{new}}$	$\displaystyle=(\bm{S}\bm{\beta}^{\mathsf{T}})\left(nI-n\bm{\beta}\bm{\Lambda}+\bm{\beta}\bm{S}\bm{\beta}^{\mathsf{T}}\right)^{-1},$		(9)
	$\displaystyle\bm{\Psi}^{\rm{new}}$	$\displaystyle=\dfrac{1}{n}\text{diag}\left(\bm{S}-\bm{\Lambda}^{\text{new}}\bm{\beta}\bm{S}\right).$		(10)

For the file-matching problem we need to compute the conditional expectation of (8) given the missing-data pattern. The observed scatter matrices in dataset A and dataset B are

	$\displaystyle\bm{P}$	$\displaystyle=\sum_{i=1}^{n_{A}}\begin{pmatrix}\bm{x}_{i}\\ \bm{y}_{i}\end{pmatrix}\begin{pmatrix}\bm{x}_{i}\\ \bm{y}_{i}\end{pmatrix}^{\mathsf{T}}=\begin{pmatrix}\bm{P}_{XX}&\bm{P}_{XY}\\ \bm{P}_{YX}&\bm{P}_{YY}\end{pmatrix},$		(11)
	$\displaystyle\quad\bm{T}$	$\displaystyle=\sum_{i=1}^{n_{B}}\begin{pmatrix}\bm{x}_{i}\\ \bm{z}_{i}\end{pmatrix}\begin{pmatrix}\bm{x}_{i}\\ \bm{z}_{i}\end{pmatrix}^{\mathsf{T}}=\begin{pmatrix}\bm{T}_{XX}&\bm{T}_{XZ}\\ \bm{T}_{ZX}&\bm{T}_{ZZ}\end{pmatrix}.$		(12)

To determine the conditional expectation of the complete-data sufficient statistic, for current parameter estimates of $\bm{\Lambda}$ and $\bm{\Psi}$ define

	$\displaystyle\bm{\omega}$	$\displaystyle=\begin{pmatrix}\bm{\omega}_{X}&\bm{\omega}_{Y}\end{pmatrix}=\begin{pmatrix}\bm{\Lambda}_{Z}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Z}\bm{\Lambda}_{Y}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Y}^{\mathsf{T}}\\ \bm{\Lambda}_{Y}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Y}\bm{\Lambda}_{Y}^{\mathsf{T}}+\bm{\Psi}_{Y}\end{pmatrix}^{-1},$
	$\displaystyle\bm{\alpha}$	$\displaystyle=\begin{pmatrix}\bm{\alpha}_{X}&\bm{\alpha}_{Z}\end{pmatrix}=\begin{pmatrix}\bm{\Lambda}_{Y}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Y}\bm{\Lambda}_{Z}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Z}^{\mathsf{T}}\\ \bm{\Lambda}_{Z}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Z}\bm{\Lambda}_{Z}^{\mathsf{T}}+\bm{\Psi}_{Z}\end{pmatrix}^{-1},$

and the conditional covariance matrices

	$\displaystyle\bm{\Sigma}_{Z\mid XY}$	$\displaystyle=\begin{pmatrix}\bm{\Lambda}_{Z}\bm{\Lambda}_{Z}^{\mathsf{T}}+\bm{\Psi}_{Z}\end{pmatrix}-\begin{pmatrix}\bm{\Lambda}_{Z}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Z}\bm{\Lambda}_{Y}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Y}^{\mathsf{T}}\\ \bm{\Lambda}_{Y}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Y}\bm{\Lambda}_{Y}^{\mathsf{T}}+\bm{\Psi}_{Y}\end{pmatrix}^{-1}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{Z}^{\mathsf{T}}\\ \bm{\Lambda}_{Y}\bm{\Lambda}_{Z}^{\mathsf{T}}\end{pmatrix},$
	$\displaystyle\bm{\Sigma}_{Y\mid XZ}$	$\displaystyle=\begin{pmatrix}\bm{\Lambda}_{Y}\bm{\Lambda}_{Y}^{\mathsf{T}}+\bm{\Psi}_{Y}\end{pmatrix}-\begin{pmatrix}\bm{\Lambda}_{Y}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Y}\bm{\Lambda}_{Z}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Z}^{\mathsf{T}}\\ \bm{\Lambda}_{Z}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Z}\bm{\Lambda}_{Z}^{\mathsf{T}}+\bm{\Psi}_{Z}\end{pmatrix}^{-1}\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{Y}^{\mathsf{T}}\\ \bm{\Lambda}_{Z}\bm{\Lambda}_{Y}^{\mathsf{T}}\end{pmatrix}.$

Using the augmented scatter matrices

	$\displaystyle\widetilde{\bm{P}}$	$\displaystyle=\begin{pmatrix}\bm{P}_{XX}&\bm{P}_{XY}&\bm{P}_{XX}\bm{\omega}_{X}^{\mathsf{T}}+\bm{P}_{XY}\bm{\omega}_{Y}^{\mathsf{T}}\\ \bm{P}_{YX}&\bm{P}_{YY}&\bm{P}_{YX}\bm{\omega}_{X}^{\mathsf{T}}+\bm{P}_{YY}\bm{\omega}_{Y}^{\mathsf{T}}\\ \bm{\omega}_{Y}\bm{P}_{YX}+\bm{\omega}_{X}\bm{P}_{XX}&\bm{\omega}_{Y}\bm{P}_{YY}+\bm{\omega}_{X}\bm{P}_{XY}&\bm{\omega}\bm{P}\bm{\omega}^{\mathsf{T}}+n_{A}\bm{\Sigma}_{Z\mid XY}\end{pmatrix},$
	$\displaystyle\widetilde{\bm{T}}$	$\displaystyle=\begin{pmatrix}\bm{T}_{XX}&\bm{T}_{XX}\bm{\alpha}_{X}^{\mathsf{T}}+\bm{T}_{XZ}\bm{\alpha}_{Z}^{\mathsf{T}}&\bm{T}_{XZ}\\ \bm{\alpha}_{Z}\bm{T}_{ZX}+\bm{\alpha}_{X}\bm{T}_{XX}&\bm{\alpha}\bm{T}\bm{\alpha}^{\mathsf{T}}+n_{B}\bm{\Sigma}_{Y\mid XZ}&\bm{\alpha}_{Z}\bm{T}_{ZZ}+\bm{\alpha}_{X}\bm{T}_{XZ}\\ \bm{T}_{ZX}&\bm{T}_{ZX}\bm{\alpha}_{X}^{\mathsf{T}}+\bm{T}_{ZZ}\bm{\alpha}_{Z}^{\mathsf{T}}&\bm{T}_{ZZ}\end{pmatrix},$

the expected value of the complete-data sufficient statistic is given by

$\displaystyle\widetilde{\bm{S}}$

$\displaystyle=\widetilde{\bm{P}}+\widetilde{\bm{T}}.$

Given current parameter estimates of ${\bm{\Lambda}},{\bm{\Psi}}$, define $\bm{\beta}\equiv\bm{\Lambda}^{\mathsf{T}}(\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi})^{-1}$. The M-step is then

	$\displaystyle\bm{\Lambda}^{\text{new}}$	$\displaystyle=(\widetilde{\bm{S}}\bm{\beta}^{\mathsf{T}})\left(n\bm{I}-n\bm{\beta}\bm{\Lambda}+\bm{\beta}\widetilde{\bm{S}}\bm{\beta}^{\mathsf{T}}\right)^{-1},$		(13)
	$\displaystyle\bm{\Psi}^{\rm{new}}$	$\displaystyle=\dfrac{1}{n}\text{diag}\left(\widetilde{\bm{S}}-\bm{\Lambda}^{\text{new}}\bm{\beta}\tilde{\bm{S}}\right).$		(14)

The EM algorithm can be run using the sample covariance matrices from dataset A and dataset B, and so it is not necessary to have the original observations. The multivariate normality assumption used to develop the algorithm is not crucial, as long as the population covariance matrix fits the factor analysis model $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi}$ (Browne, 1984).

In practice, the number of factors $q$ must be determined. The Bayesian information criterion (Schwarz, 1978) is a generic tool for model selection that has been well explored in factor analysis for selecting $q$ (Preacher et al., 2013). The Bayesian information criterion for a model $\mathcal{M}$ with likelihood function $L$ and parameter $\bm{\theta}$ is given by

$\displaystyle\text{BIC}(\mathcal{M})$

$\displaystyle=-2\log L(\bm{\widehat{\theta}})+d\log n,$

(15)

where $\widehat{\bm{\theta}}$ is the maximum likelihood estimate, $d$ is the number of free parameters in the model $\mathcal{M}$, and $n$ is the sample size (Schwarz, 1978). For factor analysis with complete-cases the log-likelihood is

$\displaystyle\log L(\bm{\Lambda},\bm{\Psi})$

$\displaystyle=-\dfrac{n}{2}\log\lvert\bm{\Sigma}\rvert-\dfrac{1}{2}\text{trace}\{\bm{\Sigma}^{-1}\bm{S}\}-\dfrac{np}{2}\log 2\pi,$

where $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi}$ and $\bm{S}$ is the sample scatter matrix (8). The Bayesian information criterion for factor analysis with complete-cases is therefore

$\displaystyle\text{BIC}(q)$

$\displaystyle=-2\log L(\widehat{\bm{\Lambda}},\widehat{\bm{\Psi}})+(qp+p-q(q-1)/2)\log n,$

where the number of free parameters is obtained by subracting the number of rotation constraints $q(q-1)/2$ from the nominal number of parameters $qp+p$.

For model selection in missing-data problems it has been suggested to use the observed-data log-likelihood in place of the complete-data log-likelihood in the standard definition of the Bayesian information criterion (15) (Ibrahim et al., 2008). This is the strategy we propose to use for determining the number of factors $q$ in the file-matching problem. The log-likelihood functions for each dataset are given by

	$\displaystyle\log L_{A}(\bm{\Sigma},\bm{\Psi})$	$\displaystyle=-\dfrac{n_{A}}{2}\log\lvert\bm{\Sigma}_{A}\rvert-\dfrac{1}{2}\text{trace}\{\bm{\Sigma}_{A}^{-1}\bm{P}\}-\dfrac{n_{A}p_{A}}{2}\log 2\pi,$
	$\displaystyle\log L_{B}(\bm{\Sigma},\bm{\Psi})$	$\displaystyle=-\dfrac{n_{B}}{2}\log\lvert\bm{\Sigma}_{B}\rvert-\dfrac{1}{2}\text{trace}\{\bm{\Sigma}_{B}^{-1}\bm{T}\}-\dfrac{n_{B}p_{B}}{2}\log 2\pi,$

where $p_{A}=p_{X}+p_{Y}$, $p_{B}=p_{X}+p_{Z}$, $\bm{P}$ is the observed scatter matrix from dataset A (11), $\bm{T}$ is the observed scatter matrix from dataset B (12) and

	$\displaystyle\bm{\Sigma}_{A}$	$\displaystyle=\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Y}^{\mathsf{T}}\\ \bm{\Lambda}_{Y}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Y}\bm{\Lambda}_{Y}^{\mathsf{T}}+\bm{\Psi}_{Y}\end{pmatrix},$
	$\displaystyle\bm{\Sigma}_{B}$	$\displaystyle=\begin{pmatrix}\bm{\Lambda}_{X}\bm{\Lambda}_{X}^{\mathsf{T}}+\bm{\Psi}_{X}&\bm{\Lambda}_{X}\bm{\Lambda}_{Z}^{\mathsf{T}}\\ \bm{\Lambda}_{Z}\bm{\Lambda}_{X}^{\mathsf{T}}&\bm{\Lambda}_{Z}\bm{\Lambda}_{Z}^{\mathsf{T}}+\bm{\Psi}_{Z}\end{pmatrix}.$

The Bayesian information criterion using the observed-data log-likelihood is then

$\displaystyle\text{BIC}(q)$

$\displaystyle=-2\{\log L_{A}(\widehat{\bm{\Lambda}},\widehat{\bm{\Psi}})+\log L_{B}(\widehat{\bm{\Lambda}},\widehat{\bm{\Psi}})\}+(qp+p-q(q-1)/2)\log n.$

In this section we compare the performance of various algorithms for the estimation of $\bm{\Sigma}_{YZ}$. Existing algorithms for low-rank matrix completion were implemented using the R package filling (You, 2020). We also considered the standard conditional independence assumption (CIA). We tested the algorithms on the following datasets.

• •

  Sachs dataset. This dataset consists of results from 11 flow cytometry experiments measuring the same set of 12 variables. The measurements were log transformed. Data from each experiment was centered, then combined and scaled so that each variable had unit variance.

• •

  Wisconsin breast cancer dataset. The original observations in the malignant class were centered and scaled so that each variable had unit variance.

• •

  Abalone dataset. The original observations in the female class were centered and scaled so that each variable had unit variance.

For each dataset we considered matching scenarios with $p_{X}$ common variables, $p_{Y}$ unique variables in dataset A, and $p_{Z}$ unique variables in dataset B. This defines partitions of the complete-dataset. Table 3 gives the file-matching scenarios considered for each dataset, as well as the number of factors $q$ that were used in the fitting of the factor analysis model. The final columns indicate whether the degrees of freedom criteria $C$ and $C_{M}$, and Assumption 2 are satisfied for the chosen number of factors $q$. Assumption 2 is violated on the Sachs dataset, but the degrees of freedom $C_{M}$ is positive.

In each simulation, variables were randomly allocated to the $\bm{X}$, $\bm{Y}$, and $\bm{Z}$ groups. We then used different algorithms to estimate $\bm{\Sigma}_{YZ}$. The mean squared error over elements of $\bm{\Sigma}_{YZ}$ was then recorded for each algorithm, that is $\lVert\widehat{\bm{\Sigma}}_{YZ}-\bm{\Sigma}_{YZ}\rVert^{2}/(p_{Y}p_{Z})$. This process was repeated 100 times. The EM algorithm was run for 2000 iterations in each simulation. We considered two different initial value settings for the EM algorithm, initialisation using the factor analysis solution from the full covariance matrix, and random initialisation of $\bm{\Lambda}$ by sampling from a standard normal distribution. For the random initialisation protocol we took one hundred samples of $\bm{\Lambda}$ followed by a short run of the EM algorithm for fifty iterations. The parameters with the highest log-likelihood after the fifty iterations were then used in the longer run. In each simulation we also fit a factor model using the complete-dataset with no missing data to provide a reference point for the goodness-of-fit of the factor model.

Figure 2 compares the results using the different algorithms. The results for factor analysis model with the favourable initial values are shown as FM, the results for the factor analysis model with random initialisation are shown as FM (random). The factor model with good initial values (FM) has the lowest median error across each dataset (Table 4). The errors on the Sachs dataset and the Wisconsin breast cancer dataset are particularly small.

The results for the factor model with random initialisation (FM (random)) are very similar to those of FM using the favourable initial values except on the Sachs dataset. On the Sachs dataset, Assumption 2 is violated and there may be multiple factor solutions. In Figure 2 (a), the errors for the factor model with random initialisation (FM (random)) are much more dispersed than the factor model with the good initial values (FM). Although there are sufficient degrees of freedom $(C_{M}\geq 0)$, it appears that the model may not be identifiable. With good initial values, the EM algorithm appears to converge to a local mode that gives a good estimate of $\Sigma_{YZ}$.

The interquartile range of the error for each algorithm over the one hundred variable permutations is given in Table 4. The factor analysis model (FM) has the smallest interquartile range on each dataset, closely followed by the factor model using random initialisations with the exception of the Sachs dataset.

We also performed some experiments to assess the behavior of the Bayesian information criterion for determining the number of factors $q$. The following datasets and variable partitions were used in the experiments

• •

  Reise dataset. This dataset is a $16\times 16$ correlation matrix for mental health items. The variables were partitioned as $p_{X}=6,p_{Y}=5$, and $p_{Z}=5$.

• •

  Harman dataset. This dataset is a $24\times 24$ correlation matrix for psychological tests. The variables were partitioned as $p_{X}=5,p_{Y}=10$, and $p_{Z}=9$.

• •

  Holzinger dataset. This dataset is a $14\times 14$ correlation matrix for mental ability scores. The variables were partitioned as $p_{X}=6,p_{Y}=4$, and $p_{Z}=4$.

• •

  Simulated from a factor model with $q=6$ and $p=100$. Elements of $\Lambda$ were sampled from a $\mathcal{N}(2,1)$ distribution. The $j$th diagonal element of $\Psi$ was set as $D_{j}^{2}$, where $D_{j}\sim\mathcal{N}(3,0.01)$. The covariance matrix $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\mathsf{T}}+\bm{\Psi}$ was then standardised to a correlation matrix. We simulated one dataset with $n_{A}=n_{B}=50$ observations and one dataset of $n_{A}=n_{B}=500$ observations. The variables were partitioned as $p_{X}=6,p_{Y}=47,p_{Z}=47$.

In each simulation, variables were randomly allocated to the $\bm{X}$, $\bm{Y}$ and $\bm{Z}$ groups with $n_{A}=n_{B}=n/2$. We then computed the BIC for a range of values for $q$ and recorded the optimal value $q^{*}$. This process was repeated 100 times. In each of the simulations we considered the range of values for $q$ such that Assumptions 1 and 2 were satisfied. Table 5 reports the number of times a $q$ factor model was chosen using the BIC over the 100 replications. The optimal number of factors according to the BIC using complete-cases is given as $q_{c}^{*}$ in Table 5. The missing-data appears to lead to the BIC acting conservatively, selecting a number of factors less than or equal to the number chosen with complete-cases. The correct number of factors $q=6$ is selected for the simulated datasets when using complete-cases. For the simulated dataset with $n=50$, the BIC selects fewer than six factors in each trial. For the simulated dataset with $n=500$ the BIC selects the correct number of factors $q=6$ in each replication. Figure 3 shows boxplots of the estimation error $\lVert\widehat{\bm{\Sigma}}_{YZ}-\bm{\Sigma}_{YZ}\rVert^{2}/(p_{Y}p_{Z})$ for different values of $q$. The errors from using the conditional independence assumption and from fitting a factor model with complete-cases are also reported for comparison. For almost all values of $q$, the median error using the factor model is lower than that of the conditional independence model.