Probabilistic partial least squares model: Identifiability, estimation   and application

Said el Bouhaddani; Hae-Won Uh; Caroline Hayward; Geurt Jongbloed,; Jeanine Houwing-Duistermaat

arXiv:1706.03597·stat.ME·March 26, 2021·J. Multivar. Anal.

Probabilistic partial least squares model: Identifiability, estimation and application

Said el Bouhaddani, Hae-Won Uh, Caroline Hayward, Geurt Jongbloed,, Jeanine Houwing-Duistermaat

PDF

TL;DR

This paper introduces Probabilistic PLS (PPLS), a new model that ensures parameter identifiability, provides a statistical inference framework, and performs well in high-dimensional, non-normal data scenarios, demonstrated through simulations and real data.

Contribution

The paper develops a probabilistic formulation of PLS with identifiable parameters and an EM algorithm for estimation, addressing limitations of traditional PLS methods.

Findings

01

PPLS parameters are identifiable up to sign.

02

Estimates are robust against non-normality.

03

Method performs well in high-dimensional settings.

Abstract

With a rapid increase in volume and complexity of data sets, there is a need for methods that can extract useful information, for example the relationship between two data sets measured for the same persons. The Partial Least Squares (PLS) method can be used for this dimension reduction task. Within life sciences, results across studies are compared and combined. Therefore, parameters need to be identifiable, which is not the case for PLS. In addition, PLS is an algorithm, while epidemiological study designs are often outcome-dependent and methods to analyze such data require a probabilistic formulation. Moreover, a probabilistic model provides a statistical framework for inference. To address these issues, we develop Probabilistic PLS (PPLS). We derive maximum likelihood estimators that satisfy the identifiability conditions by using an EM algorithm with a constrained optimization in…

Tables2

Table 1. Table 1 : Overview of the simulation scenarios. The noise level is defined as the proportion of variation in the noise matrices E 𝐸 E , F 𝐹 F and H 𝐻 H relative to the total variation in X 𝑋 X , Y 𝑌 Y and U 𝑈 U respectively.

Sample size	$N = (50, 500)$
Dimensionality	$p = q = (20, 1000)$
Noise level	$α_{n} = (0.1, 0.5)$
Distribution of $t$ , $u$ , $e$ , $f$ and $h$	${𝒩 (0, 1), t_{2}, 𝒫 (1), ℬ ℐ 𝒩 (2, 0.3)}$

Table 2. Table 2 : Proportion of correct order of loadings W 𝑊 W and C 𝐶 C across 1000 simulation replicates. These were obtained for different values of the dimensionality (high = 1000 variables, low = 20 variables), sample size (large = 500 subjects, small = 50 subjects) and noise level (high = 50% noise, low = 10% noise).

Dimensionality	Sample Size	Noise Level	Correct Ordering Proportion
low	large	low	1.000
	large	high	0.989
	small	low	0.932
	small	high	0.435
high	large	low	0.990
	large	high	0.985
	small	low	0.940
	small	high	0.665

Equations92

x = t W^{⊤} + e, y = u C^{⊤} + f, u = tB + h .

x = t W^{⊤} + e, y = u C^{⊤} + f, u = tB + h .

e \sim N (0, σ_{e}^{2} I_{p}), f \sim N (0, σ_{f}^{2} I_{q}), h \sim N (0, σ_{h}^{2} I_{r}) .

e \sim N (0, σ_{e}^{2} I_{p}), f \sim N (0, σ_{f}^{2} I_{q}), h \sim N (0, σ_{h}^{2} I_{r}) .

t \sim N (0, Σ_{t}) .

t \sim N (0, Σ_{t}) .

θ = (W, C, B, Σ_{t}, σ_{e}, σ_{f}, σ_{h}) .

θ = (W, C, B, Σ_{t}, σ_{e}, σ_{f}, σ_{h}) .

f (x, y) = (2 π)^{- (p + q) /2} ∣Σ ∣^{- \frac{1}{2}} e^{(x, y) Σ^{- 1} (x, y)^{⊤}},

f (x, y) = (2 π)^{- (p + q) /2} ∣Σ ∣^{- \frac{1}{2}} e^{(x, y) Σ^{- 1} (x, y)^{⊤}},

Σ = (Σ_{x} Σ_{y, x} Σ_{x, y} Σ_{y}) = (W Σ_{t} W^{⊤} + σ_{e}^{2} I_{p} C B Σ_{t} W^{⊤} W Σ_{t} B C^{⊤} C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q}) .

Σ = (Σ_{x} Σ_{y, x} Σ_{x, y} Σ_{y}) = (W Σ_{t} W^{⊤} + σ_{e}^{2} I_{p} C B Σ_{t} W^{⊤} W Σ_{t} B C^{⊤} C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q}) .

x = t R_{S} R_{S}^{⊤} W^{⊤} + e = j = 1 \sum r t_{j} (R_{S})_{j j}^{2} w_{j}^{⊤} + e .

x = t R_{S} R_{S}^{⊤} W^{⊤} + e = j = 1 \sum r t_{j} (R_{S})_{j j}^{2} w_{j}^{⊤} + e .

x = t W^{⊤} + e = tR R^{⊤} W^{⊤} + e,

x = t W^{⊤} + e = tR R^{⊤} W^{⊤} + e,

L (θ) = - \frac{N ( p + q )}{2} - \frac{N}{2} ln ∣Σ∣ - \frac{N}{2} tr (S Σ^{- 1})

L (θ) = - \frac{N ( p + q )}{2} - \frac{N}{2} ln ∣Σ∣ - \frac{N}{2} tr (S Σ^{- 1})

f (x, y, t, u) = f (x ∣ t) f (y ∣ u) f (u ∣ t) f (t) .

f (x, y, t, u) = f (x ∣ t) f (y ∣ u) f (u ∣ t) f (t) .

f (x, y, t, u) = f (x, y ∣ t, u) f (t, u) = f (x ∣ t, u) f (y ∣ t, u) f (t, u) .

f (x, y, t, u) = f (x, y ∣ t, u) f (t, u) = f (x ∣ t, u) f (y ∣ t, u) f (t, u) .

ln f (X ∣ T) ln f (Y ∣ U) ln f (U ∣ T) = - \frac{N p}{2 π σ _{e}^{2}} - \frac{1}{2 σ _{e}^{2}} i = 1 \sum N ∣∣ X_{i} - T_{i} W^{⊤} ∣ ∣^{2}, = - \frac{N q}{2 π σ _{f}^{2}} - \frac{1}{2 σ _{f}^{2}} i = 1 \sum N ∣∣ Y_{i} - U_{i} C^{⊤} ∣ ∣^{2}, = - \frac{N r}{2 π σ _{h}^{2}} - \frac{1}{2 σ _{h}^{2}} i = 1 \sum N ∣∣ U_{i} - T_{i} B ∣ ∣^{2} .

ln f (X ∣ T) ln f (Y ∣ U) ln f (U ∣ T) = - \frac{N p}{2 π σ _{e}^{2}} - \frac{1}{2 σ _{e}^{2}} i = 1 \sum N ∣∣ X_{i} - T_{i} W^{⊤} ∣ ∣^{2}, = - \frac{N q}{2 π σ _{f}^{2}} - \frac{1}{2 σ _{f}^{2}} i = 1 \sum N ∣∣ Y_{i} - U_{i} C^{⊤} ∣ ∣^{2}, = - \frac{N r}{2 π σ _{h}^{2}} - \frac{1}{2 σ _{h}^{2}} i = 1 \sum N ∣∣ U_{i} - T_{i} B ∣ ∣^{2} .

