BIN-3-Way-PARAFAC-PLS: A 3-Way Partial Least Squares for Binary Response

1. Introduction

The prediction of a dependent variable (y) based on a set of multidimensional independent variables (x) has attracted significant interest across various fields. A commonly used approach for this prediction is Partial Least Squares (PLS) [1,2], which was originally developed for continuous data. PLS identifies a set of linear combinations of the predictor variables that exhibit the highest correlation with the response variables.

Multiway data, also known as higher-order data or tensors, has become increasingly common in various areas of our information society, such as bioinformatics [3], chemometrics [4], neuropsychology [5] and personality research [6]. In the context of multiway data, a ’way’ represents one of the dimensions along which the data is structured. For example, in a dataset with three dimensions, these might be subjects, variables and occasions. The term ’mode’ is used to specify the type of entity associated with each dimension. Typically, the number of ways and modes are equal, as in a three-dimensional dataset where each way corresponds to a different entity type (e.g. subjects, variables and occasions). However, there are scenarios where the number of modes is less than the number of ways. An example is a two-mode, three-way dataset where one entity type occurs in more than one dimension. This occurs, for example, when objects are compared across multiple subjects, resulting in a data structure where the same objects are repeated across dimensions [7,8]. To analyze these types of data, traditional methods based on Principal Component Analysis (PCA) [9] have been extended to accommodate the complexity of higher-dimensional structures. Techniques such as Tucker decomposition [10,11,12] and PARAFAC decomposition [13,14,15] are among the most widely used. Tucker decomposition generalizes PCA by decomposing the data into a core tensor and factor matrices, with each factor matrix corresponding to one of the modes of the tensor [16]. On the other hand, PARAFAC decomposes the tensor into a sum of component tensors, where each component is a rank-1 tensor.

Given the prevalence of this type of data, Partial Least Squares (PLS) has been extended to accommodate such data structures, facilitating the prediction of $\mathbf{Y}$ from $\mathbf{X}$, where $\mathbf{X}$ and/or $\mathbf{Y}$ are three-way data arrays. Specifically, [17] first combined the PARAFAC three-way analysis method with the PLS method. This method, called multilinear PLS (N-PLS), is thus a mixture of a trilinear model (PARAFAC) and PLS [18]. Subsequently, [19] propose a modification of the above model that aims not to improve the prediction of the dependent variables, but rather to improve the modeling of the independent variables. To achieve this, they rely on the distinction between low-rank and subspace approximation inherent in the PARAFAC and Tucker3 models. More recently, Sparse N-Way Partial Least Squares with L1 penalization has been proposed [20,21], which aims to achieve lower prediction errors by filtering out variable noise, while improving the interpretability and utilization of the N-PLS method.

The methods discussed earlier are designed for scenarios where both predictors and responses are continuous. When the response matrix contains binary variables, Partial Least Squares Discriminant Analysis (PLS-DA) [22] is often employed, which essentially fits a PLS regression to a dummy or fictitious variable. [23] introduced a PLS model for a single binary response, analogous to the PLS-1 model. [24] proposed the PLS-BLR method, a Partial Least Squares method for situations where the responses are binary and the predictors are continuous. Subsequently, they extended this approach to address cases where both the predictors and the responses are binary [25], introducing additional adaptations to handle this more complex scenario. This method, known as Binary Partial Least Squares Regression (BPLSR), further broadens the applicability of Partial Least Squares methods in analyzing binary data structures. This approaches extend PLS generalized linear regression [23] from a single response variable (PLS1) to multiple responses (PLS2), introducing dimension reduction, a feature not present in the original method. In the context of genomic data, other studies [26,27] applied Iteratively Re-weighted Least Squares (IRWLS) to handle binary responses with numerical predictors.

This paper proposes an extension of the N-PLS method to deal with scenarios where there is a set of binary dependent variables alongside a three-way data set of explanatory variables. It is therefore a combination of the previously mentioned N-PLS and PLS-BLR methods. Although the proposed BIN-3-Way-PARAFAC-PLS method is currently designed and tested for three-way data arrays, the underlying principles and techniques hold the potential for extension to N-way configurations. This limitation to three-way data is explicitly addressed throughout the manuscript to ensure clarity and accuracy. In addition, the use of biplot methods is proposed to graphically assess prediction quality and to identify explanatory variables that are more strongly associated with response variables. The biplot used will be a combination of logistic biplot [28,29] and interactive biplot [30].

