Statistical matching of non-Gaussian data

TL;DR

This paper introduces computational methods for statistically matching non-Gaussian data across datasets with structured missingness, addressing limitations of nearest-neighbour imputation through model-based approaches.

Contribution

It develops feasible procedures for non-Gaussian data matching using data augmentation and identifiability constraints, improving over traditional nearest-neighbour methods.

Findings

Model-based matching addresses weaknesses of nearest-neighbour imputation.

Approach is validated on flow cytometry datasets.

Provides a robust alternative for non-Gaussian data integration.

Abstract

The statistical matching problem is a data integration problem with structured missing data. The general form involves the analysis of multiple datasets that only have a strict subset of variables jointly observed across all datasets. The simplest version involves two datasets, labelled A and B, with three variables of interest $X, Y$ and $Z$. Variables $X$ and $Y$ are observed in dataset A and variables $X$ and $Z$ are observed in dataset $B$. Statistical inference is complicated by the absence of joint $(Y, Z)$ observations. Parametric modelling can be challenging due to identifiability issues and the difficulty of parameter estimation. We develop computationally feasible procedures for the statistical matching of non-Gaussian data using suitable data augmentation schemes and identifiability constraints. Nearest-neighbour imputation is a common alternative technique due to its ease of use and generality. Nearest-neighbour matching is based on a conditional independence assumption that may be inappropriate for non-Gaussian data. The violation of the conditional independence assumption can lead to improper imputations. We compare model based approaches to nearest-neighbour imputation on a number of flow cytometry datasets and find that the model based approach can address some of the weaknesses of the nonparametric nearest-neighbour technique.

Tables

Table 1. Table 1: Missing data pattern in the statistical matching problem.

Dimensions	$𝑿$	$𝒀$	$𝒁$
Dataset A	✓	✓	✗
Dataset B	✓	✗	✓

Table 2. Table 2: Matching scenarios using the iris dataset. Clustering using the 𝑿 𝑿 \bm{X} variables is substantially easier in Scenario 2.

	$𝑿$	$𝒀$	$𝒁$
Scenario 1	Sepal Width	Petal Length	Petal Width
Scenario 2	Sepal Width, Sepal Length	Petal Length	Petal Width

Table 3. Table 3: Matching scenarios using the two cell subpopulation flow cytometry dataset. Clustering using the 𝑿 𝑿 \bm{X} variables is substantially easier in Scenario 2.

	$𝑿$	$𝒀$	$𝒁$
Scenario 1	FL1.H	FL4.H	FL2.H
Scenario 2	FL3.H	FL4.H	FL2.H

Table 4. Table 4: Estimates of the ( 𝒀 , 𝒁 ) 𝒀 𝒁 (\bm{Y},\bm{Z}) correlation ρ Y ​ Z subscript 𝜌 𝑌 𝑍 \rho_{YZ} for the two cell subpopulation flow cytometry dataset.

	Scenario 1		Scenario 2
Correlation	Component 1	Component 2	Component 1	Component 2
True sample	0.10	0.74	0.10	0.74
Nearest-neighbour	-0.04	-0.69	0.03	0.13
Parameteric Model	0.16	0.04	0.01	0.81

Equations

$$\left(\begin{array}[]{cccc}\pagecolor{black!10}\bm{x}_{1}&\pagecolor{black!10}\bm{y}_{1}&\bm{z}_{1}\\ \pagecolor{black!10}\vdots&\pagecolor{black!10}\vdots&\vdots\\ \pagecolor{black!10}\bm{x}_{n_{A}}&\pagecolor{black!10}\bm{y}_{n_{A}}&\bm{z}_{n_{A}}\\ \pagecolor{black!10}\bm{x}_{n_{A}+1}&\bm{y}_{n_{A}+1}&\pagecolor{black!10}\bm{z}_{n_{A}+1}\\ \pagecolor{black!10}\vdots&\vdots&\pagecolor{black!10}\vdots\\ \pagecolor{black!10}\bm{x}_{n_{A}+n_{B}}&\bm{y}_{n_{A}+n_{B}}&\pagecolor{black!10}\bm{z}_{n_{A}+n_{B}}\end{array}\right)$$

$$\left[\begin{array}{c}{\mathbf{X}}_i\\ {\mathbf{Y}}_i\\ {\mathbf{Z}}_i\end{array}\right]$$

$${\mathbf{Y}}_i|{\mathbf{X}}_i={\mathbf{x}}_i$$

$${\mathbf{\beta }}_Y$$

$${\mathbf{\alpha }}_Y$$

$${\Omega }_Y$$

$${\mathbf{Z}}_i|\mathbf{X}={\mathbf{x}}_i$$

$${\mathbf{\beta }}_Z$$

$${\mathbf{\alpha }}_Z$$

$${\Omega }_Z$$

$$f({\mathbf{x}}_{\text{obs}},{\mathbf{y}}_{\text{obs}},{\mathbf{z}}_{\text{obs}};\mathbf{\eta })$$

$$=\prod _{i=1}^{n}f({\mathbf{x}}_i;{\mathbf{\eta }}_X)\prod _{j=1}^{n_A}f({\mathbf{y}}_j|{\mathbf{x}}_j;{\mathbf{\eta }}_Y)\prod _{k=n_A+1}^{n}f({\mathbf{z}}_k|{\mathbf{x}}_k;{\mathbf{\eta }}_Z)$$