Q (θ) = E {L_{Comp} (θ) ∣ X, Y, θ^{'}},

Q (θ) = E {L_{Comp} (θ) ∣ X, Y, θ^{'}},

W^{k + 1} = X^{⊤} E (T ∣ X, Y, θ^{k}) (L_{W}^{⊤})^{- 1}; C^{k + 1} = Y^{⊤} E (U ∣ X, Y, θ^{k}) (L_{C}^{⊤})^{- 1}; B^{k + 1} = E (U^{⊤} T ∣ X, Y, θ^{k}) {E (T^{⊤} T ∣ X, Y, θ^{k})}^{- 1} \circ I_{r}; Σ_{t}^{k + 1} = \frac{1}{N} E (T^{⊤} T ∣ X, Y, θ^{k}) \circ I_{r}; (σ_{h}^{2})^{k + 1} = \frac{1}{N r} tr E (H^{⊤} H ∣ X, Y, θ^{k}); (σ_{e}^{2})^{k + 1} = \frac{1}{N p} tr E (E^{⊤} E ∣ X, Y, θ^{k}); (σ_{f}^{2})^{k + 1} = \frac{1}{N q} tr E (F^{⊤} F ∣ X, Y, θ^{k});

W^{k + 1} = X^{⊤} E (T ∣ X, Y, θ^{k}) (L_{W}^{⊤})^{- 1}; C^{k + 1} = Y^{⊤} E (U ∣ X, Y, θ^{k}) (L_{C}^{⊤})^{- 1}; B^{k + 1} = E (U^{⊤} T ∣ X, Y, θ^{k}) {E (T^{⊤} T ∣ X, Y, θ^{k})}^{- 1} \circ I_{r}; Σ_{t}^{k + 1} = \frac{1}{N} E (T^{⊤} T ∣ X, Y, θ^{k}) \circ I_{r}; (σ_{h}^{2})^{k + 1} = \frac{1}{N r} tr E (H^{⊤} H ∣ X, Y, θ^{k}); (σ_{e}^{2})^{k + 1} = \frac{1}{N p} tr E (E^{⊤} E ∣ X, Y, θ^{k}); (σ_{f}^{2})^{k + 1} = \frac{1}{N q} tr E (F^{⊤} F ∣ X, Y, θ^{k});

L_{W} L_{W}^{⊤} = E (T^{⊤} ∣ X, Y, θ^{k}) X X^{⊤} E (T ∣ X, Y, θ^{k}), L_{C} L_{C}^{⊤} = E (U^{⊤} ∣ X, Y, θ^{k}) Y Y^{⊤} E (U ∣ X, Y, θ^{k}) .

L_{W} L_{W}^{⊤} = E (T^{⊤} ∣ X, Y, θ^{k}) X X^{⊤} E (T ∣ X, Y, θ^{k}), L_{C} L_{C}^{⊤} = E (U^{⊤} ∣ X, Y, θ^{k}) Y Y^{⊤} E (U ∣ X, Y, θ^{k}) .

(W^{k + 1})^{⊤} W^{k + 1} = L_{W}^{- 1} \tilde{T}^{⊤} X X^{⊤} \tilde{T} (L_{W}^{⊤})^{- 1} = L_{W}^{- 1} L_{W} L_{W}^{⊤} (L_{W}^{⊤})^{- 1} = I_{r},

(W^{k + 1})^{⊤} W^{k + 1} = L_{W}^{- 1} \tilde{T}^{⊤} X X^{⊤} \tilde{T} (L_{W}^{⊤})^{- 1} = L_{W}^{- 1} L_{W} L_{W}^{⊤} (L_{W}^{⊤})^{- 1} = I_{r},

E {B (\hat{θ}) ∣ X, Y} - E {S (\hat{θ}) S (\hat{θ})^{⊤} ∣ X, Y} .

E {B (\hat{θ}) ∣ X, Y} - E {S (\hat{θ}) S (\hat{θ})^{⊤} ∣ X, Y} .

w_{j, k} = N {j; (1/2 + 1/10 j) p, 1/10 p}, c_{j, k} = N {j; (3/5 + 1/10 j) q, 1/10 q} .

w_{j, k} = N {j; (1/2 + 1/10 j) p, 1/10 p}, c_{j, k} = N {j; (3/5 + 1/10 j) q, 1/10 q} .

var (x) var (y) cov (x, y) = var (t W^{⊤} + e) = W var (t) W^{⊤} + var (e) = W Σ_{t} W^{⊤} + σ_{e}^{2} I_{p}, = var (U C^{⊤} + f) = C var (u) C^{⊤} + var (f) = C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q}, = cov (t W^{⊤} + e, u C^{⊤} + f) = W cov (t, u) C^{⊤} = W cov (t, tB) C^{⊤} = W B Σ_{t} C^{⊤} .

var (x) var (y) cov (x, y) = var (t W^{⊤} + e) = W var (t) W^{⊤} + var (e) = W Σ_{t} W^{⊤} + σ_{e}^{2} I_{p}, = var (U C^{⊤} + f) = C var (u) C^{⊤} + var (f) = C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q}, = cov (t W^{⊤} + e, u C^{⊤} + f) = W cov (t, u) C^{⊤} = W cov (t, tB) C^{⊤} = W B Σ_{t} C^{⊤} .

cov (x, t) cov (x, u) cov (y, t) cov (y, u) = cov (t W^{⊤} + e, t) = W var (t) = W Σ_{t}, = cov (t W^{⊤} + e, tB + h) = W var (t) B = W Σ_{t} B, = cov (u C^{⊤} + f, t) = C cov (tB + h, t) = C Σ_{t} B, = cov (u C^{⊤} + f, u) = C cov (tB + h, tB + h) = C (Σ_{t} B^{2} + σ_{h}^{2} I_{r}) .

cov (x, t) cov (x, u) cov (y, t) cov (y, u) = cov (t W^{⊤} + e, t) = W var (t) = W Σ_{t}, = cov (t W^{⊤} + e, tB + h) = W var (t) B = W Σ_{t} B, = cov (u C^{⊤} + f, t) = C cov (tB + h, t) = C Σ_{t} B, = cov (u C^{⊤} + f, u) = C cov (tB + h, tB + h) = C (Σ_{t} B^{2} + σ_{h}^{2} I_{r}) .

Σ = \tilde{Σ} \Leftrightarrow θ = \tilde{θ},

Σ = \tilde{Σ} \Leftrightarrow θ = \tilde{θ},

A_{1} A_{1}^{⊤} = W D^{2} W^{⊤} = \tilde{W} \tilde{D}^{2} \tilde{W}^{⊤} = A_{2} A_{2}^{⊤}; A_{1}^{⊤} A_{1} = C D^{2} C^{⊤} = \tilde{C} \tilde{D}^{2} \tilde{C}^{⊤} = A_{2}^{⊤} A_{2} .

A_{1} A_{1}^{⊤} = W D^{2} W^{⊤} = \tilde{W} \tilde{D}^{2} \tilde{W}^{⊤} = A_{2} A_{2}^{⊤}; A_{1}^{⊤} A_{1} = C D^{2} C^{⊤} = \tilde{C} \tilde{D}^{2} \tilde{C}^{⊤} = A_{2}^{⊤} A_{2} .

C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q} = C (B^{2} Σ_{t} + \tilde{σ}_{h}^{2} I_{r}) C^{⊤} + \tilde{σ}_{f}^{2} I_{q},

C (B^{2} Σ_{t} + σ_{h}^{2} I_{r}) C^{⊤} + σ_{f}^{2} I_{q} = C (B^{2} Σ_{t} + \tilde{σ}_{h}^{2} I_{r}) C^{⊤} + \tilde{σ}_{f}^{2} I_{q},

Σ_{x, y} = \tilde{Σ}_{x, y}, Σ_{x} = \tilde{Σ}_{x}, Σ_{y} = \tilde{Σ}_{y} .

Σ_{x, y} = \tilde{Σ}_{x, y}, Σ_{x} = \tilde{Σ}_{x}, Σ_{y} = \tilde{Σ}_{y} .

(Σ_{z} Σ_{x, z} Σ_{z, x} Σ_{x}) .

(Σ_{z} Σ_{x, z} Σ_{z, x} Σ_{x}) .

μ_{t} = E (t ∣ x, y, θ) = (x, y) Σ^{- 1} cov {(x, y), t}, μ_{u} = E (u ∣ x, y, θ) = (x, y) Σ^{- 1} cov {(x, y), u} .

μ_{t} = E (t ∣ x, y, θ) = (x, y) Σ^{- 1} cov {(x, y), t}, μ_{u} = E (u ∣ x, y, θ) = (x, y) Σ^{- 1} cov {(x, y), u} .

C_{T T} = E (t^{⊤} t ∣ x, y, θ) C_{U U} = E (u^{⊤} u ∣ x, y, θ) = I_{r} - cov {t, (x, y)} Σ^{- 1} cov {(x, y), t} + cov {t, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), t}, = I_{r} - cov {u, (x, y)} Σ^{- 1} cov {(x, y), u} + cov {u, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), u},

C_{T T} = E (t^{⊤} t ∣ x, y, θ) C_{U U} = E (u^{⊤} u ∣ x, y, θ) = I_{r} - cov {t, (x, y)} Σ^{- 1} cov {(x, y), t} + cov {t, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), t}, = I_{r} - cov {u, (x, y)} Σ^{- 1} cov {(x, y), u} + cov {u, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), u},

C_{U T} = E (u^{⊤} t ∣ x, y, θ) = Σ_{t} B - cov {u, (x, y)} Σ^{- 1} cov {(x, y), t} + cov {u, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), t}

C_{U T} = E (u^{⊤} t ∣ x, y, θ) = Σ_{t} B - cov {u, (x, y)} Σ^{- 1} cov {(x, y), t} + cov {u, (x, y)} Σ^{- 1} S Σ^{- 1} cov {(x, y), t}

cov {(x, y), t} = (W Σ_{t} C Σ_{t} B), cov {(x, y), u} = (W Σ_{t} B C (Σ_{t} B + σ_{h}^{2} I_{r})) .

cov {(x, y), t} = (W Σ_{t} C Σ_{t} B), cov {(x, y), u} = (W Σ_{t} B C (Σ_{t} B + σ_{h}^{2} I_{r})) .

W, σ_{e} f (x ∣ t) C, σ_{f} f (y ∣ u) B, σ_{h} f (u ∣ t) Σ_{t} f (t)

W, σ_{e} f (x ∣ t) C, σ_{f} f (y ∣ u) B, σ_{h} f (u ∣ t) Σ_{t} f (t)

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Probabilistic partial least squares model:

Identifiability, estimation and application

Said el Bouhaddani

S.el_ [email protected]

Hae-Won Uh

[email protected]

Caroline Hayward

[email protected]

Geurt Jongbloed

[email protected]

Jeanine Houwing-Duistermaat

[email protected]