The paper is structured as follows: Section 2 provides a brief overview of the PLS methods, followed by the introduction of the Bin-3-Way-PARAFAC-PLS method in Section 3, which allows PLS analysis of a matrix of binary responses alongside a three-way data array of explanatory variables. An algorithm for calculating model parameters is outlined in this section. Section 4 introduces a biplot method for visual assessment of the results obtained with Bin-3-Way-PARAFAC-PLS. Section 5 illustrates the application of the proposed method, and finally Section 6 offers concluding remarks and outlines potential future research directions.

2. A Brief Overview of PLS Models for Continuous Data

In this work, scalars are denoted by italic lowercase, vectors by bold lowercase, matrices by bold uppercase, and three-way data arrays by underlined bold uppercase. The letters I, J, K, L are used to denote the dimension of different orders. The abbreviations for the different PLS models are Bi/Tri-PLS-1/2/3, where the prefix (Bi/Tri) indicates the order of $\mathbf{X}$ and the last number (1/2/3) indicates the order of $\mathbf{Y}$, e.g. Bi-PLS3 refers to PLS with $\mathbf{X}$ as a data matrix and $\underline{\mathbf{Y}}$ as an 3-way data array.

Lets begin with a brief overview of bi-PLS1. Let be $\mathbf{X}\left(I\times J\right)$ and $\mathbf{y}\left(I\times 1\right)$ column centered (and scaled if necessary) matrices. The bi-PLS1 algorithm consists of two main steps. In the first step for each component, a rank one model of $\mathbf{X}$ and y is constructed. Then in a second step, these models are subtracted from $\mathbf{X}$ and $\mathbf{y}$, resulting in residuals. From these residuals a new set of components is obtained. It is therefore a question of finding a one-component model of $\mathbf{X}$ such that ${\widehat{x}}_{ij}=t_i{w}_j$ where $t_i$ and $w_j$ are the scores and the weights respectively. This is done by maximizing the following function

$$\underset{w}{max}cov\left(\mathbf{t},\mathbf{y}\right)$$

subject to $\mathbf{X}\mathbf{w}=\mathbf{t}\wedge ∥\mathbf{w}∥=1$ or equivalently

$$\underset{w}{max}\sum _{i=1}^I{t}_i{y}_i$$

subject to $t_i={\sum }_{j=1}^J{x}_{ij}{w}_j\wedge ∥\mathbf{w}∥=1$

It can be verified that the value of the above expression is maximal when

$$\mathbf{w}={\displaystyle \frac{{\mathbf{X}}^T\mathbf{y}}{∥{\mathbf{X}}^T\mathbf{y}∥}}$$

or, in other words, if $\mathbf{w}$ is the left singular vector of ${\mathbf{X}}^T\mathbf{y}$. Once the first component is obtained, $\mathbf{y}$ and $\mathbf{X}$ are deflated and the second component is sought following the same process.

Similarly, when $\underline{\mathbf{X}}$ is a three-way data array $\left(I\times J\times K\right)$, the goal of the algorithm is to maximize the covariance between $\underline{\mathbf{X}}$ and $\mathbf{y}$. Let $\mathbf{X}$ be the $\left(I\times JK\right)$ unfolded version [31] of $\underline{\mathbf{X}}$, then it is a matter of finding latent spaces ${\mathbf{W}}^J$$\left(J\times S\right)$ and ${\mathbf{W}}^K$$\left(K\times S\right)$ that maximise the covariance between $\mathbf{X}$ and $\mathbf{y}$, so that $\mathbf{X}$ can be written as:

$$\mathbf{X}=\mathbf{T}{\left({\mathbf{W}}^J\odot {\mathbf{W}}^K\right)}^T+{\mathbf{E}}_X=\mathbf{T}{\mathbf{P}}^T+{\mathbf{E}}_X$$

(1)