$$=f({\mathbf{x}}_{\text{obs}};{\mathbf{\eta }}_X)f({\mathbf{y}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^A;{\mathbf{\eta }}_Y)f({\mathbf{z}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^B;{\mathbf{\eta }}_Z).$$

$$f({\mathbf{x}}_{\text{obs}},{\mathbf{y}}_{\text{obs}},{\mathbf{z}}_{\text{obs}};\mathbf{\eta })=f({\mathbf{x}}_{\text{obs}};{\mathbf{\eta }}_X)f({\mathbf{y}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^A;{\mathbf{\eta }}_Y)f({\mathbf{z}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^B;{\mathbf{\eta }}_Z)$$

$${\mathbf{B}}_A$$

$${\mathbf{Y}}_A$$

$${\mathbf{Z}}_B$$

$${\hat{\Gamma }}_A$$

$${\hat{\Gamma }}_B$$

$${\hat{\Omega }}_Y$$

$${\hat{\Omega }}_Z$$

$$f({\mathbf{x}}_i,{\mathbf{y}}_i,{\mathbf{z}}_i;\mathbf{\mu },\Sigma ,\mathbf{\delta })$$

$$\left[\begin{array}{c}{\mathbf{X}}_i\\ {\mathbf{Y}}_i\\ {\mathbf{Z}}_i\end{array}\right]$$

$$\left[\begin{array}{c}{\mathbf{X}}_i\\ {\mathbf{Y}}_i\\ {\mathbf{Z}}_i\end{array}\right]|U_i=u_i$$

$$\mathbb{E}[{\mathbf{Y}}_i|{\mathbf{X}}_i={\mathbf{x}}_i,U_i=u_i]$$

$$\mathbb{E}[{\mathbf{Z}}_i|{\mathbf{X}}_i={\mathbf{x}}_i,U_i=u_i]$$

$${\mathbf{Y}}_i|{\mathbf{X}}_i={\mathbf{x}}_i,U_i=u_i$$

$${\mathbf{Z}}_i|{\mathbf{X}}_i={\mathbf{x}}_i,U_i=u_i$$

$${\mathbf{Y}}_A$$

$${\mathbf{Z}}_B$$

$${\mathbf{B}}_A$$

$$f({\mathbf{x}}_i,{\mathbf{y}}_i,{\mathbf{z}}_i,u_i;\mathbf{\theta })$$

$$=f({\mathbf{x}}_i|u_i;\mathbf{\theta })f({\mathbf{y}}_i|{\mathbf{x}}_i,u_i;\mathbf{\theta })f({\mathbf{z}}_i|u_i;\mathbf{\theta })f(u_i).$$

$$f({\mathbf{x}}_{\text{obs}},{\mathbf{y}}_{\text{obs}},{\mathbf{z}}_{\text{obs}},\mathbf{u};\mathbf{\eta })=f({\mathbf{x}}_{\text{obs}}|\mathbf{u};{\mathbf{\eta }}_X)f({\mathbf{y}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^A,\mathbf{u};{\mathbf{\eta }}_Y)f({\mathbf{z}}_{\text{obs}}|{\mathbf{x}}_{\text{obs}}^B,\mathbf{u};{\mathbf{\eta }}_Z)f(\mathbf{u}).$$

$${\mathbf{X}}_i,{\mathbf{Y}}_i,{\mathbf{Z}}_i|S_i=h$$

$$S_i$$

$${\Sigma }_{YZ}^{(h)}$$

Full text

Introduction

The statistical matching problem is a data integration scenario where the set of available datasets shares only a strict subset of variables. The pattern of missingness is such that some dimensions are wholly unobserved in each dataset (D’Orazio et al., 2006). The simplest case involves two datasets, dataset A and dataset B, and three groups of variables, labelled $\bm{X}$, $\bm{Y}$ and $\bm{Z}$. We assume the data generating process is governed by a parametric model $f(\bm{x},\bm{y},\bm{z};\bm{\theta})$, for some parameter $\bm{\theta}\in\Omega\subset\mathbb{R}^{p}$. Dataset A contains $n_{A}$ independent observations from the marginal distribution $f(\bm{x},\bm{y};\bm{\theta})$ and dataset $B$ consists of $n_{B}$ independent observations from the marginal distribution $f(\bm{x},\bm{z};\bm{\theta})$. Table 1 represents the missing data pattern. A common objective in the statistical matching problem is to impute the missing observations in each dataset so complete data techniques can be used in downstream analyses (Rässler, 2002). The statistical matching problem occurs in flow cytometry analysis due to technological limitations, and there has been recent interest in developing appropriate statistical methods for integrating cytometry datasets (Pedreira et al., 2008; O’Neill et al., 2015). Most previous work on the statistical matching problem for continuous data assumes a multivariate normal model. Flow cytometry data typically have characteristics that make a multivariate Gaussian model inappropriate. Observations are from a number of heterogeneous cell subpopulations, and observations within a cell subtype often exhibit skewness and heavy tails. Finite mixtures of skew-normal and skew-$t$ distributions have shown excellent performance in the analysis of flow cytometry data (Pyne et al., 2009; Lee and McLachlan, 2016). We study theoretical aspects of the statistical matching of non-Gaussian data and propose new computationally efficient algorithms for the statistical matching of skew-normal data, mixtures of Gaussians and mixtures of skew-normal distributions.

A central issue in the statistical matching problem is that the lack of joint observations on the $\bm{Y}$ and $\bm{Z}$ variables renders most statistical models nonidentifiable. For example, if $f(\bm{x},\bm{y},\bm{z};\bm{\theta})$ is a multivariate normal distribution, the $\bm{Y},\bm{Z}$ covariance parameters are not identifiable. In a parametric framework, it is common to adopt restrictions on the parameters so that the generative model remains identifiable. Under identification constraints it is possible to impute the missing data using a parametric modelling strategy. Let $\widehat{\bm{\theta}}$ denote the maximum likelihood estimate of $\bm{\theta}$. The missing $\bm{Z}$ values in dataset A can be sampled from $f(\bm{z}|\bm{x},\bm{y};\widehat{\bm{\theta}})$, and the missing $\bm{Y}$ values in dataset B can be sampled from $f(\bm{y}|\bm{x},\bm{z};\widehat{\bm{\theta}})$.

Parametric strategies for Gaussian data typically involve reparameterising the model so that the likelihood can then be factored into a mathematically tractable form (D’Orazio et al., 2006). Reparameterisation and factorisation are powerful strategies in missing data problems, as it can greatly simplify maximum likelihood estimation (Little and Rubin, 2002, Chapter 7). The reparameterisation for the statistical matching of Gaussian data involves expressing the observed-data likelihood in terms of a multiple response regression model. We show how the regression model specification can be extended to cover skew-normal distributions, mixtures of Gaussians, and mixtures of skew-normal distributions through appropriate data augmentation schemes. Factorisation of the complete-data log likelihood allows the EM algorithm (Dempster et al., 1977) to be used for computationally efficient parameter estimation.

An alternative method for imputation is nearest-neighbour matching, a nonparametric method based on the assumption that the $\bm{Y}$ and $\bm{Z}$ variables are conditionally independent given the $\bm{X}$ variables. Nearest-neighbour imputation has been widely applied in the statistical matching problem due to its ease of implementation and generality (D’Orazio et al., 2006). In particular, nearest-neighbour matching has been advocated for the file matching of flow cytometry data (Pedreira et al., 2008). An important finding is that the $\bm{Y}$ and $\bm{Z}$ variables are likely to be conditionally dependent given the $\bm{X}$ variables if the generative model is a mixture model or a skewed distribution. The violation of the conditional independence assumption can lead to nearest-neighbour imputation exhibiting undesirable behaviour.

The conditional independence assumption can be adopted for any data generating process, and is the dominant approach to resolve the pathological missing data pattern in the statistical matching problem. This conditional independence assumption is difficult to justify for a range of latent variable models, and model based approaches can impose more coherent constraints. We compare nearest-neighbour matching to model based imputation on a number of real datasets. The model based approach showed a large improvement over the nearest-neighbour method, however there can be issues when the identification constraints are violated. Parametric imputation strategies are a flexible and computationally feasible option for statistical matching.

Background

Notation

We assume $n_{A}$ observations in dataset A and $n_{B}$ observations in dataset B, for a total of $n=n_{A}+n_{B}$ records. These datasets are combined to form a single data matrix, where rows 1 to $n_{A}$ contain the records from dataset $A$ on the $\bm{X}$ and $\bm{Y}$ variables and rows $n_{B}$ rows contain the observations from dataset $B$ on the $\bm{X}$ and $\bm{Z}$ variables. Let $d_{X}$, $d_{Y}$ and $d_{Z}$ give the dimension of the $\bm{X},\bm{Y}$ and $\bm{Z}$ variables respectively and set $d=d_{X}+d_{Y}+d_{Z}$.

We assume the data generating process consists of $n$ independently and identically distributed observations from the parametric model $(\bm{X}_{i}^{\mathsf{T}},\bm{Y}_{i}^{\mathsf{T}},\bm{Z}_{i}^{\mathsf{T}})^{\mathsf{T}}\sim f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\theta})$ for $i=1,\ldots,n$. Let $\bm{x}_{i}$, $\bm{y}_{i}$ and $\bm{z}_{i}$ represent the realised values for the $\bm{X}_{i}$, $\bm{Y}_{i}$ and $\bm{Z}_{i}$ random vectors respectively for $i=1,\ldots,n$. Observations $\bm{z}_{i}$ are missing for $i=1,\ldots,n_{A}$. Observations $\bm{y}_{i}$ are missing for $j=n_{A}+1,n_{A}+2,\ldots,n_{A}+n_{B}$. The joint $n\times d$ data matrix of interest is represented in the display (1). Observed values are shaded, with white cells indicating missing values in the joint data matrix. Let $\bm{x}_{\text{obs}}=(\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{n})$ represent the observed $\bm{X}$ values over the $n$ records. Let $\bm{x}_{\text{obs}}^{A}=(\bm{x}_{1},\ldots,\bm{x}_{n_{A}})$ represent the $n_{A}$ observed $\bm{X}$ values on dataset $A$.

[TABLE]

Let $\bm{x}_{\text{obs}}^{B}=(\bm{x}_{n_{A}+1},\ldots,\bm{x}_{n_{A}+n_{B}})$ denote the $n_{B}$ observed $\bm{X}$ values in dataset ${B}$. Let $\bm{y}_{\text{obs}}=(\bm{y}_{1},\ldots,\bm{y}_{n_{A}})$ represent the $n_{A}$ observed values $\bm{Y}$ variables in dataset $A$, and let $\bm{z}_{\text{obs}}=(\bm{z}_{n_{A}+1},\ldots,\bm{z}_{n_{A}+n_{B}})$ denote the $n_{B}$ observed $\bm{Z}$ values in dataset $B$. Likewise, let $\bm{z}_{\text{mis}}=(\bm{z}_{1},\ldots,\bm{z}_{n_{A}})$ represent the $n_{A}$ missing $\bm{Z}$ values in dataset $A$, and let $\bm{y}_{\text{mis}}=(\bm{y}_{n_{A}+1},\ldots,\bm{y}_{n_{A}+n_{B}})$ represent the $n_{B}$ missing $\bm{Y}$ values in dataset $B$.

Gaussian data

The statistical matching problem has been explored in depth under the assumption of a multivariate Gaussian generative model. Suppose that for $i=1,\ldots,n$ we have observations from the multivariate normal distribution

[TABLE]

where the mean and covariance parameters have been partitioned in an obvious fashion. With the missing data pattern represented in (1), the only non-identifiable parameter is $\Sigma_{YZ}$. The most common identifiability constraint is that $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$.

Parameter estimation under the identification constraint $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$ has a long history in the literature (Lord, 1955; Anderson, 1957; Moriarity and Scheuren, 2003). As discussed in D’Orazio et al. (2006, Chapter 2), the multivariate normal statistical matching problem has a useful connection to linear regression modelling that can be used to obtain closed form maximum likelihood estimates. The conditional distribution of $\bm{Y}_{i}$ given $\bm{X}_{i}=\bm{x}_{i}$ can be represented as a regression model

[TABLE]

for $i=1,\ldots,n$ where

[TABLE]

The same holds for the conditional distribution of $\bm{Z}_{i}$ given $\bm{X}_{i}=\bm{x}_{i}$

[TABLE]

for $i=1,\ldots,n$ where

[TABLE]

Now let $\bm{\eta}=(\bm{\eta}_{X},\bm{\eta}_{Y},\bm{\eta}_{Z})$, where $\bm{\eta}_{X}=(\bm{\mu}_{X},\Sigma_{XX}),\bm{\eta}_{Y}=(\bm{\alpha}_{Y},\bm{\beta}_{Y},\Omega_{Y})$ and $\bm{\eta}_{Z}=(\bm{\alpha}_{Z},\bm{\beta}_{Z},\Omega_{Z})$. The benefit of the regression parameterisation is that that the likelihood for each observation factors into three components with distinct parameter blocks. We can write $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\eta})=f(\bm{x}_{i};\bm{\eta}_{X})f(\bm{y}_{i}|\bm{x}_{i};\bm{\eta}_{Y})f(\bm{z}_{i}|\bm{x}_{i};\bm{\eta}_{Z})$. The observed-data likelihood can be expressed as a product of three terms,

[TABLE]

The first likelihood block $f(\bm{x}_{\text{obs}};\bm{\eta}_{X})$ involves the $n$ observed $X$ variables. The second likelihood block $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A};\bm{\eta}_{Y})$ involves the conditional likelihood for dataset A and the third likelihood block $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}}^{B};\bm{\eta}_{Z})$ involves the conditional likelihood for dataset B. Maximisation of the observed-data likelihood