Department of Medical statistics and bioinformatics, Leiden University Medical Center, The Netherlands

Department of Biostatistics and Research Support UMC Utrecht, div. Julius Centrum, University Medical Center Utrecht, The Netherlands

Department of Applied Mathematics, Delft University of Technology, The Netherlands

Department of Statistics, University of Leeds, United Kingdom

MRC Human Genetics Unit, Institute of Genetics and Molecular Medicine, University of Edinburgh, Scotland

Abstract

With a rapid increase in volume and complexity of data sets, there is a need for methods that can extract useful information, for example the relationship between two data sets measured for the same persons. The Partial Least Squares (PLS) method can be used for this dimension reduction task. Within life sciences, results across studies are compared and combined. Therefore, parameters need to be identifiable, which is not the case for PLS. In addition, PLS is an algorithm, while epidemiological study designs are often outcome-dependent and methods to analyze such data require a probabilistic formulation. Moreover, a probabilistic model provides a statistical framework for inference. To address these issues, we develop Probabilistic PLS (PPLS). We derive maximum likelihood estimators that satisfy the identifiability conditions by using an EM algorithm with a constrained optimization in the M step. We show that the PPLS parameters are identifiable up to sign. A simulation study is conducted to study the performance of PPLS compared to existing methods. The PPLS estimates performed well in various scenarios, even in high dimensions. Most notably, the estimates seem to be robust against departures from normality. To illustrate our method, we applied it to IgG glycan data from two cohorts. Our PPLS model provided insight as well as interpretable results across the two cohorts.

keywords:

Dimension reduction , EM algorithm , Identifiability , Inference , Probabilistic partial least squares

††journal: Journal of Multivariate Analysis

http://creativecommons.org/licenses/by-nc-nd/4.0/

1 Introduction

With the exponentially growing volume of data sets, multivariate methods for reducing dimensionality are an important research area in statistics. For combining two data sets, Partial Least Squares (PLS) regression [28] is a popular dimension reduction method [1]. PLS decomposes variation in each data set in a joint part and a residual part. The joint part is a linear projection of one data set on the other that best explains the covariance between the two data sets. These projections are obtained by iterative algorithms, such as NIPALS [28].

Partial Least Squares is popular in chemometrics [3]. In this field, the focus is on development of algorithms with good prediction performance, while the underlying model is less important. For applications in life sciences, interpretation of parameter estimates is necessary to gain understanding of the underlying molecular mechanisms.

For interpretation, a model needs to be identifiable. A model is said to be unidentifiable if the model corresponds to more than one set of parameter values. For PLS, rotation of the parameters does not change the model [26]. Hence, PLS does not provide an identifiable model. By constraining the parameter space, identifiability can be obtained. This involves solving a challenging optimization problem, since PLS requires estimating a structured covariance matrix [19].

For many problems in life sciences the study design needs to be accounted for, and algorithmic approaches such as PLS cannot be applied. Hence, a probabilistic formulation is necessary. Since likelihood method provides asymptotic standard errors of parameter estimates, computer-intensive resampling procedures can be avoided.

Also for other dimension reduction techniques, probabilistic methods have been developed. In 1999, Tipping and Bishop [23] developed the Probabilistic Principal Component Analysis (PPCA), in order to deal with missing data and dependent samples. In 2005, Bach and Jordan [2] developed Probabilistic Canonical correlation analysis (PCCA). However, for both PPCA and PCCA the model parameters are not identifiable, since rotation of the parameters does not change the model [23, 2]. In addition, in 2015, simultaneous envelopes models have been developed [4] for ‘low-dimensional’ settings. Further, Probabilistic PLS Regression and Probabilistic PLS have been proposed [14, 30]. For all these approaches, the model parameters are not identifiable.

In this paper we propose the Probabilistic Partial Least Squares (PPLS) model and show that the model parameters are identifiable up to a sign. We propose to maximize the PPLS likelihood with an EM algorithm that decouples the likelihood into several factors involving distinct sets of parameters. In the M step, a constrained optimization problem is solved by using a matrix of Lagrange multipliers.

The rest of the paper is organized as follows: In Section 2 we develop the PPLS model and establish identifiability of the model parameters. We develop an efficient algorithm for estimating the PPLS parameters. In Section 3 we study the performance of the PPLS estimators via simulations. In Section 4 we illustrate the PPLS model with two data matrices from two cohorts. We finish with a discussion.

2 Model and estimation

2.1 The PPLS model

Let $x$ and $y$ be two random row-vectors of dimension $p$ and $q$ , respectively. The Probabilistic Partial Least Squares (PPLS) model describes the two random vectors in terms of a joint part and a noise part. The joint part consists of correlated latent vectors, denoted by $t$ and $u$ , while the noise part consists of isotropic normal random vectors referred to as $e$ , $f$ and $h$ . The dimension of $t$ and $u$ is denoted by $r$ . The PPLS model describing the relationship between $x$ , $y$ and the joint and noise parts is

[TABLE]

Specifically, $e=(e_{1},\ldots,e_{p})$ , $f=(f_{1},\ldots,f_{q})$ and $h=(h_{1},\ldots,h_{r})$ are independent with zero mean and referred to as noise variables. The distributions of $e$ , $f$ and $h$ are multivariate normal with positive definite covariance matrix proportional to the identity matrix,

[TABLE]

The latent vector $t=(t_{1},\ldots,t_{r})$ is an $r$ -dimensional multivariate normal vector with with zero mean and diagonal positive definite covariance matrix $\Sigma_{t}=\mathrm{diag}(\sigma^{2}_{t_{1}},\ldots,\sigma^{2}_{t_{r}})$ , so

[TABLE]

The matrix $B=\mathrm{diag}\left(b_{1},\ldots,b_{r}\right)$ is a diagonal matrix of size $r$ , containing regression coefficients of $u$ on $t$ . Finally $W$ ( $p\times r$ ) and $C$ ( $q\times r$ ) are parameter matrices, referred to as loadings. The PPLS model for the random $p$ -dimensional row-vector $x$ and random $q$ -dimensional row-vector $y$ is given in Eq. (1). Let $\theta$ be the parameters of the PPLS model, i.e.,

[TABLE]

The PPLS model and its parameters are formulated conditional on the value of the dimension of the latent space $r$ .

The PPLS model (1) assumes a multivariate normal distribution for the observable random vectors $x$ and $y$ . The covariance between $x$ and $y$ is modeled by the regression of the latent vector $u$ on $t$ . The distribution of $(x,y)$ is $\mathcal{N}\left(0,\Sigma\right)$ with density given, for ${x}\in\mathbb{R}^{p}$ and ${y}\in\mathbb{R}^{q}$ , by

[TABLE]

and covariance matrix

[TABLE]

This follows from the normality property and from computing the variances and covariances of the random vectors; see Appendix A for the details.

2.2 Identifiability of probabilistic PLS

To establish identifiability of the PPLS model, some assumptions about its parameters have to be made. First, we assume that $0<r<\min(p,q)$ . Second, we assume that the diagonal elements of $B$ are positive, $b_{k}>0$ for $k\in\{1,\ldots,r\}$ . This will not restrict the model, since $t_{k}b_{k}$ is equal to $-t_{k}b_{k}$ in distribution. To identify the order of the loading vectors, the elements of $(\sigma_{t_{k}}^{2}b_{k})_{k=1}^{r}$ are assumed to be strictly decreasing with $k$ . Finally, we assume that the loading matrices $W$ and $C$ are orthogonal, i.e., $W^{\top}W=C^{\top}C=I_{r}$ . Together with the diagonality of $\Sigma_{t}$ in (2), it implies identifiability of all parameters up to sign. This is shown in the following theorem.

Theorem 1.

Let $r$ be fixed such that $0<r<\min(p,q)$ . Let $(x,y)_{1}$ and $(x,y)_{2}$ be generated by the PPLS model (1) having covariance matrix $\Sigma_{1}$ and $\Sigma_{2}$ with underlying parameters $\theta_{1}$ and $\theta_{2}$ as defined in (3), respectively. Then $\Sigma_{1}=\Sigma_{2}$ implies that $W_{1}=W_{2}\Delta$ , $C_{1}=C_{2}\Delta$ for some diagonal matrix $\Delta$ with on the diagonal elements $\delta_{i}\in\{-1,1\}$ , for $i\in\{1,\ldots,r\}$ , and all other parameters in $\theta_{1}$ and $\theta_{2}$ are equal.

The formal proof is given in Appendix B. Identifiability up to sign can be represented by a diagonal orthogonal matrix, namely a diagonal matrix with diagonal elements in $\{-1,1\}$ . For example, taking the model for $x$ in (1), we may substitute $W$ by $WR_{S}$ and $t$ by $tR_{S}$ , where $R_{S}$ is a diagonal orthogonal matrix, and get

[TABLE]

Since $(R_{S})_{jj}^{2}=1$ and the distribution of $t_{j}$ and $-t_{j}$ is the same, the right-hand side reduces to the original model for $x$ in (1). Note that the PPLS model is not invariant under general rotation matrices. Take a general rotation matrix $R$ , then we still get

[TABLE]