where $\mathbf{T}$$\left(I\times S\right)$ is the matrix of scores in the first mode, ${\mathbf{W}}^J$ and ${\mathbf{W}}^K$ are the weights of the second and third mode respectively, ⊙ is the Khatri-Rao product (column-wise Kronecker product) [32], ${\mathbf{E}}_X\left(I\times JK\right)$ is a matrix of residuals and S is the number of components. The above decomposition corresponds to the PARAFAC model for three-way arrays. In this case the function to be maximized is

$$\begin{array}{cc}\hfill & \underset{{\mathbf{w}}^J{\mathbf{w}}^K}{max}\left(\sum _{i=1}^I{t}_i{y}_i\mid t_i=\sum _{j=1}^J\sum _{k=1}^K{x}_{ijk}{w}_j^J{w}_k^K\wedge ∥{\mathbf{w}}^J∥=∥{\mathbf{w}}^K∥=1\right)=\hfill \\ \hfill & \underset{{\mathbf{w}}^J{\mathbf{w}}^K}{max}\left(\sum _{i=1}^I\sum _{j=1}^J\sum _{k=1}^K{y}_i{x}_{ijk}{w}_j^J{w}_k^K\right)=\hfill \\ \hfill & \underset{{\mathbf{w}}^J{\mathbf{w}}^K}{max}\left(\sum _{j=1}^J\sum _{k=1}^K{z}_{jk}{w}_j^J{w}_k^K\right)\hfill \end{array}$$

where ${\mathbf{w}}^J$ and ${\mathbf{w}}^K$ can be obtained using the singular value decomposition (SVD) of $\mathbf{Z}$, being $\mathbf{Z}$ the matrix obtained as $z_{jk}={\sum }_i{y}_i{x}_{ijk}$.

Given the number of components S, the full algorithm for the Tri-PLS1 model is as follows:

Algorithm 1 Tri-PLS1 for S components

1: Center $\underline{\mathbf{X}}$ and $\mathbf{y}$ 2: ${\mathbf{y}}_{\mathbf{0}}=\mathbf{y}$ 3: for $s=1\to S$do 4:     Calculate ${\mathbf{Z}}_s$ 5:     Determine ${\mathbf{w}}_s^J$ and ${\mathbf{w}}_s^K$ by SVD of ${\mathbf{Z}}_s$ 6:     ${\mathbf{w}}_s={\mathbf{w}}_s^K\odot {\mathbf{w}}_s^J$                         ▹ Update ${\mathbf{w}}_s$ 7:     ${\mathbf{t}}_s=\mathbf{X}{\mathbf{w}}_s/{\mathbf{w}}_s^T{\mathbf{w}}_s$                                ▹ Update ${\mathbf{t}}_s$ 8:     ${\mathbf{b}}_s={\left({\mathbf{T}}^T\mathbf{T}\right)}^{-1}{\mathbf{T}}^T{\mathbf{y}}_0$; where $\mathbf{T}=[{\mathbf{t}}_1,{\mathbf{t}}_2\dots {\mathbf{t}}_s]$ 9:     $\mathbf{X}=\mathbf{X}-{\mathbf{t}}_s{\mathbf{w}}_s^T$                          ▹ Deflate: $\mathbf{X}$ 10:     $\mathbf{y}={\mathbf{y}}_0-{\mathbf{t}}_s{\mathbf{b}}_s$                         ▹ Deflate: $\mathbf{y}$ 11: end for

In case where there are several dependent variables $\mathbf{Y}$$\left(I\times L\right)$, $\mathbf{Y}$ can be expressed as:

$$\mathbf{Y}=\mathbf{U}{\mathbf{Q}}^T+{\mathbf{E}}_Y$$

(2)

where $\mathbf{U}\left(I\times S\right)$ is the matrix of scores, $\mathbf{Q}\left(L\times S\right)$ is the loading matrix and ${\mathbf{E}}_Y\left(I\times L\right)$ is a matrix of residuals.

The model is estimated in such a way that the covariance between $\mathbf{T}$ and $\mathbf{U}$ is maximized. The prediction model between $\mathbf{X}$ and $\mathbf{Y}$ can be expressed as:

$$\mathbf{U}=\mathbf{T}\mathbf{B}+{\mathbf{E}}_U$$

where matrix $\mathbf{B}\left(S\times L\right)$ represents the regression coefficients.