[TABLE]

is straightforward, as we now have three separate maximisation problems over $\bm{\eta}_{X}$, $\bm{\eta}_{Y}$ and $\bm{\eta}_{Z}$. Both datasets are used to estimate $\bm{\eta}_{X}$ using the fully observed $\bm{X}$ variables. Dataset A is used to estimate $\bm{\eta}_{Y}$ and dataset B is used to estimate $\bm{\eta}_{Z}$.

Using the regression specifications (3) and (7), we can express the likelihoods $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A};\bm{\eta}_{Y})$ and $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}}^{B};\bm{\eta}_{Z})$ as multiple response regression models. Recall the data structure in display (1). Let $\bm{Y}_{A}$ give the $n_{A}\times d_{Y}$ matrix of responses from dataset A. Row $i$ in $\bm{Y}_{A}$ is given by $\bm{y}_{i}^{\mathsf{T}}$ for $i=1,\ldots,n_{A}$. Let $\bm{Z}_{B}$ give the $n_{B}\times d_{Z}$ from dataset B. Row $i$ in $\bm{Z}_{B}$ is given by $\bm{z}_{i}^{\mathsf{T}}$ for $i=n_{A}+1,\ldots,n_{A}+n_{B}$. Let $\bm{B}_{A}$ represent the design matrix for dataset A and let $\bm{B}_{B}$ give the design matrix for dataset $B$. Specifically,

[TABLE]

Let $\Gamma_{A}$ and $\Gamma_{B}$ contain the regression coefficients for (3) and (7) respectively: $\Gamma_{A}=[\bm{\alpha}_{Y}\ \bm{\beta}_{Y}]^{\mathsf{T}}$ and $\Gamma_{B}=[\bm{\alpha}_{Z}\ \bm{\beta}_{Z}]^{\mathsf{T}}$. We will use the matrix normal distribution introduced by Dawid (1981) to specify the regression models compactly. A random $n\times p$ matrix $\bm{M}$ is said to have the matrix normal distribution $MN(\bm{I}_{n},\Sigma)$ if each row is a draw from a $p$-variate normal distribution $N(\bm{0},\Sigma)$. The conditional regression models for each dataset can be written as

[TABLE]

where $\bm{\epsilon}_{A}\sim MN(\bm{I}_{n_{A}},\Omega_{Y})$ and $\bm{\epsilon}_{B}\sim MN(\bm{I}_{n_{B}},\Omega_{Z})$. Equation (16) is a representation of the likelihood factor $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A};\bm{\eta}_{Y})$ and equation (17) is a representation of the likelihood factor $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}^{B}};\bm{\eta}_{Z})$ . The sufficient statistics for the regressions are $\bm{Y}_{A}^{\mathsf{T}}\bm{Y}_{A},\bm{Z}_{B}^{\mathsf{T}}\bm{Z}_{B},\bm{B}_{A}^{\mathsf{T}}\bm{Y}_{A},\bm{B}_{B}^{\mathsf{T}}\bm{Z}_{B},\bm{B}_{A}^{\mathsf{T}}\bm{B}_{A}$ and $\bm{B}_{B}^{\mathsf{T}}\bm{B}_{B}$. The maximum likelihood estimators of the regression parameters $\Gamma_{A}=[\bm{\alpha}_{Y}\ \bm{\beta}_{Y}]^{\mathsf{T}}$ and $\Gamma_{B}=[\bm{\alpha}_{Z}\ \bm{\beta}_{Z}]^{\mathsf{T}}$ are given by

[TABLE]

The maximum likelihood estimators of the error covariance matrices $\Omega_{Y}$ and $\Omega_{Z}$ are given by

[TABLE]

These results follow from general results on multiple response regression models, for example see Rencher and Christensen (2012, Chapter 10). The maximum likelihood estimators of the $\bm{\eta}_{X}$ parameters are given by the sample mean and the sample covariance of the $n$ observed $\bm{X}$ values $\widehat{\bm{\mu}}_{X}=\dfrac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}$, $\widehat{\Sigma}_{XX}=\dfrac{1}{n}\sum_{i=1}^{n}(\bm{x}_{i}-\widehat{\bm{\mu}}_{X})(\bm{x}_{i}-\widehat{\bm{\mu}}_{X})^{\mathsf{T}}$. Inverting the transformations defined in equations (5) to (10) gives the maximum likelihood estimates of the original parameters $\bm{\mu}$ and $\bm{\Sigma}$. We thus obtain closed form maximum likelihood estimators under the identification restriction $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$.

Extensions

Overview

A parametric approach in the statistical matching problem requires the identification of suitable identifiability constraints, and a feasible procedure for maximum likelihood estimation. One approach is to introduce the missing data $\bm{y}_{\text{obs}}$, $\bm{z}_{\text{obs}}$ into the model as latent variables. In general, we can form the complete-data likelihood $f(\bm{x}_{\text{obs}},\bm{y}_{\text{obs}},\bm{z}_{\text{obs}},\bm{y}_{\text{mis}},\bm{z}_{\text{mis}};\bm{\theta})$ and define appropriate EM iterations for parameter estimation. Lee et al. (2011) consider the statistical matching of flow cytometry data with mixtures of PCA models and take this approach. Introduction of the missing observations $\bm{y}_{\text{mis}}$ and $\bm{z}_{\text{mis}}$ into the complete data log-likelihood can be computationally demanding on flow cytometry datasets with a large number of observations. Due to the non-identifiability of the model, the fitted mixture model is only used to assign cluster labels to the observations. Lee et al. propose to use nearest-neighbour matching within each group of labelled points.

We propose a different method, where we introduce appropriate parameter constraints so that $\bm{y}_{\text{mis}}$ and $\bm{z}_{\text{mis}}$ do not need to be included in the complete-data log likelihood. This is to obtain a more computationally efficient EM algorithm, and to avoid the identifiability issues with fitting a completely unconstrained model. Secondly, we propose to impute the missing data using the constrained fitted model. The strategy involves extending the linear regression connection that was used for Gaussian data.

Skew-normal

There are many forms of the skew-normal distribution, we will work with the same multivariate version as in Pyne et al. (2009). For compact notation, let the vector $\bm{w}_{i}$ represent the vector of joint observations $\bm{w}_{i}=(\bm{x}_{i}^{\mathsf{T}},\bm{y}_{i}^{\mathsf{T}},\bm{z}_{i}^{\mathsf{T}})^{\mathsf{T}}$. The density of the skew-normal distribution is

[TABLE]

where $\Lambda=\Sigma+\bm{\delta}\bm{\delta}^{\mathsf{T}}$ and $\alpha^{\mathsf{T}}=\bm{\delta}^{\mathsf{T}}\Lambda^{-1}/(1-\bm{\delta}^{\mathsf{T}}\Lambda^{-1}\bm{\delta})$. Let $TN(\mu,\sigma^{2},a)$ denote a lower truncated normal distribution where $\mu$ and $\sigma^{2}$ give the mean and variance of the underling normal distribution and $a$ gives the lower truncation bound. We say $U\sim TN(\mu,\sigma^{2},a)$ if $U\overset{d}{=}[W|W>a]$ where $W\sim N(\mu,\sigma^{2})$. The model (18) has a hierarchical representation,

[TABLE]

where ${U}_{i}\sim TN(0,1,0)$ and $\bm{V}_{i}\sim N(\bm{0},\Sigma)$, where $\bm{\mu}$ and $\Sigma$ are partitioned as in (2). The only non-identifiable parameter is $\Sigma_{YZ}$. The regression model specification also enables efficient inference for the skew-normal distribution under the identification restriction $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. The statistical matching problem for skew-normal data can be expressed in terms of a multiple response regression model with latent unobserved variables in the conditional mean function. Using the hierarchical model (19), conditional on knowing $U_{i}=u_{i}$, it holds that,

[TABLE]

where $\bm{V}_{i}=N(\bm{0},\Sigma)$, for $i=1,\ldots,n$. Let $\bm{\beta}_{Y}=\Sigma_{YX}\Sigma_{XX}^{-1}$ and $\bm{\beta}_{Z}=\Sigma_{ZX}\Sigma_{XX}^{-1}$. The conditional mean of $\bm{Y}_{i}$ given $\bm{X}_{i}$ and ${U}_{i}$, and the conditional mean of $\bm{Z}_{i}$ given $\bm{X}_{i}$ and $U_{i}$ are given respectively by:

[TABLE]