since $RR^{\top}=I_{r}$ . Inspecting the covariance of $TR$ we see that $\mathrm{cov}(TR)=R^{\top}\Sigma_{t}R$ , which is not diagonal if $R$ is not diagonal, and violates the PPLS model assumption on $\Sigma_{t}$ in Eq. (2).

2.3 Estimating the parameters

Unlike the iterative PLS methods, we propose a simultaneous approach for estimating the parameters, while taking the constraints in the PPLS model into account. Given the number of PPLS components, $r$ , the log likelihood of an independent and identically distributed (iid) sample $(X,Y)=\{(X_{1},Y_{1})^{\top},\ldots,(X_{N},Y_{N})^{\top}\}^{\top}$ of size $N$ from $(x,y)$ is

[TABLE]

with $S=N^{-1}\sum_{i=1}^{N}(X_{i},Y_{i})^{\top}(X_{i},Y_{i})$ and $\Sigma$ as in Eq. (4). To ensure empirical identifiability, we assume that $r<N$ . Note that the data dimensionality $p$ and $q$ may be larger than $N$ . For estimation of $\theta$ , maximum likelihood is used.

The log likelihood (5) depends in a non-linear way on the theoretical covariance matrix $\Sigma$ , which contains the loadings and variances. Optimizing this function directly is a non-trivial task, especially in high dimensions (i.e. when $p$ and $q$ are large). However, the PPLS model allows for a more simple (but iterative) optimization approach. Indeed, the maximum likelihood estimates for $\theta$ are a least squares type solution if the latent variables $t$ and $u$ are observed, as the model for $x$ and $y$ in (1) involves known $t$ and $u$ . In contrast, knowing $\theta$ allows for reconstruction of $t$ and $u$ by computing their conditional means given $x$ and $y$ . Alternating these two scenarios is actually an Expectation-Maximization (EM) [5] algorithm, with observed data $(x,y)$ and missing data $(t,u)$ .

The EM algorithm

The joint distribution of the complete data $(x,y,t,u)$ can (with abuse of notation) be decomposed as

[TABLE]

This follows from

[TABLE]

The second equation is implied by the fact that $x$ and $y$ are independent given $t$ and $u$ . The first two factors in the right-hand side can be rewritten as $f(x|t,u)=f(x|t)$ and $f(y|t,u)=f(y|u)$ , since $x$ and $u$ are independent given $t$ , and $y$ and $t$ are independent given $u$ . The last factor can be rewritten as $f(u|t)f(t)$ , yielding Eq. (6). The logarithm of the first three factors in the product in (6) can be written as

[TABLE]

Denote by $L_{\rm Comp}=\ln f(X,Y,T,U)$ the complete data log-likelihood, and define

[TABLE]

where the expectation is taken conditional on the observed $X$ and $Y$ , and $\theta^{\prime}$ is a fixed current estimate of the parameters. By optimizing $Q$ over all allowed $\theta$ , we get a non-negative increase in the observed log-likelihood $L$ . Moreover, by iterating this process of taking the expectation and maximizing over $\theta$ , the estimates in general converge to a stationary point or, in particular, a (possibly local) maximum of $L$ [5, 29]. The expectation step calculates the conditional expectation of the missing data given the observed data given by $Q(\theta)$ , which may in general involve intractable integration. However, for the exponential family, in particular the multivariate normal family, the complete likelihood depends on the complete data only via the sufficient statistics (called $t(\mathbf{x})$ in [5]), which are given in terms of the first and second moments of the complete data for the multivariate normal distribution. Computing $Q(\theta)$ implies computing the expected first and second moment of the latent variables: $\mathrm{E}\left(T|X,Y,\theta\right)$ , $\mathrm{E}\left(T^{\top}T|X,Y,\theta\right)$ , $\mathrm{E}\left(U|X,Y,\theta\right)$ , $\mathrm{E}\left(U^{\top}U|X,X,\theta\right)$ and $\mathrm{E}\left(U^{\top}T|X,Y,\theta\right)$ ; see Appendix C for details. Moreover, the decomposition in (6) allows for optimization of $\mathrm{E}\{\ln f(X|T)\}$ , $\mathrm{E}\{\ln f(Y|U)\}$ and $\mathrm{E}\{\ln f(U|T)\}$ separately, while only considering parameters involved in each factor. Maximizing $Q$ over $\theta$ yields parameter estimates for the next iteration in the EM algorithm. This leads us to the following theorem.

Theorem 2.

Let $X$ and $Y$ be an observed data sample of size $N$ , generated according to the PPLS model (1). Let $r$ be fixed such that $0<r<\min(N,p,q)$ . The parameters in $\theta$ can be estimated with an EM algorithm, yielding the following iterative scheme in $k$ with given starting values for $k=0$ :

[TABLE]

where $L_{W}$ and $L_{C}$ are such that

[TABLE]

The proof for Theorem 2 and the expressions for the conditional expectations are given in Appendix C. Note the dependency of $W^{k+1}$ and $C^{k+1}$ on the matrices $L_{W}$ and $L_{C}$ . These matrices ensure orthogonality of $W^{k+1}$ and $C^{k+1}$ in each iteration:

[TABLE]

where $\tilde{T}=\mathrm{E}(T|X,Y,\theta^{k})$ . The exact forms of $L_{W}$ and $L_{C}$ are not unique. Two choices are the eigenvectors of $\mathrm{E}(T^{\top}|X,Y,\theta^{k})\,X\,X^{\top}\,\mathrm{E}(T|X,Y,\theta^{k})$ and the lower triangular matrix of $X^{\top}\,\mathrm{E}(T|X,Y,\theta^{k})$ in the Cholesky decomposition. Note that these two orthogonalization matrices are straightforward to calculate with standard linear algebra tools. Since the PPLS model is identifiable, all choices for $L_{W}$ and $L_{C}$ will lead to the same optimum as the iteration number $k$ tends to infinity.

Standard errors for PPLS

Asymptotic standard errors for maximum likelihood estimators are found by inverting the observed Fisher information matrix. Following the reasoning of [16], the observed information may be given by

[TABLE]

Here $S(\hat{\theta})=\nabla\lambda(\hat{\theta})$ and $B(\theta)=-\nabla^{2}\lambda(\hat{\theta})$ are the gradient and negative of the second derivative of the log likelihood $\lambda(\theta)$ respectively evaluated in the MLE $\hat{\theta}$ . The explicit form of the asymptotic covariance matrix of $w_{k}$ is given in Appendix D. The square root of the diagonal elements are the asymptotic standard errors for the corresponding loading estimates.

Finding the number of components $r$

Available approaches to determine the number of PPLS components $r$ are minimizing a cross-validated loss function [9], visually inspecting eigenvalues of a covariance matrix [17], and selecting the number of components needed to achieve a certain proportion of variance explained by the components. In this paper we apply the last approach.

The PLS and PPLS algorithms are available as R packages at github.com/selbouhaddani under repository OmicsPLS and PPLS, respectively.

3 Simulation study

To evaluate the performance of the PPLS estimates, a simulation study was conducted. The aim was (1) to investigate the performance of PPLS for various scenarios, (2) to evaluate robustness of the PPLS estimates against departures from the normality assumption, (3) to compare the performance of the loading estimates with other probabilistic approaches, and (4) to compare the asymptotic PPLS standard errors with the bootstrap standard errors.

The simulated data were generated according to the PPLS model (1). The number of components was chosen to be $3$ , both in the data generation and in the estimation. We considered combinations of small and large sample size ( $N\in\{50,500\}$ ), low and high dimensionality ( $p\in\{20,1000\}$ ), and small and large proportion of noise (denoted by $\alpha_{n}\in\{0.1,0.5\}$ ). The robustness of PPLS was evaluated by considering four different continuous and discrete distributions for the latent variables $t$ , $u$ , $e$ , $f$ and $h$ ; we chose the normal distribution, the $t$ distribution with two degrees of freedom, the Poisson distribution with rate 1, and the Binomial distribution with two trials and success probability 0.25. These distributions cover a wide range of characteristics typically observed in omics data: heavy tailed, skewed and/or discrete. The latent variables were scaled to have mean zero and variances as specified below. All scenarios are summarized in Table 1.

The true loading values per component were generated from the normal density function with parameters $\mu$ and $\sigma$ , denoted by $N(x;\mu,\sigma^{2})$ , as follows

[TABLE]

The second columns in $W$ and $C$ were orthonormalized with respect to the first columns, and the third columns were orthonormalized with respect to the first two columns; we used a Gram–Schmidt procedure for both operations. The elements of the diagonal matrix $B$ were set to $b_{k}=\mathrm{e}^{\ln(1.5)-3(k-1)/10}=(1.5,1.11,0.82)$ , for $\Sigma_{t}$ we chose $\sigma_{t_{k}}=\mathrm{e}^{-(k-1)/10}=(1,0.90,0.82)$ .

For comparing the parameter estimates with the true values $\theta$ , we computed the bias and the variance of the estimates. To deal with the identifiability up to sign, we multiplied each estimated loading vector by $-1$ if the inner product of the estimated loading vector and the true loading vector was negative. Moreover, we swapped columns in $W$ and $C$ to maintain the same ordering as the ordering in the true loadings. This was done to avoid inflation of the bias or variance due to a wrong sign or ordering of the individual components.