The algorithm that implements the Tri-PLS2 model for S components is as follows:

Algorithm 2:Tri-PLS2 for S components

1: Center $\underline{\mathbf{X}}$ and $\mathbf{Y}$ and unfold $\underline{\mathbf{X}}$ as $\mathbf{X}$ 2: for $s=1\to S$ do 3:     Let ${\mathbf{u}}_s$ equal a column in $\mathbf{Y}$                         ▹ Init: ${\mathbf{u}}_s$ 4:     repeat 5:         Calculate ${\mathbf{Z}}_s$ using $\mathbf{X}$ and ${\mathbf{u}}_s$ 6:         Determine ${\mathbf{w}}_s^J$ and ${\mathbf{w}}_s^K$ by SVD of ${\mathbf{Z}}_s$ 7:         ${\mathbf{w}}_{\mathbf{s}}={\mathbf{w}}_s^K\odot {\mathbf{w}}_s^J$                                ▹ Update: ${\mathbf{w}}_s$ 8:         ${\mathbf{t}}_{\mathbf{s}}=\mathbf{X}{\mathbf{w}}_s/{\mathbf{w}}_s^T{\mathbf{w}}_s$                                ▹ Update: ${\mathbf{t}}_s$ 9:         ${\mathbf{q}}_s={\mathbf{Y}}^T{\mathbf{t}}_s/∥{\mathbf{Y}}^T{\mathbf{t}}_s∥$                                ▹ Update: ${\mathbf{q}}_s$ 10:         ${\mathbf{u}}_s=\mathbf{Y}{\mathbf{q}}_s$                         ▹ Update: ${\mathbf{u}}_s$ 11:     until convergence 12:     ${\mathbf{b}}_s={\left({\mathbf{T}}^T\mathbf{T}\right)}^{-1}{\mathbf{T}}^T{\mathbf{u}}_s$; where $\mathbf{T}=[{\mathbf{t}}_1,{\mathbf{t}}_2\dots {\mathbf{t}}_s]$ 13:     $\mathbf{X}=\mathbf{X}-{\mathbf{t}}_s{\mathbf{w}}_s^T$                                ▹ Deflate: $\mathbf{X}$ 14:     $\mathbf{Y}=\mathbf{Y}-\mathbf{T}{\mathbf{b}}_s{\mathbf{q}}_s^T$                                ▹ Deflate: $\mathbf{Y}$ 15: end for

3. Bin-3-Way-PARAFAC-PLS

The previous linear model is not valid and a logit transformation must be performed if the response matrix is binary. In this case, the expected values of the binary response are probabilities and the logit function is used as the link function. Let be $E\left[\mathbf{Y}\right]$ = $\Pi$, then

$$logit(\Pi )=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{U}{\mathbf{Q}}^T$$

(3)

This equation is a generalization of Equation 2 except now it is needed a vector ${\mathbf{q}}_{\mathbf{0}}$ with intercepts for each variable, given that binary matrix cannot be centered in the same way as continuous matrix. Each probability ${\pi }_{il}$ can be expressed as

$${\pi }_{il}={\displaystyle \frac{e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}{1+e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}}$$

(4)

The procedure for estimating the parameters of the model is similar to the one described in [24]. The cost function in this case is as follows:

$$L=\sum _{i=1}^I\sum _{l=1}^L\left[-y_{il}log\left({\pi }_{il}\right)-(1-y_{il})log\left(1-{\pi }_{il}\right)\right]$$

(5)

The objective is to estimate the parameters $\mathbf{U}$, $\mathbf{Q}$ and ${\mathbf{q}}_0$ that minimize the value of the above cost function. There are no closed-form solutions for the optimization problem; therefore, an iterative approach is employed, where we generate a sequence of progressively smaller values for the loss function with each iteration. The gradient method is applied recursively, solving for one component at a time while keeping the others fixed. The update for each parameter is as follows:

$$q_{l0}=q_{l0}-\alpha {\displaystyle \frac{\partial L}{\partial q_{l0}}},u_{ls}=u_{ls}-\alpha {\displaystyle \frac{\partial L}{\partial u_{ls}}},q_{ls}=q_{ls}-\alpha {\displaystyle \frac{\partial L}{\partial q_{ls}}},$$