We have that $\text{var}(\bm{Y}_{i}|\bm{X}_{i},U_{i})=\Sigma_{YY}-{\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XY}$ and $\text{var}(\bm{Z}_{i}|\bm{X}_{i},U_{i})=\Sigma_{ZZ}-{\Sigma}_{ZX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. The augmented likelihood for dataset A and the augmented likelihood for dataset B can be expressed as conditional regression models. Let $\bm{\alpha}_{Y}=\bm{\mu}_{Y}-\bm{\beta}_{Y}\bm{\mu}_{X},\Omega_{Y}=\Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}$ and $\bm{\lambda}_{Y}=\bm{\delta}_{Y}-\bm{\beta}_{Y}\bm{\delta}_{X}$. Collecting terms in (21), the distribution of $\bm{Y}_{i}$ given $\bm{X}_{i}$ and the latent scaling variable $U_{i}$ can be represented as the regression model:

[TABLE]

Let $\bm{\beta}_{Z}=\Sigma_{ZX}\Sigma_{XX}^{-1},\bm{\alpha}_{Z}=\bm{\mu}_{Z}-\bm{\beta}_{Z}\bm{\mu}_{X},\Omega_{Z}=\Sigma_{ZZ}-\Sigma_{ZX}\Sigma_{XX}^{-1}\Sigma_{XZ}$ and $\bm{\lambda}_{Z}=\bm{\delta}_{Z}-\bm{\beta}_{Z}\bm{\delta}_{X}$. Similarly, collecting terms in (22) the conditional distribution of $\bm{Z}_{i}$ given $\bm{X}_{i}$ and the latent scaling variable $U_{i}$ can be represented as a regression model:

[TABLE]

Let $\Gamma_{A}$ and $\Gamma_{B}$ contain the regression parameters for (23) and (24) respectively: $\Gamma_{A}=[\bm{\alpha}_{Y}\ \bm{\lambda}_{Y}\ \bm{\beta}_{Y}]^{\mathsf{T}}$ and $\Gamma_{B}=[\bm{\alpha}_{Z}\ \bm{\lambda}_{Z}\ \bm{\beta}_{Z}]^{\mathsf{T}}$. The regressions for dataset A and dataset B can be written as

[TABLE]

where $\bm{\epsilon}_{A}\sim MN(\bm{I}_{n_{A}},\Omega_{Y})$ and $\bm{\epsilon}_{B}\sim MN(\bm{I}_{n_{B}},\Omega_{Z})$. The design matrices for the regressions now include the latent $U_{i}$ terms. The complete-data design matrices are given by

[TABLE]

The complete-data sufficient statistics are again $\bm{Y}_{A}^{\mathsf{T}}\bm{Y}_{A},\bm{Z}_{B}^{\mathsf{T}}\bm{Z}_{B},\bm{B}_{A}^{\mathsf{T}}\bm{Y}_{A},\bm{B}_{B}^{\mathsf{T}}\bm{Z}_{B},\bm{B}_{A}^{\mathsf{T}}\bm{B}_{A}$ and $\bm{B}_{B}^{\mathsf{T}}\bm{B}_{B}$. Suppose that we impose the identification restriction $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. This corresponds to a conditional independence assumption on the components of the latent random variable $\bm{V}$ in the hierarchical skew-normal model (19). Under the the identification restriction, the $\bm{Y}$ and $\bm{Z}$ variables are conditionally independent given $\bm{X}$ and the latent $U$. The augmented data likelihood for a single observation then has the structure

[TABLE]

Using this property we can show the complete-data likelihood has a similar structure to the Gaussian likelihood (14). Let $\bm{\eta}=(\bm{\eta}_{X},\bm{\eta}_{Y},\bm{\eta}_{Z})$ represent the parameters for each likelihood block. We have $\bm{\eta}_{X}=(\bm{\mu}_{X},\Sigma_{XX},\bm{\delta}_{X})$, $\bm{\eta}_{Y}=(\bm{\alpha}_{Y},\bm{\beta}_{Y},\Omega_{Y},\bm{\lambda}_{Y})$ and $\bm{\eta}_{Z}=(\bm{\alpha}_{Z},\bm{\beta}_{Z},\Omega_{Z},\bm{\lambda}_{Z})$. Let $\bm{u}$ contain the latent skewing variable for each observation, so $\bm{u}=(u_{1},\ldots,u_{n})$. Using the conditional regression models, the complete-data likelihood can be expressed as

[TABLE]

Equation (25) is a representation of the likelihood factor $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A},\bm{u};\bm{\eta}_{Y})$ and equation (26) is a representation of the likelihood factor $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}}^{B},\bm{u};\bm{\eta}_{Z})$. The E-step involves calculating the expected value of the sufficient statistics given the current parameter estimates, and can be carried out using the results in Pyne et al. (2009). The M-step is simplified due to choice of parameterisation for the regression model. There are three separate maximisation tasks over $\bm{\eta}_{X}$, $\bm{\eta}_{Y}$ and $\bm{\eta}_{Z}$ that have closed form solutions. The key point is that we do not have to introduce the missing observations $\bm{y}_{\text{mis}}$ and $\bm{z}_{\text{mis}}$ to obtain a workable EM algorithm under the identification restriction $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. The algorithm is given in Web Appendix A.

Mixture of Gaussians

We can also obtain an efficient EM algorithm for statistical matching of Gaussian mixture data using the conditional regression specifications and appropriate identification restrictions. The matching problem for Gaussian mixture data can be represented as a mixture of regression problem (Quandt, 1972). Recall the hierarchical model for the mixture model discussed in Section 4.2. For independent observations $i=1,\ldots,n$:

[TABLE]

Here $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\theta}_{h})$ denotes a multivariate normal density with parameters $\bm{\theta}_{h}=(\bm{\mu}^{(h)},\Sigma^{(h)})$. The vector $\bm{\mu}^{(h)}$ gives the mean vector for component $h$ and the matrix $\Sigma^{(h)}$ represent the covariance matrix of the $h$th component distribution for $h=1,\ldots,g$. We can partition each component mean $\bm{\mu}^{(h)}$ and each component covariance matrix $\Sigma^{(h)}$ as in (2). An identification restriction is then

[TABLE]

Consider the augmented likelihood for a single observation given the latent cluster indicator $S_{i}$. Conditional on the latent cluster indicator $S_{i}$, we can again use the regression model specification from Section 2.2. Specifically,

[TABLE]

where $\bm{\beta}_{Y}^{(h)}=\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}$, $\bm{\alpha}_{Y}^{(h)}=\bm{\mu}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\mu}_{X}^{(h)}$ and $\Omega_{Y}^{(h)}=\Sigma_{YY}^{(h)}-\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XY}^{(h)}$. The conditional distribution of $\bm{Z}_{i}$ given $\bm{X}_{i}$ obeys a similar regression equation:

[TABLE]

where $\bm{\beta}_{Z}^{(h)}=\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}$, $\bm{\alpha}_{Z}^{(h)}=\bm{\mu}_{Z}^{(h)}-\bm{\beta}_{Z}\bm{\mu}_{X}$ and $\Omega_{Z}=\Sigma_{ZZ}-\Sigma_{ZX}\Sigma_{XX}^{-1}\Sigma_{XZ}$. To describe the complete-data log likelihood first let $\bm{s}$ represent the vectors of latent cluster indicators, so $\bm{s}=(s_{1},\ldots,s_{n})$, where $s_{i}\in\left\{1,\ldots,g\right\}$ for $i=1,\ldots,n$. Let $\bm{\pi}=(\pi_{1},\ldots,\pi_{g})$ be a vector containing the mixing proportions. We can show that complete-data likelihood has a similar structure to what we obtained in for the Gaussian likelihood (13). For Gaussian mixture data, the complete-data likelihood has the structure

[TABLE]

The complete-data likelihood is the product of four components. The first likelihood block $f(\bm{x}_{\text{obs}}|\bm{s};\bm{\eta}_{X})$ is a standard Gaussian likelihood. The second and third likelihood blocks $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A},\bm{s};\bm{\eta}_{Y})$ and $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}}^{B},\bm{s};\bm{\eta}_{Z})$ correspond to regression likelihoods. The final likelihood contribution from $f(\bm{s};\bm{\pi})$ is the usual multinomial likelihood seen in a finite mixture. The E-step and M-step for $f(\bm{x}_{\text{obs}}|\bm{s};\bm{\eta}_{X})$ follow from standard results on Gaussian mixture models. The E-step and M-steps for $f(\bm{y}_{\text{obs}}|\bm{x}_{\text{obs}}^{A},\bm{s};\bm{\eta}_{Y})$ and $f(\bm{z}_{\text{obs}}|\bm{x}_{\text{obs}}^{B},\bm{s};\bm{\eta}_{Z})$ follow from results on mixtures of regression models given in Jones and McLachlan (1992). The details are given in Web Appendix B. We avoid introducing the missing observations $\bm{y}_{\text{mis}}$ and $\bm{z}_{\text{mis}}$ into the complete-data log likelihood by using the component-wise restrictions (30).

Mixtures of skew-normal

The statistical matching of data from a mixture of skew-normal distributions as a mixture of regressions problem, with latent variables in the design matrix. For independent observations $i=1,\ldots,n$ we have that

[TABLE]