PPLS estimates were compared to PLS estimates (with orthogonal loadings, see [20] for an overview) for all scenarios above. For comparing PPLS with PPCA and PCCA, we constructed a ‘null model’, i.e., $B=0$ , as well as $B\neq 0$ . We used the same scenarios as above, but we only considered the normal distribution.

Regarding standard errors for PPLS loadings, we compared asymptotic standard errors (as in Section 2) and bootstrap standard errors [27]. One set of two data matrices $X$ and $Y$ was simulated from a PPLS model with $p=q=20$ normally distributed variables. Based on these data, asymptotic and bootstrap standard errors were calculated. The number of bootstrap replicates was 1000. Furthermore, simulation-based standard errors for the loadings (based on standard deviations over 1000 data matrices drawn from the PPLS model used to generate the original data) were included as reference. Low and high noise levels ( $\alpha_{n}=0.1$ resp. $\alpha_{n}=0.5$ ), and small, large and ‘extra large’ sample sizes ( $N=50$ , $N=500$ and $N=5000$ , respectively) were considered. In the ‘extra large’ sample size scenario, no simulation-based reference was calculated. The PPLS estimation algorithm was considered to be converged when either the log-likelihood increment was below $10^{-6}$ , or $10^{4}$ EM steps were made. For each scenario, $1000$ replicates are used.

3.1 Results

Results for the loadings

The biases and variances of the estimated first component $W_{1}$ for the low dimensional case for normally distributed latent variables are graphically depicted in Figure 1. A black dot represents the average estimated PPLS loading value across 1000 simulations, whereas the width of the black dashed vertical line equals two times the standard deviation across 1000 simulations. The red star and red dashed vertical line represent the average loading value and twice the standard deviation for the PLS estimates. The true loading values are represented by a step function with steps at each index $j\in\{1,\ldots,p\}$ . Results for other components and scenarios are included in the Online Supplement.

Comparing the estimates for the first loading component $W_{1}$ , a better performance of PPLS compared to PLS was observed in terms of bias. In all scenarios the bias of the PPLS estimators were about the same as or less than the bias of the PLS estimators. Both estimators showed larger bias towards zero for higher absolute loading values. The biases decreased with a larger sample size and lower noise level. The biases of both estimators were very similar across different distributions. In the scenario where there is 50% noise and few (50) samples the variance of the PPLS estimators tended to be slightly larger than the variances of the PLS estimators when the true loading values were larger. This was observed across all distributions. The variances of the PPLS estimates were about the same or lower than the PPLS estimates in all other scenarios, where either the noise level was less or more samples were available. For both PPLS and PLS estimators the variances tended to increase with higher loading values. The variances decreased with larger sample size and lower noise level. The variances of bots estimators were very similar across different distributions. For the loading component $C_{1}$ and their PLS and PPLS estimators the same conclusions were obtained.

For the second loading component $W_{2}$ (shown in the Online Supplement), the biases of the PPLS loading estimates were as good as, and often better than the PLS loading estimates, especially at lower values. In the scenarios of 50% noise and a small sample size ( $N=50$ ) the bias was slightly larger for PPLS estimators compared to PLS estimates when the loading values were larger. Both estimators showed larger bias towards zero for higher loading values. The biases decreased with a larger sample size and lower noise level. For all distributions, the biases of both estimators were very similar. The variances of the PPLS estimators were as good as or lower than the PLS estimators, except in the scenario in which both the noise level was high (50%) and the sample size was small (50). In this scenario the variances of the PPLS estimators were still lower if the true loading values were close to zero, but higher for higher loading values. For both PPLS and PLS estimators the variances tended to increase with higher loading values. The variances decreased with larger sample size and lower noise level. The variances of both estimators were very similar across different distributions. For the loading component $C_{2}$ and their PLS and PPLS estimators the same conclusions were obtained.

For the third loading components $W_{3}$ and $C_{3}$ (shown in the Online Supplement), the same observations were made as for the first loading components $W_{1}$ and $C_{1}$ , both for the biases as for the variances.

For the high and extra high-dimensional case, the same results were obtained for the loadings $W$ and $C$ . See the Online Supplement for more details.

With regard to the comparison of PPLS with PPCA and PCCA, PPLS performed better than PCCA and similar to PPCA in most scenarios. Details are given in the Online Supplement.

Results for the variance parameters

The performance of the estimators of the variance parameters $B$ , $\sigma_{t}$ , $\sigma_{e}$ , $\sigma_{f}$ and $\sigma_{h}$ were also evaluated, the results are shown in Figure 2. We did not compare with PLS as these model parameters are not present in the PLS framework. For sake of comparison, we calculated the relative biases and variances of the estimates with respect to the true corresponding parameter value. The biases of the PPLS estimators for all variance parameters were very small for large sample size ( $N=500$ ), regardless of the noise. For small sample size ( $N=50$ ), the first two diagonal elements of $B$ and $\Sigma_{t}$ were overestimated, while the last component was underestimated. The noise parameters $\sigma_{e}$ and $\sigma_{f}$ were underestimated in these scenarios, while the estimator for $\sigma_{h}$ was unbiased, except in combination with a low noise level (10%). The relative biases decreased slightly with lower noise level, except for the earlier mentioned $\sigma_{h}$ , and decreased more with larger sample size. The relative variances of the estimators of $B$ , $\Sigma_{t}$ and $\sigma_{h}$ were larger than the variances of the estimators of $\sigma_{e}$ and $\sigma_{f}$ . For $B$ , there was a slight increase in variance across the three components. The variances decreased slightly with lower noise and more with larger sample size. The variances slightly decreased in the scenario of high dimensionality and high noise level. The same observations were made across the different distributions.

Ordering of the loadings

We compared the ordering of the true loadings $W$ and $C$ with the ordering of the estimated loadings. This provides a proportion across the 1000 simulation replicates in which the ordering matched. In Table 2, the proportion of correct orderings of $W$ for the scenario with normally distributed latent variables is shown for different scenarios. It can be seen that the proportion of correct orderings tends to be lower with smaller sample size and with higher noise level. Moreover, if the sample size is small, the proportion of correct orderings is much lower with higher noise. A higher dimensionality has a slightly negative impact on the correct ordering proportion when the sample size is larger, but a positive impact in the small sample size scenario. Especially, when also the noise level is high, this can be considerable. The same observations were made for the other distributions. Exactly the same proportions were observed for the loadings $C$ .

Comparison of PPLS standard errors

The results for low noise level are shown in Figure 3. In all scenarios, the asymptotic standard errors were smaller than the bootstrap standard errors for nearly all loading elements. In particular, for high loading values the difference between asymptotic and bootstrap standard errors tended to be large. This difference decreased with larger sample size: In the ‘extra large’ sample size, the bootstrap and asymptotic standard errors had similar magnitude. Similar observations were made for other distributions. For details, see the Online Supplement.

4 Data analysis

To illustrate the Probabilistic Partial Least Squares model, we apply it to IgG glycan datasets. Glycans, in particular IgG glycans, play an important role in the innate immune system, as well as in cell signaling. IgG2 glycans are less abundant than IgG1 glycans and more difficult to measure. Therefore, by using the relationships between IgG1 and IgG2 glycans, the characteristics of IgG2 can be better estimated. Hence, we will use IgG1 as $X$ matrix, and IgG2 as $Y$ matrix.

In total, 40 IgG glycans were measured, of which $p=20$ are of subclass IgG1 and $q=20$ are of subclass IgG2. These data were measured in two cohorts (CROATIA_Korcula with 951 participants and CROATIA_Vis with 796 participants) [12]. The data matrices containing IgG1 and IgG2 glycan variables are denoted by $X_{m}$ and $Y_{m}$ , with $m\in\{1,2\}$ , where $m=1$ corresponds to CROATIA_Korcula and $m=2$ corresponds to CROATIA_Vis. We apply PPLS to IgG1 and IgG2 glycans in each cohort separately and compare results.

In Figure 4, heatmaps of the correlations within and between the IgG1 and IgG2 glycans are shown, from which it is clear that there are many highly positive correlations between and within IgG1 and IgG2 in each data set. The RV coefficient [18], which generalizes the squared Pearson correlation coefficient to two matrices, was about 0.60 and 0.45 for CROATIA_Korcula and CROATIA_Vis cohorts respectively.

To determine the number of latent variables to use, we considered the total amount of variance explained by the latent space relative to the total amount of variation in the data: $||T_{m}||/||X_{m}||$ and $||U_{m}||/||Y_{m}||$ for $m\in\{1,2\}$ . By using four components, at least 89% of the total variation in each of the matrices $X_{1}$ , $X_{2}$ , $Y_{1}$ and $Y_{2}$ was accounted for.