(6)

for a choice of $\alpha$.

The corresponding gradients are:

$${\displaystyle \frac{\partial L}{\partial q_{l0}}}=\sum _{i=1}^I({\pi }_{il}-y_{il}),{\displaystyle \frac{\partial L}{\partial u_{is}}}=\sum _{l=1}^L{q}_{ls}({\pi }_{il}-y_{il}),{\displaystyle \frac{\partial L}{\partial q_{ls}}}=\sum _{i=1}^I{u}_{is}({\pi }_{il}-y_{il})$$

(7)

The algorithm will iteratively compute parameters for rows ${\mathbf{u}}_s=(u_{1s},\dots ,u_{Is})$ and columns ${\mathbf{q}}_s=(q_{1s},\dots ,q_{Ls})$ for each dimension s, keeping fixed the parameters obtained for each previous dimension, with the aim of deriving uncorrelated components. Previously, the constants ${\mathbf{q}}_0=(q_{10},\dots ,q_{L0})$ must be determined separately.

In logistic regression, a common issue known as separation can arise when the data are perfectly separable by a linear combination of the predictors. In such cases, the maximum likelihood estimator fails to exist and diverges to infinity [33]. Even in situations of quasi-separation, where the separation is not complete, the maximum likelihood estimator can become highly unstable.

To address this problem, a typical approach is to employ a penalized likelihood method [34]. This technique incorporates a penalty term into the likelihood function, promoting smaller coefficient values and improving stability. In this work, we apply the Ridge penalty [35] to mitigate the effects of separation.

Combining the Tri-PLS2 algorithm with the previously discussed method for obtaining the components of a binary matrix, we would have the following algorithm:

Algorithm 3: Bin-Tri-PLS2 for S components

1: Choose $\alpha$ 2: ${\mathbf{q}}_{\left(0\right)}=random$                         ▹ Init: ${\mathbf{q}}_0$ 3: repeat 4:     $q_{l0}\leftarrow q_{l0}-\alpha {\displaystyle \sum _{i=1}^I}({\pi }_{il}-y_{il}),(l=1,\dots ,L)$                         ▹ Update: ${\mathbf{q}}_0$ 5:     ${\pi }_{il}\leftarrow {\displaystyle \frac{e^{q_{l0}}}{1+e^{q_{l0}}}};(i=1,\dots ,I;l=1,\dots ,L)$                         ▹ Update: $\Pi$ 6: until${\mathbf{q}}_0$ do not change 7: for$s=1\to S$do 8:     ${\mathbf{u}}_s\leftarrow random$                         ▹ Init: ${\mathbf{u}}_s$ 9:     repeat 10:         ${\mathbf{w}}_s^T$ = ${\mathbf{u}}_{\mathbf{s}}^T$$\mathbf{X}$/ ${\mathbf{u}}_{\mathbf{s}}^T$${\mathbf{u}}_{\mathbf{s}}$ 11:         Build ${\mathbf{Z}}_s$ by refolding ${\mathbf{w}}_s$ 12:         Determine ${\mathbf{w}}_s^J$ and ${\mathbf{w}}_s^K$ by SVD of ${\mathbf{Z}}_s$ 13:         ${\mathbf{w}}_{\mathbf{s}}={\mathbf{w}}_s^K\odot {\mathbf{w}}_s^J$                         ▹ Update: ${\mathbf{w}}_s$ 14:         ${\mathbf{t}}_{\mathbf{s}}=\mathbf{X}{\mathbf{w}}_s/{\mathbf{w}}_s^T{\mathbf{w}}_s$                         ▹ Update: ${\mathbf{t}}_s$ 15:         repeat 16:            $q_{ls}\leftarrow q_{ls}-\alpha {\displaystyle \sum _{i=1}^I}{u}_{is}({\pi }_{il}-y_{il})$                         ▹ Update: ${\mathbf{q}}_s$ 17:            ${\pi }_{il}={\displaystyle \frac{e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}{1+e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}}$                         ▹ Update: $\Pi$ 18:         until ${\mathbf{q}}_s$ do not change 19:         repeat 20:            $u_{is}\leftarrow u_{is}-\alpha {\displaystyle \sum _{l=1}^L}{q}_{ls}({\pi }_{il}-y_{il})$                         ▹ Update: ${\mathbf{u}}_s$ 21:            ${\pi }_{il}={\displaystyle \frac{e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}{1+e^{\left(q_{l0}+{\sum }_s{u}_{is}{q}_{ls}\right)}}}$                         ▹ Update: $\Pi$ 22:         until ${\mathbf{u}}_s$ do not change 23:         $L\leftarrow {\displaystyle \sum _{i=1}^I}{\displaystyle \sum _{l=1}^L}\left[-y_{il}log\left({\pi }_{il}\right)-(1-y_{il})log\left(1-{\pi }_{il}\right)\right]$ 24:     until L does not change 25:     $\mathbf{X}\leftarrow \mathbf{X}-{\mathbf{t}}_s{\mathbf{w}}_s^T$                         ▹ Deflate: $\mathbf{X}$ 26: end for