Here $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\theta}_{h})$ denotes a skew-normal density with parameters $\bm{\theta}_{h}=(\bm{\mu}^{(h)},\Sigma^{(h)},\bm{\delta}_{h})$. The vector $\bm{\mu}^{(h)}$ gives the mean vector for component $h$, the matrix $\Sigma^{(h)}$ represents the scale matrix for component $h$ and $\bm{\delta}^{(h)}$ gives the skewness vector for component $h$ for $h=1,\ldots,g$. We can partition each component mean $\bm{\mu}^{(h)}$ and skewness vector $\bm{\delta}^{(h)}$ as in (19). Each component scale matrix $\Sigma^{(h)}$ can be partitioned as in (2). An identification restriction is then

[TABLE]

Introducing the latent $U_{i}$ and $S_{i}$ into the augmented likelihood, the complete-data model can be expressed in terms of conditional regression specifications. Let $\bm{\beta}_{Y}^{(h)}=\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1},\bm{\alpha}_{Y}^{(h)}=\bm{\mu}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\mu}_{X}^{(h)},\Omega_{Y}^{(h)}=\Sigma_{YY}^{(h)}-\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XY}^{(h)}$ and $\bm{\lambda}_{Y}^{(h)}=\bm{\delta}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\delta}_{X}^{(h)}$.

[TABLE]

Let $\bm{\beta}_{Z}^{(h)}=\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1},\bm{\alpha}_{Z}^{(h)}=\bm{\mu}_{Z}^{(h)}-\bm{\beta}_{Z}^{(h)}\bm{\mu}_{X}^{(h)},\Omega_{Z}^{(h)}=\Sigma_{ZZ}^{(h)}-\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XZ}^{(h)}$ and $\bm{\lambda}_{Z}^{(h)}=\bm{\delta}_{Z}^{(h)}-\bm{\beta}_{Z}^{(h)}\bm{\delta}_{X}^{(h)}$. Similarly, collecting terms in (22) the conditional distribution of $\bm{Z}_{i}$ given $\bm{X}_{i}$ and the latent scaling variable $U_{i}$ can be represented as a regression model

[TABLE]

The complete-data likelihood has the structure

[TABLE]

The first likelihood block $f(\bm{x}_{\text{obs}}|\bm{s},\bm{u};\bm{\eta}_{X})$ corresponds to a skew-normal likelihood. The second and third likelihood blocks are regression likelihoods with the latent $u$ variables in the design matrices. We can again define an EM algorithm for parameter estimation under the identification constraint (33). The E-step is carried out using existing results for the skew-normal distribution. The M-step is again simplified by the choice of parameterisation. We have three separate maximisation tasks over $\bm{\eta}_{X}$, $\bm{\eta}_{Y}$ and $\bm{\eta}_{Z}$. The algorithm is given in full in Web Appendix C. Once again we do not introduce the missing observations $\bm{y}_{\text{mis}}$ and $\bm{z}_{\text{mis}}$ into complete-data likelihood by using the identification restriction (33).

Nearest-neighbour matching

Nearest-neighbour imputation is a popular alternative to parametric methods in data fusion problems (Aluja-Banet et al., 2007; Saporta, 2002). Nearest-neighbour methods rely on the conditional independence assumption, that is $\bm{Y}$ and $\bm{Z}$ are conditionally independent given $\bm{X}$. The nearest-neighbour method matches observations in dataset A and dataset B based on the Euclidean distance measured using the common $\bm{X}$ dimensions. Missing values are imputed by taking values from the nearest-neighbour in the donor set. More formally, the missing $\bm{Z}$ values for observation $i$ in dataset $A$ are set as

[TABLE]

for $i=1,\ldots,n_{A}$. Similarly, the missing $\bm{Y}$ values for observation $i$ in dataset B are set as

[TABLE]

for $i=n_{A}+1,\ldots,n_{A}+n_{B}$. As discussed by Rässler (2002), each of the observations in imputed dataset can be viewed as exchangeable draws from some distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})$. The distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})$ is useful to characterise the behaviour of the nearest-neighbour method. Marella et al. (2008) consider the asymptotic form of the nearest-neighbour imputation distribution and show that as the size of the donor set tends to infinity,

[TABLE]

In equation (37) $f(\bm{x}_{i};\bm{\theta}),f(\bm{y}_{i}|\bm{x}_{i};\bm{\theta})$ and $f(\bm{z}_{i}|\bm{x}_{i};\bm{\theta})$ are the marginal and conditional distributions from the true generative model $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\theta})$. The nearest-neighbour strategy will only produce the correct joint distribution if $\bm{Y}$ and $\bm{Z}$ are conditionally independent given $\bm{X}$. For multivariate-normal data, the asymptotic form of the nearest-neighbour imputation scheme is

[TABLE]

The distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})$ is equivalent to the distribution recovered by the maximum likelihood parametric approach under the identification constraint $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. The asymptotic equivalence between nearest-neighbour matching and a parametric methods may not necessarily hold for non-Gaussian data.

It can be difficult to justify the conditional independence assumption for non-Gaussian data as there may be no parameter $\bm{\theta}^{*}$ such that $f(\bm{x},\bm{y},\bm{z};\bm{\theta}^{*})=f(\bm{x};\bm{\theta}^{*})f(\bm{y}|\bm{x};\bm{\theta}^{*})f(\bm{z}|\bm{x};\bm{\theta}^{*})$. This problem is likely to arise when there is some latent structure in the generative model that links the $\bm{Y}$ and $\bm{Z}$ variables. In these situations, the nearest-neighbour imputation distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})=f(\bm{x};\bm{\theta})f(\bm{y}|\bm{x}_{i};\bm{\theta})f(\bm{z}|\bm{x};\bm{\theta})$ may not be a good approximation for the true generative model $f(\bm{x},\bm{y},\bm{z};\bm{\theta})$. Even in situations where the generative model is identifiable, the nearest-neighbour scheme may fail to produce statistically sound imputations. We examine this issue for skew-normal distribution and Gaussian mixtures.

Skew-normal

In general, the imputation of skewed data can be challenging relative to the multivariate normal case (Sterne et al., 2009), and the statistical matching of skew-normal data presents some difficulties. Assuming that all variables have some skewness, so all elements of $\bm{\delta}$ are nonzero, it is not possible for $\bm{Y}$ and $\bm{Z}$ to be conditionally independent given $\bm{X}$ (Azzalini and Capitanio, 1999). Contrary to the normal case, there is no parameter constraint on $\bm{\mu},\Sigma,\bm{\delta}$ such that $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\mu},\Sigma,\bm{\delta})=f(\bm{x}_{i};\bm{\mu},\Sigma,\bm{\delta})f(\bm{y}_{i}|\bm{x}_{i};\bm{\mu},\Sigma,\bm{\delta})f(\bm{z}_{i}|\bm{x}_{i};\bm{\mu},\Sigma,\bm{\delta})$. The absence of a conditional independence constraint poses problems for nearest-neighbour imputation. To illustrate, we generated data from a three-dimensional skew-normal distribution. The parameters were set as $\bm{\mu}=\bm{0}$, $\Sigma=\bm{I}$ and $\bm{\delta}=(1,3,5)^{\mathsf{T}}$. The first, second and third dimensions were labelled as $\bm{X}$, $\bm{Y}$ and $\bm{Z}$ respectively. We generated two datasets with $n_{A}=n_{B}=5000$. Panels (a) and (b) in Figure 1 compare the true joint $(\bm{Y},\bm{Z})$ distribution to the nearest-neighbour asymptotic distribution. The nearest-neighbour distribution shows a much weaker linear association between the $\bm{Y}$ and $\bm{Z}$ variables compared to the true distribution. The contours in panel (b) are more square shaped than the skewed ellipse in panel (a). We applied nearest-neighbour imputation to the simulated dataset to check the correspondence with the theoretical distribution in (b). Panel (c) shows a smoothed scatter plot of the imputed data along with contours from a kernel density estimate as red dashed lines. The contours of the density estimate closely resemble the theoretical contours in (b), and are again more square shaped than the contours of the true $(\bm{Y},\bm{Z})$ distribution in (a).

Finite mixture models

Finite mixture models are well suited for the statistical analysis of flow cytometry data. To describe the model, let $\bm{\Psi}$ contain the mixing proportions $\pi_{1},\ldots,\pi_{g}$ and parameters $\bm{\theta}_{1},\ldots,\bm{\theta}_{g}$ for each component. The $j$th component distribution can be represented as $f(\bm{x},\bm{y},\bm{z};\bm{\theta}_{j})$ for $j=1,\ldots,g$. The distribution function of the $g$-component mixture model is given by

[TABLE]

It will be helpful to conceptualise the mixture model in terms of a latent cluster indicator $S_{i}\sim\text{Categorical}(\pi_{1},\ldots,\pi_{g})$ for each observation $i=1,\ldots,n$. The generative model has the hierarchical representation:

[TABLE]

It is unlikely for the conditional independence assumption to be appropriate when the generative model is a finite mixture. Given a mixture model with well-separated components, it is not possible to find an appropriate restrictions on the mixture parameters $\bm{\Psi}$ such that $f(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i};\bm{\Psi})=f(\bm{x}_{i};\bm{\Psi})f(\bm{y}_{i}|\bm{x}_{i};\bm{\Psi})f(\bm{z}_{i}|\bm{x}_{i};\bm{\Psi})$. Conditional on $\bm{X}$, the remaining $\bm{Y}$ and $\bm{Z}$ variables are almost certainly dependent due to the latent cluster indicator $S$. The violation of the conditional independence assumption for finite mixture models means that the nearest-neighbour imputation may have undesirable behaviour. The tendency for nearest-neighbour imputation to produce spurious clusters has been demonstrated empirically in statistical matching problems in flow cytometry analysis (Lee et al., 2011; O’Neill et al., 2015).