For both cohorts, we fitted the PPLS models using $r=4$ latent components. The amount of overlap in each cohort, estimated by $\mathrm{tr}\hat{\Sigma}_{x,y}/\mathrm{tr}\hat{\Sigma}_{y}$ given in (4), was 58% and 46% for CROATIA_Korcula and CROATIA_Vis cohorts, respectively. The PPLS loadings are inspected to identify which IgG glycans contribute most to this overlap. The estimated IgG1 loadings $w_{j,k}$ , $j\in\{1,\ldots,p\}$ and $k\in\{1,\ldots,4\}$ , for both cohorts and both subclasses are depicted in Figure 5. The first joint component is proportional to the average glycan, as all glycans get about the same loading value. The second joint component involves especially G0 and G2 glycan subtypes, in which they are negatively correlated. Inspection of the loading values for the third component shows contibutions of fucosylated and non-fucosylated glycan subtypes. In the fourth component a pattern of positive and negative loading values is visible regarding the presence and absence of bisecting GlcNAc, respectively. The large loading value for G1NS is remarkable. The same conclusions hold for IgG2, as the estimated loading values were very similar. It is interesting to note that the observed patterns within components potentially reflect enzymatic synthesis where monosaccharides are added to glycan substrates [22]. Additionally, similar patterns are seen reflecting the inflammatory characteristics of glycans in aging and several different diseases [13]. Finally, the observed loading patterns were strikingly similar for both cohorts.

5 Discussion

We proposed PPLS to model the covariance matrix of two data sets. Maximum likelihood estimators for the model parameters were derived by solving a constrained optimization problem, and the parameter loadings were shown to be identifiable up to sign. This property ensures that PPLS estimates are comparable across several studies.

Our simulation study showed that the PPLS estimators had good performance and lower bias compared to PLS. Most notably, the performance of PPLS was robust to misspecification of the distribution of the variables. A smaller sample size and high noise level had a negative effect on the accuracy of the estimates, but large loading values were still non-zero. Also, compared to Probabilistic CCA estimates, the PPLS estimates were less biased and more efficient. For high-dimensional data, PCCA estimates have larger bias and higher variance. This is likely to be caused by the unstable inverse sample covariance matrix calculated when using PCCA. Moreover, if the number of variables is larger than the sample size, PCCA estimates cannot be obtained. Therefore, especially in omics data analysis, PPLS provides more robust findings.

As an illustration of the PPLS model, we analyzed IgG glycomics data from two cohorts. The high correlations in the data (Figure 4) and the use of multiple cohorts illustrate the applicability of PPLS to facilitate combination of results derived from different experimental settings. We found that the estimated loading values were almost identical across the two cohorts (Figure 5).

When multiple cohorts are available, a meta-analysis on the parameter estimates can be performed. In ordinary regression models, this has been addressed for both low-dimensional [6] and high-dimensional [10] design matrices. When there is no access to all data, PPLS estimates can be combined by using standard meta-analysis approaches [6]. Such an approach requires that the PPLS parameter estimates are identifiable and asymptotically normally distributed. For the PLS framework, several approaches to combine estimates across cohorts were developed when there is access to all data. A group-PLS approach was considered [15] to incorporate several groups of subjects in the model. The authors showed that under certain assumptions this approach provided better predictions than a model without group effects. However their model is not identifiable and requires $N>p$ . Another method is based on weighted least squares to combine data from different studies with potentially different covariates [11]. An alternative method, when access to data is possible, is to estimate joint parts between the data sets and the studies simultaneously. This yields a joint space with variables that have high loading values in most studies. For example, in [25], a non-probabilistic approach is pursued in a least squares estimation method using PCA. Performing data integration across studies, while taking into account uncertainties within each study, is one of our topics for future research, and will lead to more powerful analysis of the IgG glycans across cohorts.

To assess the statistical significance of loadings, the probabilistic framework provides alternative approaches to jackknifing and bootstrapping [27]. The observed Fisher information matrix can be used to estimate standard errors for individual loading parameters. For small sample sizes, bootstrap approaches appears to better reflect the uncertainty of the parameters. For large enough sample sizes, the asymptotic standard errors are close to the simulation-based standard errors. Typically, in epidemiological studies, the sample size is large enough to use asymptotic standard errors.

In this paper we ignored the fact that different biological ‘omics’ measurements have different error structures. An extension of Partial Least Squares was proposed to correct for systematic variation (variation induced by latent variables uncorrelated to the other data set) in the data sets, named Two-Way Orthogonal PLS (O2PLS) [8, 24]. Such an extension can be pursued for PPLS by adding for both $X$ and $Y$ in (1) a set of independent latent variables multiplied by their loading parameters. We are currently working on exploring the possibilities of a Probabilistic O2PLS for heterogeneous data sets.

Acknowledgments

The authors would like to thank the Editor-in-Chief, the Associate Editor and the referees for their valuable comments and suggestions. This work has been supported by the European Union’s Seventh Framework Programme (FP7-Health-F5-2012) under grant agreement number 305280 (MIMOmics). The CROATIA_Vis and CROATIA_Korcula studies were funded by grants from the Medical Research Council (UK), European Commission Framework 6 project EUROSPAN (Contract No. LSHG-CT-2006-018947), FP7 contract BBMRI-LPC (grant No. 313010), Croatian Science Foundation (grant 8875) and the Republic of Croatia Ministry of Science, Education and Sports (216-1080315-0302). We would like to acknowledge the staff of several institutions in Croatia that supported the field work, including but not limited to The University of Split and Zagreb Medical Schools, Institute for Anthropological Research in Zagreb and the Croatian Institute for Public Health. Glycome analysis was supported by the European Commission HighGlycan (contract No. 278535), MIMOmics (contract No. 305280), HTP-GlycoMet (contract No. 324400), IntegraLife (contract No. 315997). The IgG glycan data are available upon request.

Appendix A Variances and covariances

The covariance matrix of $(x,y)$ is given in (4). First note that $\mathrm{var}(u)=\mathrm{var}(tB+h)=B^{2}\Sigma_{t}+\sigma_{h}^{2}I_{r}$ , then compute

[TABLE]

The covariances between the observed and latent variables are as follows

[TABLE]

See, e.g., [21] for more details.

Appendix B Identifiability of PPLS

For establishing identifiability of the PPLS model, we need to prove that if the distribution of $(x,y)$ is given, there is only one corresponding set of parameters yielding this distribution. Since $(x,y)$ follows a zero mean normal distribution, identifiability is equivalent to

[TABLE]

where $\Sigma,\tilde{\Sigma}$ is defined, through $\theta,\tilde{\theta}$ , in (4). The following lemma will be very useful in establishing identifiability.

Lemma 1.

(Singular Value Decomposition).*

Let $W,\tilde{W}$ be $p\times r$ and $C,\tilde{C}$ be $q\times r$ , all orthogonal matrices. Let $D,\tilde{D}$ be $r\times r$ diagonal with $r$ distinct positive elements on the diagonal. Then $WDC^{\top}=\tilde{W}\tilde{D}\tilde{C}^{\top}$ (B.1) implies $W=\tilde{W}\Delta$ , $C=\tilde{C}\Delta$ for some diagonal matrix $\Delta$ of size $r\times r$ with on the diagonal elements $\delta_{i}\in\{-1,1\}$ and $D=\tilde{D}$ .*

Proof.

Let $A_{1}=WDC^{\top}$ and $A_{2}=\tilde{W}\tilde{D}\tilde{C}^{\top}$ . Consider $A_{i}A_{i}^{\top}$ and $A_{i}^{\top}A_{i}$ , $i\in\{1,2\}$ . The assertion (B.1) then implies the following.

[TABLE]

Note that both $WD^{2}W^{\top}$ and $\tilde{W}\tilde{D}^{2}\tilde{W}^{\top}$ are eigenvalue decompositions, as $D^{2}$ and $\tilde{D}^{2}$ are diagonal and $W,\tilde{W}$ and $C,\tilde{C}$ are orthogonal. The spectral theorem for matrices [7] then implies that whenever the elements in $D^{2},\tilde{D}^{2}$ are distinct, the corresponding columns in $W,\tilde{W}$ and $C,\tilde{C}$ are equal up to multiplication with the same sign. We thus get $W=\tilde{W}\Delta$ , $C=\tilde{C}\Delta$ and $D=\tilde{D}$ . ∎

Using this lemma, we show identifiability of the off-diagonal block part of the covariance matrix as given in Eq. (4).

Lemma 2.

If for matrices $W,\tilde{W}$ , $C,\tilde{C}$ and diagonal $B,\tilde{B}$ and $\Sigma_{t},\tilde{\Sigma}_{t}$ , given as in the PPLS model, $W\Sigma_{t}BC^{\top}=\tilde{W}\tilde{\Sigma}_{t}\tilde{B}\tilde{C}^{\top}$ , then $W=\tilde{W}\Delta$ , $C=\tilde{C}\Delta$ and $\Sigma_{t}B=\tilde{\Sigma}_{t}\tilde{B}$ .

Proof.

Applying Lemma 1 with $D=\Sigma_{t}B$ and $\tilde{D}=\tilde{\Sigma}_{t}\tilde{B}$ gives the desired result, since $\Sigma_{t}B$ and $\tilde{\Sigma}_{t}\tilde{B}$ are diagonal matrices with distinct ordered elements. ∎

Given $\Sigma_{x,y}$ we can identify $W$ and $C$ up to sign and the product $\Sigma_{t}B$ . We now show that in particular also the individual parameters $\Sigma_{t}$ and $B$ are identified from the upper diagonal block matrix $\Sigma_{x}$ .

Lemma 3.