The regression coefficients of the response variables on the observed variables can be obtained as follows:

$$logit(\Pi )=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{U}{\mathbf{Q}}^T=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}\mathbf{B}{\mathbf{Q}}^T=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{X}\mathbf{P}\mathbf{B}{\mathbf{Q}}^T$$

(8)

then

$$logit(\Pi )=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{X}\mathbf{C}$$

(9)

with

$$\mathbf{C}=\mathbf{P}\mathbf{B}{\mathbf{Q}}^T$$

(10)

i.e., $\mathbf{C}$ are the regression coefficients relative to the observed variables.

4. Bin-3-Way-PARAFAC-PLS Biplot

In this section, a biplot representation is proposed for the graphical visualization of the results obtained by the Bin-3-Way-PARAFAC-PLS algorithm. It begins with a brief review of the biplot for continuous data to introduce the biplot for binary data, known as the logistic biplot. The combination of the PLS-biplot [36] with the logistic biplot [28,29] and the interactive biplot [30] for three-way data will result in the proposed Bin-3-Way-PARAFAC-PLS biplot.

4.1. Classical Biplot

The biplot was originally introduced by [37] and extended by [38] to visualise relationships between variables and observations in a multivariate data set. It has since been widely adopted in fields as diverse as ecology, chemometrics and the social sciences because it provides a comprehensive graphical representation of both observations and variables simultaneously in a low-dimensional space.

It is a well-known fact that any $I\times J$ matrix $\mathbf{X}$ with rank R can be decomposed into two matrices, $\mathbf{A}$ of size $I\times R$ and $\mathbf{B}$ of size $J\times R$, such that $\mathbf{X}=\mathbf{A}{\mathbf{B}}^T$. This decomposition enables the assignment of vectors (row markers) $a_1,\dots ,a_I$ to each of the I rows of $\mathbf{X}$ and vectors (column markers) $b_1,\dots ,b_J$ to each column of $\mathbf{X}$, which provides a representation of $\mathbf{X}$ using $I+J$ vectors in an R-dimensional space. In the case of matrices with rank two, this representation occurs in a two-dimensional space (plane), but for higher-rank matrices, an approximation in a reduced-dimensional space (either 2 or 3 dimensions) is obtained, such that $\mathbf{X}\approx \mathbf{A}{\mathbf{B}}^T$.

This factorization, usually achieved by Singular Value Decomposition (SVD), allows the representation of the matrix $\mathbf{X}$ in a space of reduced dimension such that $x_{ij}\approx {\mathbf{a}}_i^T{\mathbf{b}}_j$. Geometrically, this means that the projection of a row marker onto the vector representing a variable approximates the value of the individual on that variable and allows the individuals to be ordered by that variable. It is also possible to interpret the angles formed by two variables or the distances between individuals. The angles between variables in a biplot represent the correlations or associations between them. If two variables form a small angle, this indicates a strong positive correlation between them, while a large angle indicates a weaker or negative correlation. The distances between individuals (observations) in a biplot reflect their similarities or dissimilarities based on the values of the variables. Closer distances indicate greater similarity, while further distances indicate greater dissimilarity.

4.2. Logistic Biplot