From the analysis of the asymptotic model, spurious clusters emerge when the nearest-neighbour match originates from a different mixture component to the query point $\bm{x}_{i}$. The probability of observing spurious clusters is related to how informative the common $\bm{X}$ variables are for classification. If the clusters are poorly separated using the $\bm{X}$ variables, then there is a high probability of mismatching observations across clusters and giving improper imputations.

To illustrate, we generated data from a equally weighted two component Gaussian mixture model. The first, second and third dimensions were labelled as the $\bm{X}$, $\bm{Y}$ and $\bm{Z}$ variables respectively. Component 1 had a mean of $(-0.1,0,0)^{\mathsf{T}}$ and component 2 had a mean of $(0.1,1,1)^{\mathsf{T}}$. Each component had covariance matrix $0.01\bm{I}$. The components are well separated in the marginal $\bm{Y}$ and $\bm{Z}$ dimensions and there is a large amount of overlap in the $\bm{X}$ dimension. We generated $n_{A}=n_{B}=500$ observations in each dataset. Panels (a) and (b) in Figure 2 compare the true distribution on the $\bm{Y}$ and $\bm{Z}$ variables to the asymptotic nearest neighbour distribution. The nearest-neighbour distribution shows four clusters in the joint distribution when there should only be two. This is because the $\bm{X}$ variables contain limited information for clustering. Panel (c) shows the results of applying nearest-neighbour matching to the simulated dataset. The introduction of two spurious clusters is consistent with the predicted behaviour from panel (b).

Examples

Skew-normal

We return to the skew-normal example in Section 4.1. We simulated another two datasets with $n_{A}=n_{B}=500$. The parameters were again set as $\bm{\mu}=\bm{0}$, $\Sigma=\bm{I}$ and $\bm{\delta}=(1,3,5)^{\mathsf{T}}$. The first, second and third dimensions were labelled as $\bm{X}$, $\bm{Y}$ and $\bm{Z}$ respectively. We applied nearest-neighbour matching and parametric imputation using the constraint $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. Parameters were estimated using maximum likelihood. The results are shown in Figure 3. In this example, there is a large difference between the nearest-neighbour imputations (b) and the parametric imputations (c). The true sample correlation between the $Y$ and $Z$ variables is $\rho_{YZ}=0.83$. Nearest-neighbour imputation gives an underestimate, with $\widehat{\rho}_{YZ}=0.14$. Parametric imputation gives much better estimate $\widehat{\rho}_{YZ}=0.85$. The $(Y,Z)$ distribution is recoverable in this scenario, as the generative model satisfies the identification restriction $\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$.

Data application

Iris dataset

We analysed a Anderson’s iris dataset (Anderson, 1935; Fisher, 1936) using a Gaussian mixture model. Figure 4 shows a pairs plot of the dataset with colour and shape giving the species labels.

We considered two matching scenarios, under the assumption that the species labels were unavailable. These are described in Table 2. Nearest-neighbour matching was compared to parametric imputation using the same identification constraints in Section 4.2. An important difference between the two scenarios is that in Scenario 1, the $\bm{X}$ single variable sepal with does not give enough information to reliably classify the observations into species. In Scenario 2, the sepal length and sepal width together give enough information to reliably classify the observations into species.

Figure 5 shows the results for Scenario 1. The asymptotic analysis showed that when the $\bm{X}$ variables are uninformative for clustering, we expect nearest-neighbour imputation to produce spurious clusters. In this example, we encounter mismatched species labels across datasets A and B when using sepal width to find nearest-neighbours. The parametric strategy does not create spurious clusters. The restrictions on the model based clustering approach encode a different set of assumptions that are more reasonable for a Gaussian mixture. As such, the parametric model imputed data shows more fidelity with the true data.

Figure 6 shows the results for Scenario 2. In this case nearest-neighbour imputation does not produce any spurious clusters. Spurious clusters are not expected here as the addition of sepal length into the common $\bm{X}$ group gives enough information to cluster the observations. The additional information should greatly reduce the number of mismatched species labels when finding nearest-neighbours across datasets. Parametric imputation gives very similar results to the nearest-neighbour method in Scenario 2. In both scenarios the imputed data shows weaker correlations than in the true dataset, this is because the component-wise identification restriction $\Sigma_{YZ}^{(h)}={\Sigma}_{YX}^{(h)}[{\Sigma}_{XX}^{(h)}]^{-1}{\Sigma}_{XZ}^{(h)}$ for $h=1,\ldots,g$ is perhaps not appropriate for this dataset.

Flow cytometry data (single cell subpopulation)

We analysed a subset of flow cytometry data from a study on graft vs host disease (Brinkman et al., 2007). See the supplementary material for more information on the dataset. The dataset is displayed in Figure 7. The data subset has $n=1162$ observations on $p=3$ markers. The dataset was split into two datasets of $n_{A}=n_{B}=681$ observations. We took FL1.H, FL4.H and FL3.H as the $\bm{X}$, $\bm{Y}$ and $\bm{Z}$ variables respectively. We applied nearest-neighbour matching and parametric imputation with a skew-normal model using the identification restriction described in Section 3.2. Figure 8 compares the imputed data to the actual data. Nearest-neighbour imputation gives different results than the parametric approach. The nearest-neighbour imputed data seems to exhibit the broad shape of the asymptotic imputation distribution that was studied in Section 4.1. If the skew-normal model is appropriate we expect to see weaker joint associations between the $\bm{Y}$ and $\bm{Z}$ variables in the imputed dataset than in the true dataset. It appears that the parametric model does a better job of recovering the general shape of the $(\bm{Y},\bm{Z})$ joint distribution. The true sample correlation between the $(\bm{Y},\bm{Z})$ variables is $\rho_{YZ}=0.39$. Parametric imputation gives a better estimate of the correlation ($\widehat{\rho}_{YZ}=0.44$) compared to nearest-neighbour imputation ($\widehat{\rho}_{YZ}=0.03$).

Flow cytometry data (two cell subpopulations)

We also consider another dataset from the graft vs host study. We extracted data on two cell subpopulations that were not considered in the previous data example. The data subset has $n=1321$ observations on $p=3$ markers. The dataset was split into two datasets of $n_{A}=660$, $n_{B}=661$ observations, and is plotted in Figure 9. We considered two different matching scenarios described in Table 3.

We applied nearest-neighbour matching and parametric imputation with a skew-normal mixture model in each scenario. Figure 10 compares the imputed data to the actual and data in Scenario one. The matching variable FL3.H is not sufficiently informative for clustering in Scenario 1, and nearest-neighbour imputation introduces spurious clusters. Focusing on panel (b), we see that the majority of points from the second cell subpopulation (red diamonds) are no longer located in the top right corner of the plot. Nearest-neighbour matching shifts most observations in cell subpopulation 2 to either the top-left or bottom-right of the plot. This is a serious distortion of the joint relationship between the FL4.H and FL2.H variables that is present in the original dataset. Looking at panel (c) we see that parameteric imputation does not suffer from this problem. We recover the correct general location of each cell cluster in the joint space of the FL4.H and FL2.H variables. Although the location of the second cell subpopulation is recovered well, we do not preserve the correct orientation of the second cell subpopulation.

Table 4 reports sample correlations for each cell subpopulation in Scenario one. Nearest-neighbour imputation gives poor estimates of the subpopulation correlations due to the mismatching of observations across clusters. In this scenario the parametric imputation strategy underestimates the $\bm{Y},\bm{Z}$ correlation for both groups.

Figure 11 compares the imputed data to the actual and data in Scenario two. Groupwise correlation estimates are also reported in Table 4. In Scenario 2, the matching variable FL3.H is able to separate the two cell subpopulations. As such, nearest-neighbour matching produces the correct number of clusters, however the correlation estimate in component two ($\widehat{\rho}_{YZ}=0.13$) is much lower than the true correlation in the source dataset ($\rho_{YZ}=0.74$). Parametric imputation also produces two clusters, but seems to preserve more of the distributional shape than the nearest-neighbour approach. Nearest-neighbour matching gives very small correlation estimates for each component. Parameteric imputation gives a good estimate of the correlation in subpopulation 2, but underestimates the correlation in subpopulation 1. It is interesting to compare the results in the two scenarios. Both parametric imputation and nearest-neighbour matching produces different correlation estimates in each scenario. The quality of results appears to be sensitive to the choice of matching variable. The matching variable in Scenario 1 is not informative for clustering and this appears to impact the faithfulness of the imputed data for both methods.

Conclusion

In the statistical matching problem, the identifiability of $f(\bm{x},\bm{y},\bm{z};\bm{\theta})$ does not necessarily imply the conditional independence of $\bm{Y}$ and $\bm{Z}$ given $\bm{X}$. As such, a parametric imputation strategy can give different results to nearest-neighbour imputation. Although nearest-neighbour imputation is a non-parameteric method, it is not assumption free. The nearest-neighbour method makes a conditional independence assumption on the observed random variables that may not be appropriate if the generative model has some latent structure.