If for matrix $W$ , diagonal matrices $\Sigma_{t}$ and $\tilde{\Sigma}_{t}$ and positive numbers $\sigma_{e}^{2},\tilde{\sigma}_{e}^{2}$ , given as in the PPLS model, $W\Sigma_{t}W^{\top}+\sigma_{e}^{2}I_{p}=W\tilde{\Sigma}_{t}W^{\top}+\tilde{\sigma}_{e}^{2}I_{p}$ (B.2), then $\sigma_{e}=\tilde{\sigma}_{e}$ and $\Sigma_{t}=\tilde{\Sigma}_{t}$ .

Proof.

Suppose (B.2) holds. Since $r<p$ and $p>1$ , one can find a unit vector $w_{\perp}$ such that $W^{\top}w_{\perp}=0$ . Multiplying with such vector yields $\sigma_{e}^{2}w_{\perp}=\tilde{\sigma}_{e}^{2}w_{\perp}$ . Multiplying again with $w^{\top}_{\perp}$ yields $\sigma_{e}^{2}=\tilde{\sigma}_{e}^{2}$ . It follows that we can identify $\sigma_{e}^{2}$ . We can now reduce (B.2) to $W\Sigma_{t}W^{\top}=W\tilde{\Sigma}_{t}W^{\top}$ . Pre-multiplying with $W^{\top}$ and post-multiplying with $W$ on both sides yields $\Sigma_{t}=\tilde{\Sigma}_{t}$ . ∎

We have seen in Theorem 2 that we can identify $\Sigma_{t}B$ . Since we identified $\Sigma_{t}$ we get identifiability of $B$ . The remaining parameters $\sigma_{h}^{2}$ and $\sigma_{f}^{2}$ are now shown to be identified using the lower block diagonal $\Sigma_{y}$ .

Lemma 4.

If for matrices $C$ , $B$ , $\Sigma_{t}$ , $\sigma_{f}^{2},\tilde{\sigma}_{f}^{2}$ and $\sigma_{h}^{2},\tilde{\sigma}_{h}^{2}$ , given as in the PPLS model, the assertion $\Sigma_{y}=\tilde{\Sigma}_{y}$ holds, i.e.,

[TABLE]

then $\sigma_{f}^{2}=\tilde{\sigma}_{f}^{2}$ and $\sigma_{h}^{2}=\tilde{\sigma}_{h}^{2}$ .

Proof.

In Theorem 3 take $W$ equal to $C$ , $\sigma_{e}^{2}$ equal to $\sigma_{f}^{2}$ , $\tilde{\sigma}_{e}^{2}$ equal to $\tilde{\sigma}_{f}^{2}$ , and the diagonal covariance matrices $\Sigma_{t}$ and $\tilde{\Sigma}_{t}$ equal to $\Sigma_{t}B^{2}+\sigma_{h}^{2}I_{p}$ and $\Sigma_{t}B^{2}+\tilde{\sigma}_{h}^{2}I_{p}$ . We find that we can identify $\Sigma_{t}B^{2}+\sigma_{h}^{2}$ and $\sigma_{f}^{2}$ . Since we already identified $\Sigma_{t}$ and $B$ , we have also identifiability of $\sigma_{h}^{2}$ . ∎

We conclude with the proof of Theorem 1.

Proof.

Suppose $\Sigma=\tilde{\Sigma}$ . This is true if and only if

[TABLE]

Applying Lemma 2 to the first equation, we identify $W$ and $C$ up to sign. Considering Lemma 3 together with Lemma 2, the second equation implies identifiability of $\Sigma_{t}$ , $B$ and $\sigma_{e}$ . The three Lemmas 2, 3 and 4 together with the last equation imply identifiability of $\sigma_{h}$ and $\sigma_{f}$ .

∎

Appendix C An Expectation-Maximization algorithm for PPLS

To obtain parameter estimates in the PPLS model, maximum likelihood is used. The EM algorithm is an iterative procedure for maximizing the observed log-likelihood (5) and consists of an Expectation step and a Maximization step. The following Lemma is convenient to make the expectation step explicit.

Lemma 5.

Let the pair $(z,x)$ be jointly multivariate normal row vectors with zero mean and covariance matrix

[TABLE]

Then $z|x$ is normally distributed with conditional mean $\mathrm{E}\left(z|x\right)=x\,\Sigma_{x}^{-1}\,\Sigma_{x,z}$ , and conditional covariance matrix $\mathrm{var}\left(z|x\right)=\Sigma_{z}-\Sigma_{z,x}\,\Sigma_{x}^{-1}\,\Sigma_{x,z}$ . Secondly, if $z=(t,u)$ , $\mathrm{cov}(t,x)=\Sigma_{t,x}$ and $\mathrm{cov}(x,u)=\Sigma_{x,u}$ , then the conditional covariance between $t$ and $u$ is $\mathrm{cov}(t,u|x)=\mathrm{cov}(t,u)-\Sigma_{t,x}\,\Sigma_{x}^{-1}\,\Sigma_{x,u}$ .

Proof.

The proof for the first part of the Lemma is found in [21]. The second part follows from the off diagonal blocks of $\mathrm{var}(z|x)$ . ∎

Expectation

The conditional first moments can be obtained by applying Lemma 5 while substituting $t$ or $u$ for $z$ and $(x,y)$ for $x$ .

[TABLE]

The same substitution can be made for the conditional second moments. Using $\mathrm{E}(a^{\top}b|z)=\mathrm{cov}(a,b|z)+\mathrm{E}(a|z)^{\top}\mathrm{E}(b|z)$ , we get

[TABLE]

where $S$ is the biased sample covariance matrix of $(x,y)$ . The conditional cross term equals

[TABLE]

The covariances are given by

[TABLE]

Although the the conditional expectations involve random variables and parameters, in the maximization step the calculated quantities are considered fixed and known.

Maximization

The maximization step involves maximizing the complete-data likelihood (6), we have seen that it can be decomposed in distinct factors. This allows optimization of the expected complete data likelihood to be split into four sub-maximizations, given by the individual factors and their respective parameters in the following annotated product:

[TABLE]

Moreover, it will become apparent that each parameter within each component can be decoupled, yielding a separate maximization per component per parameter. We focus on the part of $f(x|t)$ that depends on $W$ , which is given by

[TABLE]

To take into account the constraints on $W$ , namely $W^{\top}W=I_{r}$ , we introduce a matrix of Lagrange multipliers $\Lambda$ . We get as objective function

[TABLE]

Differentiating with respect to $W$ yields $2X^{\top}\mu_{t}-2WC_{TT}-2W\Lambda=2W\left(C_{TT}+\Lambda\right)-2X^{\top}\mu_{t}$ . One may choose $\Lambda$ so that $C_{TT}+\Lambda$ is invertible. In a maximum $W$ then satisfies $W=X^{\top}\mu_{t}\left(C_{TT}+\Lambda\right)^{-1}$ . We want to find a $\Lambda$ such that the constraint holds, i.e.,

[TABLE]

The last identity can be recognized as a Cholesky or Eigenvalue decomposition.

[TABLE]

with $L_{t}$ the lower triangular matrix of a Cholesky decomposition of $\mu_{t}^{\top}XX^{\top}\mu_{t}$ . Note that $L_{t}$ exists, since $\mu_{t}^{\top}XX^{\top}\mu_{t}$ is always positive semi-definite. Choosing $\Lambda=L_{t}^{\top}-C_{TT}$ , we get as update $W=X^{\top}\mu_{t}(L_{t}^{\top})^{-1}$ . Following the same reasoning, we obtain for the $f(Y|U)$ part $C=Y^{\top}\mu_{u}(L_{u}^{\top})^{-1}$ , where $L_{u}$ is the lower triangular matrix from the Cholesky decomposition of $\mu_{u}^{\top}YY^{\top}\mu_{u}$ .

The parameter $B$ involves maximizing $\ln f(U|T)$ , which is given by

[TABLE]

Taking the conditional expectation with respect to $(x,y)$ yields $-\mathrm{tr}\,\mathrm{E}(U^{\top}U-2U^{\top}TB+BT^{\top}TB|X,Y)$ . Differentiating with respect to $B$ and equating to the zero matrix yields

[TABLE]

To incorporate the constraint that $B$ should be diagonal, we set the diagonal elements to zero, yielding

[TABLE]

with $\circ$ the element-wise (Hadamard) product operator.

For the covariance matrix of $\Sigma_{t}$ , we consider $\ln f(T)$ which is given by

[TABLE]

After taking the conditional expectation of the last expression, it can be differentiated with respect to $\Sigma_{t}^{-1}$ , which yields

[TABLE]

The last Hadamart product ensures $\Sigma_{t}$ is diagonal.

To maximize over $\sigma_{e}^{2}$ , we consider $\ln f(X|T)$ and note that $E=X-TW^{\top}$ . Then $\ln f(X|T)$ is given by

[TABLE]

After taking the conditional expectation of the last expression, we differentiate it with respect to $\sigma_{e}^{-2}$ , yielding

[TABLE]

The same derivation can be applied to $\ln f(y|u)$ and $\ln f(u|t)$ to find

[TABLE]

Appendix D Asymptotic standard errors for PPLS loadings

Using notation as in [16] we define

[TABLE]

to be the part of the log likelihood depending on $w_{k}$ . We calculate the following first and second derivatives.

[TABLE]

We obtain

[TABLE]