Let $\mathbf{X}$ be a binary matrix and $\Pi =E\left[\mathbf{X}\right]$, we have that

$$logit(\Pi )=\mathbf{1}{\mathbf{b}}_0^T+\mathbf{A}{\mathbf{B}}^T$$

Apart from the constant vector ${\mathbf{b}}_0$, which is necessary because the binary matrix cannot be centred, this equation gives a biplot on a logit scale. The geometry of the logistic biplot is similar to that of the continuous biplot, with the distinction that it predicts probabilities. In this case it is necessary that the column vectors are accompanied by scales indicating the probabilities. Sometimes reduced prediction scales are used to simplify the presentation, where a dot is used to predict a probability of $0.5$ and an arrow is used to predict a probability of $0.75$. This approach not only reveals the direction in which the probability increases but also offers insights into discriminant power, with shorter arrows typically signifying lower discriminant power of the associated variable.

4.3. Interactive Biplot

If $\underline{\mathbf{X}}$ is a three-way data array $\left(I\times J\times K\right)$, it is possible to do a decomposition on it, let’s say a PARAFAC decomposition, so that the matricization of the first mode ${\mathbf{X}}_A$ can be written as ${\mathbf{X}}_A=\mathbf{A}{\left(\mathbf{C}\odot \mathbf{B}\right)}^T$. From this decomposition, it is possible to construct a biplot where the row and column markers are ${\mathbf{A}}^{*}=\mathbf{A}$ and ${\mathbf{B}}^{*}=\left(\mathbf{C}\otimes \mathbf{B}\right)$. In this biplot, each possible combination of variables from two of the modes is represented by a single marker. For example, as the second and third modes are combined, there will be $J\times K$ column markers. The interpretation of interactive biplots mirrors that of two-dimensional biplots, as they also provide a (close) approximation of the values within the data matrix. As the number of column markers, $J\times K$, can be large, making a single biplot of all of them can be unclear. Therefore, one of the modes is fixed and as many biplot representations are made as the size of that fixed mode. For example, if the second mode is fixed, there would be J biplot representations, each containing K variables.

4.4. Bin-3-Way-PARAFAC-PLS Biplot (Triplot)

The factorization of the matrix $\mathbf{X}$ in Equation 1 gives a biplot where the row and column markers are the rows of the matrices $\mathbf{T}$ and $\mathbf{P}=\left({\mathbf{W}}^J\odot {\mathbf{W}}^K\right)$ respectively. Furthermore, a logistic biplot of the binary matrix $\mathbf{Y}$ can be constructed from Equation 3, where the row markers also correspond to the rows of $\mathbf{T}$. Combining the two, there are three sets of markers: the individual scores in $\mathbf{T}$, the markers for the response variables in $\mathbf{Q}$ and ${\mathbf{q}}_0$, and the markers for the predictors in $\mathbf{P}$. The combined display of all these labels is called a triplot. In this case, individuals can be projected onto vectors corresponding to binary response variables and continuous explanatory variables, giving expected probabilities and values in each of the continuous variables.

In summary, the triplot allows the analysis of relationships between explanatory variables, the identification of which explanatory variables are more influential in predicting response variables, and the approximation of individual values in both types of variables. In the case of binary variables, this approximation represents the probability that the individual has such a response.

5. Illustrative Application

An important goal of psychology is to investigate the relationship between individuals’ reactions to specific situations and their inherent traits or dispositions [39]. It is plausible to assume that people’s reactions to certain scenarios can provide information about their personality traits. Therefore, we will perform a BIN-3-Way-PARAFAC-PLS analysis on a data set that includes, on the one hand, individuals’ reactions to certain situations and, on the other hand, their personality traits.

The hypothetical data set (which is a artificial data set) consists of a data matrix of 8 persons x 7 emotions x 6 situations, which captures the emotions people feel when faced with the situations considered, and a binary data matrix of 8 persons x 9 dispositions, which captures whether or not they exhibit certain personality traits. The situations considered are "Argue with someone", "Partner leaves you", "Someone lies about you", "Give a bad speech", "Fail an exam" and "Write a bad paper". The possible responses are "Other anger", "Shame", "Love", "Grief", "Fear", "Guilt" and "Self-anger". Finally, the personality traits considered are "Fear of rejection", "Kindness", "Importance of others’ judgements", "Altruism", "Neuroticism", "Being strict with oneself", "Low self-esteem", "Conscientiousness" and "Depression". The continuous data set has been appropriately centered and normalized. The aim of the analysis is to predict people’s personality traits from the way they react to certain situations, using the proposed Bin-Tri-PLS2 method.