We showed that a model based approach to the statistical matching problem for skew-normal data and Gaussian mixture data can be implemented using data augmentation and the EM algorithm. By reparameterising the complete-data likelihood we can establish a relationship between the statistical matching problem and mixtures of regression models. The reparameterisation is important to give closed form E and M steps. Additionally, we can impose identification constraints on the model to avoid introducing the missing $\bm{Y}$ and $\bm{Z}$ observations into the complete-data log likelihood. This is to minimise the computational expense of each EM iteration. The statistical matching of non-Gaussian data can be significantly more computationally demanding than non-Gaussian data, and this is an important consideration when working with large flow cytometry datasets.

The artefacts introduced by the nearest-neighbour method can lead to potential issues in downstream analysis, as the imputer’s model and the analyst’s model should be compatible (Rubin, 1996; Meng, 1994). For example, if nearest-neighbour matching introduces spurious clusters, a downstream analyst using mixture models may overestimate the number of subpopulations in the dataset. The conditional independence assumption is not a necessity in order to perform statistical matching, and as such it may be worth exploring alternative non-parametric or semi-parametric imputation methods.

The statistical matching of mixture models is particularly interesting as the group separation on the common $\bm{X}$ variables appears to strongly influence parameter estimation. In some situations, investigators may be able to engineer the data collection process such that informative variables are assigned to the common $\bm{X}$ group. This situation may arise in survey research where respondents are asked different sets of questions, with a smaller set of common questions given to every subject. The selection of common questions $\bm{X}$ is interesting from an experimental design point of view.

The statistical matching problem is a challenging data integration scenario that can require specialised algorithms for missing data imputation. We have found that model based approaches can give different results to nonparameteric imputation schemes. However, the pathological nature of the statistical matching problem does limit the accuracy of any imputation scheme. The model based approach can introduce some systemic bias into the imputed data when there is model misspecification. The identification restrictions adopted in this work are all inspired by the conditional independence constraint for the Gaussian distribution $\Sigma_{YZ}=\Sigma_{YZ}={\Sigma}_{YX}{\Sigma}_{XX}^{-1}{\Sigma}_{XZ}$. This assumption is quite strong, and is not testable. We have seen in the data applications that there is no guarantee that it will hold on real data. As such, an important future research direction is the identification of alternative parameter constraints that are testable.

Supporting Information

Additional supporting information may be found online in the Supporting Information section at the end of the article.

Appendix A Appendix

Conditional skew-normal

The conditional distributions $f(\bm{y}_{i}|\bm{x}_{i};\bm{\theta})$ and $f(\bm{z}_{i}|\bm{x}_{i};\bm{\theta})$ belong to a generalised family of skew-normal distributions that also have a latent variable representation (Areallano-Valle and Azzalini, 2006; Arellano-Valle et al., 2006). Let $TN(\mu,\sigma^{2},a)$ denote a lower truncated normal distribution where $\mu$ and $\sigma^{2}$ give the mean and variance of the underling normal distribution and $a$ gives the lower truncation bound. We say $U\sim TN(\mu,\sigma^{2},a)$ if $U\overset{d}{=}[W|W>a]$ where $W\sim N(\mu,\sigma^{2})$. The conditional distribution of $\bm{Y}$ and $\bm{Z}$ given $\bm{X}$ can be represented as

[TABLE]

where ${U}_{X}$ has a truncated normal distribution $TN(\tau_{X},\gamma_{X},0)$, and $\bm{V}_{Y|X}$ and $\bm{V}_{Z|X}$ are jointly normally distributed. The location parameters are given by

[TABLE]

The latent $U_{X}$ is distributed as a truncated normal $TN(\tau_{X},\gamma_{X},0)$ random variable, where

[TABLE]

The latent $\bm{V}_{Y|X}$ and $\bm{V}_{Z|X}$ have a normal distribution,

[TABLE]

The conditional covariance matrix is

[TABLE]

The conditional skewness parameters are

[TABLE]

Asymptotic analysis of skew-normal

The hierarchical representation of a skew-normal random variable is useful to analyse the asymptotic behaviour of the nearest-neighbour method for skew-normal data. The asymptotic form of the nearest-neighbour imputation distribution (37) will involve the true conditional distributions $f(\bm{y}_{i}|\bm{x}_{i};\bm{\theta})$ and $f(\bm{z}_{i}|\bm{x}_{i};\bm{\theta})$. The true marginal $(\bm{Y}_{i},\bm{Z}_{i})$ distribution can be expressed through a hierarchical model

[TABLE]

We can also characterise the asymptotic nearest-neighbour imputation distribution through a hierarchical model. Asymptotically, the effective model used by nearest-neighbour matching is that conditional on $\bm{X}_{i}=\bm{x}_{i}$,

[TABLE]

Comparing to (A.6) to (A.2) we see that an important difference between the nearest-neighbour imputation distribution and the true marginal distribution is the use of two independent latent scaling variables ${U}_{X}^{(1)}$ and ${U}_{X}^{(2)}$ instead of the single latent variable ${U}_{X}$ that appears in the true model. A second difference is that $\bm{Z}_{i}$ and $\bm{Y}_{i}$ are sampled independently conditional on $\bm{X}_{i}$. Levels (A.4) and (A.5) are independent conditional on the lower layers. The conditional covariance $\Sigma_{YZ|X}$ that appears in the true generative process (A.2) does not influence the nearest-neighbour imputation distribution. With these systematic differences, it is reasonable to expect that nearest-neighbour imputation will underestimate the level of dependence between the $\bm{Y}$ and $\bm{Z}$ variables. If we believe that the skew-normal model (19) is appropriate, it is difficult to justify the use of nearest-neighbour matching as the asymptotic behaviour is not consistent with our initial assumptions.

Asymptotic analysis of Gaussian mixture

The true joint distribution on the $(\bm{Y}_{i},\bm{Z}_{i})$ variables has the hierarchical representation:

[TABLE]

where for $j=1,\ldots,g$:

[TABLE]

Although conventional to include the cluster indicator $S_{i}$ at the bottom of the hierarchy, for the purposes of our analysis it is beneficial to sample it conditional on $\bm{X}_{i}=\bm{x}_{i}$. We again analyse the asymptotic behaviour nearest-neighbour matching by determining the form of the asymptotic imputation distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})$ (37). The true conditional distributions $f(\bm{y}_{i}|\bm{x}_{i};\bm{\theta})$ and $f(\bm{z}_{i}|\bm{x}_{i};\bm{\theta})$ can also be expressed as finite mixture models. In full,

[TABLE]

where the weights ${\tau}_{1},\ldots,\tau_{g}$ are given by (A.11). Asymptotically speaking, the nearest-neighbour imputation distribution $g(\bm{x}_{i},\bm{y}_{i},\bm{z}_{i})$ has the hierarchical representation:

[TABLE]

An important difference between the nearest-neighbour imputation distribution and the true conditional model is the introduction of a second latent variable $T_{i}$ (A.14). If the sampled $T_{i}$ is not equal to the sampled $S_{i}$ at level (A.14), the $\bm{Y}_{i}$ and $\bm{Z}_{i}$ random variables will not be drawn from the same component distribution in levels (A.13) and (A.12). The probability of a mismatch is related to how reliably the $\bm{X}$ variables can be used to classify observations, entering through the $\tau_{j}$ terms (A.11). The imputed joint distribution $g(\bm{y},\bm{z})$ can introduce spurious clusters if the $\bm{X}$ variables are not useful for clustering the observations.

Web Appendix:Introduction

Here we give the E-step and M-step for the EM algorithms for the statistical matching of non-Gaussian data. Recall the hierarchical model for the skew normal distribution:

[TABLE]

where ${U}_{i}\sim TN(0,1,0)$ and

[TABLE]

The combined data matrix from dataset $A$ and dataset $B$ is partitioned as

[TABLE]

where shaded cells are observed and blank cells are missing.

Web Appendix A: Skew-normal distribution

E-step

Using the hierarchical model for the skew normal distribution (19), define the following conditional expectations given the current parameter estimates:

[TABLE]

Let

[TABLE]

Furthermore define the constants

[TABLE]

Let $\phi(\cdot)$ denote the standard normal density function and let $\Phi(\cdot)$ denote the cumulative distribution function of the standard normal distribution. We have the closed form expressions for the expectations (Pyne et al., 2009).

[TABLE]

M-step

For $i=1,\ldots,n$ let

[TABLE]

The maximisation of the $X$ parameters follows from existing results on the skew normal distribution (Pyne et al., 2009):

[TABLE]

Recall that $\Gamma_{A}$ and $\Gamma_{B}$ contain the regression parameters for the conditional regression models: $\Gamma_{A}=[\bm{\alpha}_{Y}\ \bm{\lambda}_{Y}\ \bm{\beta}_{Y}]^{\mathsf{T}}$ and $\Gamma_{B}=[\bm{\alpha}_{Z}\ \bm{\lambda}_{Z}\ \bm{\beta}_{Z}]^{\mathsf{T}}$. The regressions for dataset A and dataset B can be written as

[TABLE]