Here $\mu_{k}=\mathrm{E}(t_{k}|X,Y)$ and $\sigma_{k}=\mathrm{E}(t_{k}^{\top}t_{k}|X,Y)$ . For explicit expressions of these expectations, see Appendix C. For the second derivative we get $\mathrm{E}\{B(w_{k})|X,Y\}=\sigma_{k}^{2}I_{p}/\sigma_{e}^{2}$ . The observed Fisher information is now

[TABLE]

and the asymptotic covariance matrix of $w_{k}$ is $-I_{\rm obs}^{-1}$ . The square root of the diagonal elements are the standard errors of the corresponding loading elements.

References

[1]

H. Abdi, Partial least squares regression and projection on latent structure regression (PLS Regression), Wiley Interdiscip. Rev. Comput. Stat. 2 (2010) 97–106.

[2]

F.R. Bach, M. Jordan, A Probabilistic Interpretation of Canonical Correlation Analysis, Dept. Stat. Univ. California, Berkeley, CA, Tech. Rep 688 (2005) 1–11.

[3]

A.L. Boulesteix, K. Strimmer, Partial least squares: A versatile tool for the analysis of high-dimensional genomic data, Br. Bioinform 8 (2007) 32–44.

[4]

R.D. Cook, X. Zhang, Simultaneous envelopes for multivariate linear regression, Technometrics 57 (2015) 11–25.

[5]

A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc. Ser. B 39 (1977) 1–38.

[6]

R. DerSimonian, N. Laird, Meta-analysis in clinical trials, Control. Clin. Trials 7 (1986) 177–188.

[7]

M. Eaton, Multivariate Statistics: A Vector Space Approach, Wiley, New York, 1983.

[8]

S. el Bouhaddani, J. Houwing-Duistermaat, P. Salo, M. Perola, G. Jongbloed, H.W. Uh, Evaluation of O2PLS in Omics data integration, BMC Bioinformatics 17.

[9]

S. Geisser, Predictive Inference, Taylor & Francis, London, 1993.

[10]

Q. He, H.H. Zhang, C.L. Avery, D.Y. Lin, Sparse meta-analysis with high-dimensional data, Biostatistics 17 (2016) 205–220.

[11]

X. Huang, W. Pan, X. Han, Y. Chen, L.W. Miller, J. Hall, Borrowing information from relevant microarray studies for sample classification using weighted partial least squares, Comput. Biol. Chem. 29 (2005) 204–211.

[12]

G. Lauc, J.E. Huffman, M. Pučić, L. Zgaga, B. Adamczyk, A. Mužinić, M. Novokmet, O. Polašek, O. Gornik, J. Krištić, T. Keser, V. Vitart, B. Scheijen, H.W. Uh, M. Molokhia, A.L. Patrick, P. McKeigue, I. Kolčić, I.K. Lukić, O. Swann, F.N. van Leeuwen, L.R. Ruhaak, J.J. Houwing-Duistermaat, P.E. Slagboom, M. Beekman, A.J.M. de Craen, A.M. Deelder, Q. Zeng, W. Wang, N.D. Hastie, U. Gyllensten, J.F. Wilson, M. Wuhrer, A.F. Wright, P.M. Rudd, C. Hayward, Y. Aulchenko, H. Campbell, I. Rudan, Loci Associated with N-Glycosylation of Human Immunoglobulin G Show Pleiotropy with Autoimmune Diseases and Haematological Cancers, PLoS Genet. 9 (2013) e1003225.

[13]

G. Lauc, M. Pezer, I. Rudan, H. Campbell, Mechanisms of disease: The human N-glycome, Biochim. Biophys. Acta 1860 (2016) 1574–82.

[14]

S. Li, J.O. Nyagilo, D.P. Dave, W. Wang, B. Zhang, J. Gao, Probabilistic partial least squares regression for quantitative analysis of Raman spectra, Int. J. Data Min. Bioinform. 11 (2015) 223–243.

[15]

Y. Li, P. Uden, D. von Rosen, A two-step PLS inspired method for linear prediction with group effect, Sankhyā Indian J. Stat. 75 (2013) 96–117.

[16]

T.A. Louis, Finding the observed information matrix when using the EM algorithm, J. R. Stat. Soc. Ser. B 44 (1982) 226–233.

[17]

K.V. Mardia, J.T. Kent, J.M. Bibby, Multivariate analysis, Academic Press, London, 1979.

[18]

P. Robert, Y. Escoufier, A unifying tool for linear multivariate statistical methods: The rv-coefficient, J. R. Stat. Soc. Ser. C 25 (1976) 257–265.

[19]

B. Roś, F. Bijma, J.C. de Munck, M.C.M. de Gunst, Existence and uniqueness of the maximum likelihood estimator for models with a Kronecker product covariance structure, J. Multivariate Anal. 143 (2016) 345–361.

[20]

R. Rosipal, N. Krämer, Overview and Recent Advances in Partial Least Squares, in: Proc. 2005 Int. Conf. Subspace, Latent Struct. Featur. Sel., Vol. 3940 of SLSFS’05, Springer, New York, 2006, pp. 34–51.

[21]

G.A.F. Seber, A.J. Lee, Linear Regression Analysis, 2nd ed., Wiley, Hoboken, NJ, 2003.

[22]

N. Taniguchi, K. Honke, M. Fukuda, H. Narimatsu, Y. Yamaguchi, T. Angata (Eds.), Handbook of Glycosyltransferases and Related Genes, Springer Japan, Tokyo, 2014.

[23]

M.E. Tipping, C.M. Bishop, Probabilistic principal component analysis, J. R. Stat. Soc. Ser. B 61 (1999) 611–622.

[24]

J. Trygg, S. Wold, O2-PLS, a two-block (X–Y) latent variable regression (LVR) method with an integral OSC filter, J. Chemom. 17 (2003) 53–64.

[25]

K. Van Deun, A.K. Smilde, M.J. van der Werf, H.A.L. Kiers, I. Van Mechelen, A structured overview of simultaneous component based data integration, BMC Bioinformatics 10 (2009) 246.

[26]

H. Wang, Q. Liu, Y. Tu, Interpretation of partial least-squares regression models with VARIMAX rotation, Comput. Stat. Data Anal. 48 (2005) 207–219.

[27]

R. Wehrens, W.E. van der Linden, Bootstrapping principal component regression models, J. Chemom. 11 (1997) 157–171.

[28]

H. Wold, Nonlinear iterative partial least squares (NIPALS) modelling: some current developments, in: Multivar. Anal. III (Proc. Third Internat. Symp. Wright State Univ., Dayton, Ohio, 1972), Academic Press, New York, 1973, pp. 383–407.

[29]

C.F.J. Wu, On the convergence properties of the EM algorithm, Ann. Stat. 11 (1983) 95–103.

[30]

J. Zheng, Z. Song, Z. Ge, Probabilistic learning of partial least squares regression model: Theory and industrial applications, Chemom. Intell. Lab. Syst. 158 (2016) 80–90.

Bibliography30

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] H. Abdi, Partial least squares regression and projection on latent structure regression (PLS Regression), Wiley Interdiscip. Rev. Comput. Stat. 2 (2010) 97–106.
2[2] F.R. Bach, M. Jordan, A Probabilistic Interpretation of Canonical Correlation Analysis, Dept. Stat. Univ. California, Berkeley, CA, Tech. Rep 688 (2005) 1–11.
3[3] A.L. Boulesteix, K. Strimmer, Partial least squares: A versatile tool for the analysis of high-dimensional genomic data, Br. Bioinform 8 (2007) 32–44.
4[4] R.D. Cook, X. Zhang, Simultaneous envelopes for multivariate linear regression, Technometrics 57 (2015) 11–25.
5[5] A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc. Ser. B 39 (1977) 1–38.
6[6] R. Der Simonian, N. Laird, Meta-analysis in clinical trials, Control. Clin. Trials 7 (1986) 177–188.
7[7] M. Eaton, Multivariate Statistics: A Vector Space Approach, Wiley, New York, 1983.
8[8] S. el Bouhaddani, J. Houwing-Duistermaat, P. Salo, M. Perola, G. Jongbloed, H.W. Uh, Evaluation of O 2PLS in Omics data integration, BMC Bioinformatics 17.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Probabilistic partial least squares model:

Abstract

keywords:

1 Introduction

2 Model and estimation

2.1 The PPLS model

2.2 Identifiability of probabilistic PLS

Theorem 1**.**

2.3 Estimating the parameters

The EM algorithm

Theorem 2**.**

Standard errors for PPLS

Finding the number of components rrr

3 Simulation study

3.1 Results

Results for the loadings

Results for the variance parameters

Ordering of the loadings

Comparison of PPLS standard errors

4 Data analysis

5 Discussion

Acknowledgments

Appendix A Variances and covariances

Appendix B Identifiability of PPLS

Lemma 1**.**

Proof.

Lemma 2**.**

Proof.

Lemma 3**.**

Proof.

Lemma 4**.**

Proof.

Proof.

Appendix C An Expectation-Maximization algorithm for PPLS

Lemma 5**.**

Proof.

Expectation

Maximization

Appendix D Asymptotic standard errors for PPLS loadings

References

Theorem 1.

Theorem 2.

Finding the number of components $r$

Lemma 1.

Lemma 2.

Lemma 3.

Lemma 4.

Lemma 5.