The variance of the predictors explained by the decomposition is shown in Table 1.

It can be seen that the first two dimensions explain 67.01% of the variance of the predictors.

Observing the scores for the situation component (Table 2), the first component appears to represent an interpersonal aspect, as it is strongly associated with situations involving conflict or interactions with others, such as arguments, being falsely accused, or a partner leaving. On the other hand, the second component seems to reflect an intrapersonal dimension, as it is closely related to situations that involve personal failure or self-reflection, such as giving a poor presentation, failing an exam, or submitting a subpar paper.

Regarding the response component scores (Table 3), emotions such as other anger and shame have high negative loadings on the first component, while love has a high positive loading on the same component. On the other hand, guilt and self-anger score highly on the second component with negative values. Sorrow shows high negative loadings on both components.

Table 4 presents several fit measures for binary response variables, including a test comparing the model to the null model, three pseudo R-squared values, and the percentage of correctly classified cases. The Bin-3-Way-PARAFAC-PLS model achieved an overall classification accuracy of 81.94%.

It is quite clear that the Bin-3-Way-PARAFAC-PLS biplot is very difficult to interpret because of the number of arrows ($6\times 7$). It may be a good idea to create a different biplot for each situation including the responses and the personality traits. It is important to note that in this biplot, the probabilities of binary variables are depicted with a dot representing a probability of 0.5, and an arrow indicating a probability of 0.75.

Figure 1 shows the triplot of situations, responses, and dispositions, where the disposition and responses are represented as red and black vectors, respectively, and the individuals are shown as green points. It can be seen that individuals 5, 6, 7, and 8 share very similar personality traits ("Depression", "Conscientiousness", "Low self-esteem", "Being strict with oneself", and "Neuroticism"). These traits are also present in person 4, but to a lesser extent. Furthermore, it is evident that these personality traits are highly correlated. The remaining individuals are characterized by being kind, altruistic, attentive to the judgments of others, and fearful of rejection. Similarly, these personality traits are closely related to each other.

When analysing the reactions to the different situations, similarities can be found between the situations "Giving a bad speech", "Failing a test" and "Writing a bad paper". In all cases, people who react to these situations by feeling "Sorrow", "Self anger" or "Guilt" are characterised by the following personality traits: "Depression", "Conscientiousness", "Low self esteem" and "Being strict to oneself". On the other hand, the situations "Argue with someone", "Partner leaves you" and "Someone lies about you" are similar in terms of response. People who experience "Other anger", "Shame" and "Sorrow" in these situations are people who are characterised by their kindness, by the importance they place on other people’s judgements, and by their fear of rejection.

6. Conclusion

In this paper, it has been presented a generalization of the N-PLS model for the case where the response matrix is binary, called Bin-3-Way-PARAFAC-PLS. To the best of our knowledge, there is currently no model in the literature that deals with this specific situation. In addition, we have provided an algorithm to efficiently compute all parameters of the model.

For the visualization and interpretation of the results obtained from the proposed model, we have suggested the use of a biplot method, which facilitates the understanding of relationships between explanatory variables and response variables.

In order to validate the effectiveness and usefulness of the proposed model, we applied it to a dataset in the context of personality trait analysis. The results obtained demonstrate the ability of the model to provide significant insights in this type of analysis. It is worth noting that both the data used and the developed code are available.

While the BIN-3-Way-PARAFAC-PLS method is currently limited to three-way data arrays, future research could explore its generalization to N-way data arrays. Extending the approach to handle N-way structures would significantly broaden its applicability and versatility. Another promising line of future investigation is the use of the Tucker3 model instead of the PARAFAC model for the decomposition of the three-way data array, which could potentially improve the model’s flexibility and performance. Additionally, exploring the situation where the response binary data is a three-way data array or considering the possibility of having two multi-way binary arrays would further expand the model’s applicability.