where $\bm{\epsilon}_{A}\sim MN(\bm{I}_{n_{A}},\Omega_{Y})$ and $\bm{\epsilon}_{B}\sim MN(\bm{I}_{n_{B}},\Omega_{Z})$. The design matrices for the regressions now include the latent $U_{i}$ terms. The complete-data design matrices are given by

[TABLE]

We need to calculate the conditional expectations $\mathbb{E}[\bm{B}^{\mathsf{T}}_{A}\bm{B}_{A}]$ and $\mathbb{E}[\bm{B}^{\mathsf{T}}_{B}\bm{B}_{B}]$ given the current parameter estimates and the observed data. For $i=1,\ldots,n$ let

[TABLE]

From standard results on multiple outcome regression models(Rencher and Christensen, 2012), we have a closed form M-step for the regression parameters

[TABLE]

For $i=1,\ldots,n$ let

[TABLE]

The M-step for the error covariance matrices in the regressions (A.17) and (A.18) is as follows:

[TABLE]

The estimates in terms of the original parameters of the skew-normal density can be obtained by substituting into the following equations:

[TABLE]

Web Appendix B: Mixtures of Gaussians

E-step

For $h=1,\ldots,g$, define the posterior class probabilities as

[TABLE]

M-step

For the maximisation of the $X$ parameters we follow standard results for Gaussian mixture models.

[TABLE]

Recall the conditional regression specification for dataset A,

[TABLE]

where $\bm{\beta}_{Y}^{(h)}=\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}$, $\bm{\alpha}_{Y}^{(h)}=\bm{\mu}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\mu}_{X}^{(h)}$ and $\Omega_{Y}^{(h)}=\Sigma_{YY}^{(h)}-\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XY}^{(h)}$. A similar model applies for dataset B,

[TABLE]

where $\bm{\beta}_{Z}^{(h)}=\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}$, $\bm{\alpha}_{Z}^{(h)}=\bm{\mu}_{Z}^{(h)}-\bm{\beta}_{Z}\bm{\mu}_{X}$ and $\Omega_{Z}=\Sigma_{ZZ}-\Sigma_{ZX}\Sigma_{XX}^{-1}\Sigma_{XZ}$. The maximisation for the conditional regression parameters then follows from existing results on mixtures of regression models (Jones and McLachlan, 1992). For $i=1,\ldots,n$ define:

[TABLE]

The $M$-step for the regression coefficients is given by:

[TABLE]

For the residual error variances in (A.20) and (A.21) we have

[TABLE]

For $h-1,\ldots,g$, the maximum likelihood estimates of the original parameters can then be obtained by substituting into the following formulae:

[TABLE]

Web Appendix C: Mixtures of Skew-normal distributions

E-step

For $h=1,\ldots,g$, define the posterior class probabilities as

$\displaystyle\tau_{i,h}$ $\displaystyle=\begin{cases}f\left(\begin{bmatrix}\bm{x}_{i}\\ \bm{y}_{i}\end{bmatrix};\begin{bmatrix}\bm{\mu}_{X}^{(h)}\\ \bm{\mu}_{Y}^{(h)}\end{bmatrix},\begin{bmatrix}\Sigma_{XX}^{(h)}&\Sigma_{XY}^{(h)}\\ \Sigma_{YX}^{(h)}&\Sigma_{YY}^{(h)}\end{bmatrix},\begin{bmatrix}\bm{\delta}_{X}^{(h)}\\ \bm{\delta}_{Y}^{(h)}\end{bmatrix}\right)\left\{\sum_{r=1}^{g}f\left(\begin{bmatrix}\bm{x}_{i}\\ \bm{y}_{i}\end{bmatrix};\begin{bmatrix}\bm{\mu}_{X}^{(r)}\\ \bm{\mu}_{Y}^{(r)}\end{bmatrix},\begin{bmatrix}\Sigma_{XX}^{(r)}&\Sigma_{XY}^{(r)}\\ \Sigma_{YX}^{(r)}&\Sigma_{YY}^{(r)}\end{bmatrix},\begin{bmatrix}\bm{\delta}_{X}^{(r)}\\ \bm{\delta}_{Y}^{(r)}\end{bmatrix}\right)\right\}^{-1}&\mbox{if }i=1,\ldots,n_{A}\\ \\ f\left(\begin{bmatrix}\bm{x}_{i}\\ \bm{z}_{i}\end{bmatrix};\begin{bmatrix}\bm{\mu}_{X}^{(h)}\\ \bm{\mu}_{Z}^{(h)}\end{bmatrix},\begin{bmatrix}\Sigma_{XX}^{(h)}&\Sigma_{XZ}^{(h)}\\ \Sigma_{ZX}^{(h)}&\Sigma_{ZZ}^{(h)}\end{bmatrix},\begin{bmatrix}\bm{\delta}_{X}^{(h)}\\ \bm{\delta}_{Z}^{(h)}\end{bmatrix}\right)\left\{\sum_{r=1}^{g}f\left(\begin{bmatrix}\bm{x}_{i}\\ \bm{z}_{i}\end{bmatrix};\begin{bmatrix}\bm{\mu}_{X}^{(r)}\\ \bm{\mu}_{Z}^{(r)}\end{bmatrix},\begin{bmatrix}\Sigma_{XX}^{(r)}&\Sigma_{XZ}^{(r)}\\ \Sigma_{ZX}^{(r)}&\Sigma_{ZZ}^{(r)}\end{bmatrix},\begin{bmatrix}\bm{\delta}_{X}^{(r)}\\ \bm{\delta}_{Z}^{(r)}\end{bmatrix}\right)\right\}^{-1}&\mbox{if }i=n_{A}+1,\ldots,n\end{cases}$

Where $f(\bm{v};\bm{\mu},\Sigma,\bm{\delta})$ denotes the skew normal density with parameters $\bm{\mu}$, $\Sigma$ and $\bm{\delta}$ evaluated at the point $\bm{v}$. Using the hierarchical model for the skew normal distribution (A.16), define the following conditional expectations given

[TABLE]

For $h=1,\ldots,g$, define

[TABLE]

Furthermore define the constants

[TABLE]

Let $\phi(\cdot)$ denote the standard normal density function and let $\Phi(\cdot)$ denote the cumulative distribution function of the standard normal distribution. We have the closed form expressions for the expectations.

[TABLE]

M-step

For $i=1,\ldots,n$ let

[TABLE]

The maximisation of the $X$ parameters follows from existing results on mixtures of skew normal distributions (Pyne et al., 2009). For $h=1,\ldots,g$:

[TABLE]

The complete-data model can be expressed in terms of conditional regression specifications. Let $\bm{\beta}_{Y}^{(h)}=\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1},\bm{\alpha}_{Y}^{(h)}=\bm{\mu}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\mu}_{X}^{(h)},\Omega_{Y}^{(h)}=\Sigma_{YY}^{(h)}-\Sigma_{YX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XY}^{(h)}$ and $\bm{\lambda}_{Y}^{(h)}=\bm{\delta}_{Y}^{(h)}-\bm{\beta}_{Y}^{(h)}\bm{\delta}_{X}^{(h)}$.

[TABLE]

Let $\bm{\beta}_{Z}^{(h)}=\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1},\bm{\alpha}_{Z}^{(h)}=\bm{\mu}_{Z}^{(h)}-\bm{\beta}_{Z}^{(h)}\bm{\mu}_{X}^{(h)},\Omega_{Z}^{(h)}=\Sigma_{ZZ}^{(h)}-\Sigma_{ZX}^{(h)}[\Sigma_{XX}^{(h)}]^{-1}\Sigma_{XZ}^{(h)}$ and $\bm{\lambda}_{Z}^{(h)}=\bm{\delta}_{Z}^{(h)}-\bm{\beta}_{Z}^{(h)}\bm{\delta}_{X}^{(h)}$. Similarly the conditional distribution of $\bm{Z}_{i}$ given $\bm{X}_{i}$ and the latent scaling variable $U_{i}$ can be represented as a regression model

[TABLE]

We need to calculate the conditional expectations of the predictors in the regressions given the current parameter estimates and the observed data. For $i=1,\ldots,nn$ and $h=1,\ldots,g$ let

[TABLE]

Recall that $\Gamma_{A}^{(h)}$ and $\Gamma_{B}^{(h)}$ contain the regression parameters for component $h$, so $\Gamma_{A}^{(h)}=[\bm{\alpha}_{Y}^{(h)}\ \bm{\lambda}_{Y}^{(h)}\ \bm{\beta}_{Y}^{(h)}]^{\mathsf{T}}$ and $\Gamma_{B}^{(h)}=[\bm{\alpha}_{Z}^{(h)}\ \bm{\lambda}_{Z}^{(h)}\ \bm{\beta}_{Z}]^{\mathsf{T}}$ for $h=1,\ldots,g$. We have a closed form M-step for the $\bm{Y}$ regression parameters in (A.22) and the $\bm{Z}$ regression parameters parameters in (A.23):

[TABLE]

For $h=1,\ldots,g$ define:

[TABLE]

The M-step for the error covariance matrices in the regressions (34) and (35) is as follows:

[TABLE]

The estimates in the original parameterisation for the skew-normal model can be obtained using the following results. For $h=1,\ldots,g$ we have the relationships:

[TABLE]