Partial Least Squares Regression for Binary Data

1. Introduction

Due to the great technological advances nowadays, it is increasingly common to find datasets that contain a large number of variables that have been measured in a given sample. The set of variables can range from quantitative to binary, nominal or ordinal and, very often, a combination of several types.

The rapid growth of data science has led to an increasing demand for statistical techniques that can handle large and complex datasets.

In this context, traditional univariate statistical techniques, repeated for each variable, do not capture the structure of the data, specially the relationships among variables, so it will be necessary to use multivariate techniques that offer a more powerful approach. These methods focus on the study of relationships among variables, similarities and differences among individuals and the variables responsible for them, treating all available information simultaneously.

Probably the most used multivariate statistical technique is Principal Components Analysis (PCA). Pearson [1] proposed PCA over 120 years ago, which was later formalized by Hotelling [2], Hotelling [3]. PCA is a dimensionality reduction technique that can be used to summarize the information contained in a large dataset. PCA works by finding a set of orthogonal axes, called principal components, that capture the most variation in the data. It was not until a few years later, in 1971, that K.R. Gabriel [4] proposed a simultaneous graphical representation of the rows and columns of a data matrix, related to PCA, which he called biplot.

This type of methods will allow for visualizing the structure of the data, avoiding making use only of a list of p-values associated with the performed analysis. Biplots can also allow the visualization of more complex multivariate models as Canonical Variate Analysis or Multivariate Analysis of Variance (MANOVA) [5,6,7].

All of these techniques were developed for continuous variables, but when the measured variables are of a categorical type, either binary, nominal or ordinal, these representations are not adequate.

There are several dimension reduction methods for categorical data, for example Item Response Theory (see for example [8]) or Factor Analysis for categorical data (see for example [9]), but none of these techniques had an associated biplot representation.

We will focus here on biplot representations based on logistic responses that are closely related to the mentioned techniques. It is worth highlighting the work of Vicente-Villardon et al. [10], Demey et al. [11] who have developed a biplot for a single matrix of binary data, Hernández-Sánchez and Vicente-Villardón [12] for nominal data, Vicente-Villardón and Hernández Sánchez [13] for ordinal data and Vicente-Villardón and Hernández-Sánchez [14] who have developed a biplot for mixed data that contains several types of variables. These proposals are based on the generalized bivariate models proposed by Gabriel [15].

In many practical situations, we have two or more data matrices, with symmetric or asymmetric roles. In this work, we focus on the particular case in which we have two data matrices with asymmetric roles, a set of predictor variables and a set of response variables. Generally, the goal when we have this situation is to find models that allow predicting the responses from the predictors. In both data sets, it is possible to have variables of all the types mentioned above.

For two multivariate data matrices with one playing the role of responses, one of the most used multivariate statistical techniques is Partial Least Squares (PLS) [16,17]. PLS is a regression technique that can be used to predict a set of response variables from a set of predictor variables. PLS works by finding a set of linear combinations of the predictor variables that are most correlated with the response variables. The classical PLS models were developed mainly for continuous data and obtain dimension reductions in both, the responses and the predictors, in a way similar to PCA with the difference that the focus is on prediction rather than variance maximization. A biplot representation for PLS that can help in the interpretation of the results was proposed by Oyedele and Lubbe [18].

More recently, Vicente-Gonzalez and Vicente-Villardon [19] proposed a PLS method when all the responses are binary and predictors are continuous. This extended PLS generalized linear regression [20] for a single response variable (PLS1) to several responses (PLS2) in a non-trivial way because it imply a dimension reduction not present in the former work. This proposal goes with a biplot for both matrices, a continuous one for the predictors and a logistic one for the binary responses.

When the two matrices are binary none of the mentioned procedures are suitable. In this paper we propose a PLS (PLS2) procedure when both matrices, responses and predictors, have binary variables. We use the logistic biplots, mentioned before, to make a triplot representation that combines individuals and variables of the two initial matrices.

The database that initially motivated the development of this method is studying the relation among origin of 9 types of corn strains and the presence or absence of genes. The data is composed of two binary matrices, the first one containing indicators of the 9 classes and the second a set of binary genetic polymorphisms.

In Section 2 we will describe Partial Least Squares and the generalization to two sets of binary variables; Section 3 describes biplots and its construction for the PLS models and Section 4 presents an example to illustrate the proposed models. Conclusions are in Section 5.

2. Partial Least Squares

Let ${\mathbf{X}}_{I\times J}$ be a matrix of predictors with J variables and ${\mathbf{Y}}_{I\times K}$ be a matrix of responses with K variables, both measured in I individuals.

When we have only a few predictors (significantly fewer than the number of individuals) that are not strongly correlated, multivariate linear regression (MLR) is typically employed to predict the responses. However, when the predictors are highly collinear or when there is a large number of them, traditional linear regression models may be inadequate.A notable example occurs when the number of observations is significantly smaller than the number of predictors, resulting in the absence of parameter estimates for MLR. This situation commonly arises in genomic studies, where a set of gene expressions is used to predict tumor types.

Partial Least Squares (PLS) [16,17] is one of the most recognized and widely applied methods for modeling the relationships between two datasets, particularly when both sets of variables are continuous and the MLR is not suitable.

Traditionally, this technique is used for prediction, but it lacks a graphical representation to help interpret the relationships between the two sets. When the responses are continuous, Oyedele and Lubbe [18] introduced an associated biplot representation, although earlier, less formal versions exist, such as in [21]. More recently, PLS-Biplot has been applied to team effectiveness [22].

When the response matrix contains binary variables, Partial Least Squares Discriminant Analysis (PLS-DA) [23] is often employed, which essentially fits a PLS regression to a dummy or fictitious variable. Bastien et al. [20] introduced a PLS model for a single binary response, analogous to the PLS-1 model. More recently, Vicente-Gonzalez and Vicente-Villardon [19] proposed a non-trivial extension for multiple binary variables, incorporating dimension reduction for binary data and resulting in a PLS-2 model, which serves as an alternative to PLS-DA.

In regression, to address the discriminant problem, Logistic Regression (LR) can be used instead of Multivariate Linear Regression, allowing the binary nature of the responses to be captured. While this is often a valid approach, LR faces the same limitations as MLR regarding the number of individuals, variables, and collinearity. Additionally, the separation problem may also arise.

Currently, no PLS alternative exists for situations where both the predictor and response matrices are binary.

The aim of this section is to extend PLSR to cases where both the predictor and response matrices consist of binary variables, by using logistic functions instead of linear ones. This approach, referred to as Binary Partial Least Squares (BPLS), enables us to account for the binary nature of the responses and addresses the limitations of traditional PLS methods when dealing with binary data.

This generalization is based on the adaptation of the NIPALS algorithm for binary data proposed in previous works by Vicente-Gonzalez and Vicente-Villardon [19]. The adaptation consists of replacing the linear regression step in the NIPALS algorithm by a logistic regression step. This allows the procedure to capture the binary nature of the data. Previously we make an outline of the classical PLS-2 model.

2.1. Partial Least Squares Regression: NIPALS Algorithm

PLS seeks to identify linear combinations of continuous predictors $\mathbf{X}$ that optimally predict a set of continuous responses $\mathbf{Y}$. Here, we outline the classical NIPALS algorithm [24]. Both sets are usually column-centred and probably standardized. PLS operates by extracting latent variables that capture the most significant variance and covariance between predictors and responses. It uses iterative extraction of components to build a predictive model that can handle collinearity and high-dimensional data effectively.

The latent model for predictors $\mathbf{X}$, in S dimensions, is

$$\mathbf{X}=\mathbf{T}{\mathbf{P}}^T+\mathbf{E}=\widehat{\mathbf{X}}+\mathbf{E},$$

(1)

where ${\mathbf{P}}_{J\times S}$ are variable loadings, ${\mathbf{T}}_{I\times S}$ are scores for individuals and ${\mathbf{E}}_{I\times J}$ is a matrix of residuals.

In the same way $\mathbf{Y}$ the latent model

$$\mathbf{Y}=\mathbf{T}{\mathbf{Q}}^T+\mathbf{F}=\widehat{\mathbf{Y}}+\mathbf{F},$$

(2)

where ${\mathbf{Q}}_{K\times S}$ are loadings, ${\mathbf{T}}_{I\times S}$ are scores and ${\mathbf{F}}_{I\times K}$ are residuals.

Observe that both models share the matrix of scores that best predict the responses. The NIPALS algorithm to obtain the parameters can be described as follows:

Algorithm 1 (NIPALS Algorithm)

1: procedureNIPALS($\mathbf{X},\mathbf{Y},S$)   2:     for $s=1\to S$ do   3:         ${\mathbf{t}}_{\left(s\right)}\leftarrow {\mathbf{y}}_{\left(j\right)}$ for some j. Or any other choice.           ▹ Init: ${\mathbf{t}}_{\left(s\right)}$   4:         repeat  5:            ${\mathbf{p}}_{\left(s\right)}\leftarrow {\mathbf{X}}^T{\mathbf{t}}_{\left(s\right)}/\Vert {\mathbf{X}}^T{\mathbf{t}}_{\left(s\right)}\Vert$                  ▹ Update: ${\mathbf{p}}_{\left(s\right)}$   6:            ${\mathbf{t}}_{\left(s\right)}\leftarrow {\mathbf{Xp}}_{\left(s\right)}$                       ▹ Update: ${\mathbf{t}}_{\left(s\right)}$   7:            ${\mathbf{q}}_{\left(s\right)}\leftarrow {\mathbf{Y}}^T{\mathbf{t}}_{\left(s\right)}/\Vert {\mathbf{Y}}^T{\mathbf{t}}_{\left(s\right)}\Vert$                     ▹ Update: ${\mathbf{q}}_{\left(s\right)}$   8:            ${\mathbf{t}}_{\left(s\right)}\leftarrow {\mathbf{Yq}}_{\left(s\right)}$                      ▹ Update: ${\mathbf{t}}_{\left(s\right)}$    9:         until ${\mathbf{t}}_{\left(s\right)}$ does not change   10:         $\mathbf{X}\leftarrow \mathbf{X}-{\mathbf{t}}_{\left(s\right)}{\mathbf{p}}_{\left(s\right)}^T$                    ▹ Update: $\mathbf{X}$    11:         $\mathbf{Y}\leftarrow \mathbf{Y}-{\mathbf{t}}_{\left(s\right)}{\mathbf{q}}_{\left(s\right)}^T$                     ▹ Update: $\mathbf{Y}$  12:     end for   return $\mathbf{T}=\left[{\mathbf{t}}_{\left(1\right)},\cdots ,{\mathbf{t}}_{\left(S\right)}\right]$, $\mathbf{P}=\left[{\mathbf{p}}_{\left(1\right)},\cdots ,{\mathbf{p}}_{\left(S\right)}\right]$ and $\mathbf{Q}=\left[{\mathbf{q}}_{\left(1\right)},\cdots ,{\mathbf{q}}_{\left(S\right)}\right]$ 13: end procedure

Note that by removing the step where predictors are updated in the NIPALS algorithm, we would end up with the principal components for the responses. The procedure alternatively updates the scores in matrix $\mathbf{T}$ using the information of responses and predictors.

Here, ${\mathbf{P}}^T\mathbf{P}=\mathbf{I}$ and ${\mathbf{Q}}^T\mathbf{Q}=\mathbf{I}$, and $\mathbf{T}$ contains the scores from $\mathbf{X}$ that best explain the set of responses. We have the relationship:

$$\mathbf{T}=\mathbf{X}\mathbf{P}.$$

(3)

A potential issue with this and subsequent algorithms is that they may have multiple starting points, leading to different solutions. The algorithm generates sequences of decreasing residual sum of squares values and eventually converges to at least a local minimum.

If we need the regression coefficients to express the $\mathbf{Y}$ variables as functions of the predictors $\mathbf{X}$, we have

$$\widehat{\mathbf{Y}}=\mathbf{X}\mathbf{B}.$$

(4)

Considering that $\widehat{\mathbf{Y}}=\mathbf{T}{\mathbf{Q}}^T$, we have

$$\widehat{\mathbf{Y}}=\mathbf{T}{\mathbf{Q}}^T=\mathbf{X}\mathbf{P}{\mathbf{Q}}^T.$$

(5)

Thus, the regression coefficients in terms of the original variables are

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

(6)

We have to note that the updating steps in the Algorithm 1 minimize the sum os squares of the residuals that can also be understood as a lost function.

$$\sum _{i=1}^I\sum _{j=1}^J{(y_{ij}-{\widehat{y}}_{ij})}^2=\sum _{i=1}^I\sum _{j=1}^J{(y_{ij}-\sum _{s=1}^S{\mathbf{t}}_{is}^t{\mathbf{q}}_{js})}^2=tr\left({(\mathbf{Y}-\widehat{\mathbf{Y}})}^T(\mathbf{Y}-\widehat{\mathbf{Y}})\right)=tr\left({(\mathbf{Y}-\mathbf{T}{\mathbf{Q}}^T)}^T(\mathbf{Y}-\mathbf{T}{\mathbf{Q}}^T)\right)$$

(7)

for the responses and in a similar way for the predictors.

2.2. PLS for Binary Responses

When responses are binary, the updating steps of the NIPALS algorithm (steps 7 and 8) are inadequate because the implicit linear relationship it assumes is not adequate. This is similar to how linear regression is unsuitable for binary responses, with logistic regression being the more appropriate method. The difference in step 11 makes also no sense.

Recently, Vicente-Gonzalez and Vicente-Villardon [19] generalized the NIPALS algorithm to handle a set of binary responses by incorporating a logistic relationship with the latent variables.

The expected values for the responses are now probabilities, and we use the logit function as the link function. Let $\widehat{\mathbf{Y}}={\Pi }^y$, then we have:

$$\mathrm{logit}\left({\Pi }^y\right)=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T.$$

(8)

This equation generalizes equation (2), except that we need to include a vector ${\mathbf{q}}_0$ with intercepts for each variable, as binary variables cannot be centred in the same way as continuous. Each probability ${\pi }_{ik}^y$ can be expressed as

$${\pi }_{ik}^y={\displaystyle \frac{e^{\left(q_{k0}+{\sum }_{s=1}^S{t}_{is}{q}_{ks}\right)}}{1+e^{\left(q_{k0}+{\sum }_{s=1}^S{t}_{is}{q}_{ks}\right)}}}.$$

(9)

Or in matrix form, if we have

$${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T.$$

(10)

then

$${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y}),$$

(11)

where the operations apply to each element of the matrices.

The updating steps of the responses have to be changed for binary responses. Rather than minimizing the cost function in Equation (7) we minimize the cost defined as

$$L^y=\sum _{i=1}^I\sum _{k=1}^K\left[-y_{ik}log\left({\pi }_{ik}^y\right)-(1-y_{ik})log\left(1-{\pi }_{ik}^y\right)\right].$$

(12)

We interpret here the function as a cost to minimize rather than a likelihood to maximize. We look for the parameters $\mathbf{T}$, $\mathbf{Q}$ and ${\mathbf{q}}_0$ that minimize the cost function. There are not closed-form solutions for the optimization problem so an iterative algorithm in which, we obtain a sequence decreasing values of the lost function at each iteration, is used. We will use the gradient method in a recursive way calculating one component at each time fixing the previous one. The update for each parameter would be as follows:

$$\begin{array}{c}\hfill q_{j0}\Leftarrow q_{j0}-\alpha {\displaystyle \frac{\partial L^y}{\partial q_{j0}}}=q_{j0}-\alpha \sum _{i=1}^I({\pi }_{ij}^y-y_{ij});\end{array}$$

(13)

$$\begin{array}{c}\hfill t_{is}\Leftarrow t_{is}-\alpha {\displaystyle \frac{\partial L^y}{\partial t_{is}}}=t_{is}-\alpha \sum _{j=1}^J{q}_{js}({\pi }_{ij}^y-y_{ij});\end{array}$$

(14)

$$\begin{array}{c}\hfill q_{js}\Leftarrow q_{js}-\alpha {\displaystyle \frac{\partial L^y}{\partial q_{js}}}=q_{js}-\alpha \sum _{i=1}^I{t}_{is}({\pi }_{ij}^y-y_{ik}),\end{array}$$

(15)

for some $\alpha$. Or, in matrix form

$$\begin{array}{c}\hfill {\mathbf{q}}_0\Leftarrow {\mathbf{q}}_0-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{1}}_I;\end{array}$$

(16)

$$\begin{array}{c}\hfill {\mathbf{t}}_{\left(s\right)}\Leftarrow {\mathbf{t}}_{\left(s\right)}-\alpha ({\Pi }^y-\mathbf{Y}){\mathbf{q}}_{\left(s\right)};\end{array}$$

(17)

$$\begin{array}{c}\hfill {\mathbf{q}}_{\left(s\right)}\Leftarrow {\mathbf{q}}_{\left(s\right)}-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{t}}_{\left(s\right)},\end{array}$$

(18)

where ${\mathbf{t}}_{\left(s\right)}=(t_{1s},\cdots ,t_{Is})$ and ${\mathbf{q}}_{\left(s\right)}=(q_{1s},\cdots ,q_{Js})$. Then, the adaptation of the algorithm is straightforward. We present a recursive procedure that calculates one component at a time while keeping the previously computed components fixed.

Algorithm 2 (NIPALS-Binary responses)

1: procedureNIPALSBin($\mathbf{X},\mathbf{Y},S,\alpha$, $ϵ$)   2:     Set $\mathbf{T}=\left[\right]$, $\mathbf{P}=\left[\right]$ and $\mathbf{Q}=\left[\right]$ (empty matrices)   3:     ${\mathbf{q}}_0={\displaystyle \frac{1}{I}}{\mathbf{Y}}^T{\mathbf{1}}_I$ or any other good guess.                   ▹ Init: ${\mathbf{q}}_0$   4:     repeat   5:         ${\Pi }^y=e^{\mathbf{1}{\mathbf{q}}_{\mathbf{0}}^{\mathbf{T}}}÷(1+e^{\mathbf{1}{\mathbf{q}}_0^T})$   6:         ${\mathbf{q}}_0\Leftarrow {\mathbf{q}}_0-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{1}}_I$                     ▹ Update: ${\mathbf{q}}_0$   7:     until ${\mathbf{q}}_0$ does not change  8:     for $s=1\to S$ do    9:         Init: ${\mathbf{t}}_{\left(s\right)}\leftarrow random$ or any other good guess.    10:         Init: ${\mathbf{q}}_{\left(s\right)}\leftarrow random$ or any other good guess.   11:         repeat   12:            ${\mathbf{p}}_{\left(s\right)}\leftarrow {\mathbf{X}}^T{\mathbf{t}}_{\left(s\right)}/\Vert {\mathbf{X}}^T{\mathbf{t}}_{\left(s\right)}\Vert$                     ▹ Update: ${\mathbf{p}}_{\left(s\right)}$   13:            ${\mathbf{t}}_{\left(s\right)}\leftarrow {\mathbf{Xp}}_{\left(s\right)}$                         ▹ Update: ${\mathbf{t}}_{\left(s\right)}$   14:            $\mathbf{Q}=[\mathbf{Q};{\mathbf{q}}_{\left(s\right)}]$, $\mathbf{P}=[\mathbf{P};{\mathbf{p}}_{\left(s\right)}]$ and $\mathbf{T}=[\mathbf{T};{\mathbf{t}}_{\left(s\right)}]$ ▹ Add column ${\mathbf{q}}_{\left(s\right)}$ to $\mathbf{Q}$, ${\mathbf{p}}_{\left(s\right)}$ to $\mathbf{P}$ and ${\mathbf{t}}_{\left(s\right)}$ to $\mathbf{T}$   15:            ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   16:            ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   17:            $L_{old}^y\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{Y}\odot {\Pi }^y)+((1-\mathbf{Y})\odot (1-{\Pi }^y))){\mathbf{1}}_J$   18:            repeat   19:                ${\mathbf{q}}_{\left(s\right)}\Leftarrow {\mathbf{q}}_{\left(s\right)}-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{t}}_{\left(s\right)}$   20:                ${\mathbf{q}}_s\leftarrow {\mathbf{q}}_{\left(s\right)}$    ▹ Change the $s-th$ column ${\mathbf{q}}_s$ of $\mathbf{Q}$ for the new guess ${\mathbf{q}}_{\left(s\right)}$   21:                ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   22:                ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   23:            until ${\mathbf{q}}_{\left(s\right)}$ does not change   24:            repeat   25:                ${\mathbf{t}}_{\left(s\right)}\Leftarrow {\mathbf{t}}_{\left(s\right)}-\alpha ({\Pi }^y-\mathbf{Y}){\mathbf{q}}_{\left(s\right)}$   26:                ${\mathbf{t}}_s\leftarrow {\mathbf{t}}_{\left(s\right)}$     ▹ Change the $s-th$ column ${\mathbf{t}}_s$ of $\mathbf{T}$ for the new guess ${\mathbf{t}}_{\left(s\right)}$   27:                ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   28:                ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   29:            until ${\mathbf{t}}_{\left(s\right)}$ does not change   30:            $L_{new}^y\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{Y}\odot {\Pi }^y)+((1-\mathbf{Y})\odot (1-{\Pi }^y))){\mathbf{1}}_J$   31:         until $\mid L_{old}^y-L_{new}^y\mid <ϵ$   32:         $\mathbf{X}\leftarrow \mathbf{X}-{\mathbf{t}}_{\left(s\right)}{\mathbf{p}}_{\left(s\right)}^T$                      ▹ Update: $\mathbf{X}$   33:     end for   return ${\mathbf{q}}_0$, $\mathbf{T}=\left[{\mathbf{t}}_{\left(1\right)},\cdots ,{\mathbf{t}}_{\left(S\right)}\right]$, $\mathbf{P}=\left[{\mathbf{p}}_{\left(1\right)},\cdots ,{\mathbf{p}}_{\left(S\right)}\right]$ and $\mathbf{Q}=\left[{\mathbf{q}}_{\left(1\right)},\cdots ,{\mathbf{q}}_{\left(S\right)}\right]$   34: end procedure

In practice, the choice of $\alpha$ can be avoided by using a pre-programmed optimization routine. This algorithm is implemented in the package MultBiplotR [25] developed in the R language [26].

In the same way as for the continuous case, the regression coefficients in terms of the original variables are

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

(19)

Algorithm 2 is essentially the same as algorithm 1 except that we have made the necessary adaptations to cope with binary responses. Now, the relation between the manifest and latent variables is logistic rather than linear. Observe that is also a linear relation on the logits, i. e. a generalized linear model.

For a direct extension of this algorithm, see[27], implemented in the package BiplotML [28].

2.3. Partial Least Squares for Two Binary Matrices

Now both ${\mathbf{X}}_{I\times J}$ and ${\mathbf{Y}}_{I\times K}$ contain binary variables measured in the same I individuals. The goal is still predicting $\mathbf{Y}$ from $\mathbf{X}$ and describing the common structure of predictors and responses.

We denote the expected values of $\mathbf{X}$ as $E\left[\mathbf{X}\right]={\Pi }^x$ and the expected values of $\mathbf{Y}$ as $E\left[\mathbf{Y}\right]={\Pi }^y$. Then we have

$$\mathrm{logit}\left({\Pi }^{\mathrm{y}}\right)=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T;$$

(20)

and

$$\mathrm{logit}\left({\Pi }^{\mathrm{x}}\right)=\mathbf{1}{\mathbf{p}}_0^T+\mathbf{T}{\mathbf{P}}^T.$$

(21)

These Equations are generalizations of the decompositions of the matrices $\mathbf{X}$ and $\mathbf{Y}$ for the continuous case in Equations 1 and 2. The only difference is that we need to add the intersection vector for each variable (${\mathbf{p}}_0$ and ${\mathbf{q}}_0$ respectively). This is because, unlike the continuous case, the variables cannot be centered beforehand.

The relation among observed and latent variables, can be written as:

$${\pi }_{ik}^y={\displaystyle \frac{e^{\left(q_{k0}+{\sum }_{s=1}^S{t}_{is}{q}_{ks}\right)}}{1+e^{\left(q_{k0}+{\sum }_{s=1}^S{t}_{is}{q}_{ks}\right)}}};{\pi }_{ij}^x={\displaystyle \frac{e^{\left(p_{j0}+{\sum }_{s=1}^S{t}_{is}{p}_{js}\right)}}{1+e^{\left(p_{j0}+{\sum }_{s=1}^S{t}_{is}{p}_{js}\right)}}}.$$

(22)

That is, a logistic rather than linear relation.

Now, we have to design an algorithm to calculate the model parameters. The procedure alternatively updates the common scores in $\mathbf{T}$ using the information of responses and predictors. We use the gradient method and the lost functions as follows

$$L_y=\sum _{i=1}^I\sum _{k=1}^K\left[-y_{ik}log\left({\pi }_{ik}^y\right)-(1-y_{ik})log\left(1-{\pi }_{ik}^y\right)\right],$$

(23)

and

$$L_x=\sum _{i=1}^I\sum _{j=1}^J\left[-x_{ij}log\left({\pi }_{ij}^x\right)-(1-x_{ij})log\left(1-{\pi }_{ij}^x\right)\right].$$

(24)

We will update the scores using a gradient descent procedures as in Equations (16), (17) and (18) for the responses. The extension for predictors in $\mathbf{T}$ is straightforward.

$$\begin{array}{c}\hfill {\mathbf{p}}_0\Leftarrow {\mathbf{p}}_0-\alpha {({\Pi }^x-\mathbf{X})}^T{\mathbf{1}}_I\end{array}$$

(25)

$$\begin{array}{c}\hfill {\mathbf{t}}_{\left(s\right)}\Leftarrow {\mathbf{t}}_{\left(s\right)}-\alpha ({\Pi }^x-\mathbf{X}){\mathbf{p}}_{\left(s\right)}\end{array}$$

(26)

$$\begin{array}{c}\hfill {\mathbf{p}}_{\left(s\right)}\Leftarrow {\mathbf{p}}_{\left(s\right)}-\alpha {({\Pi }^x-\mathbf{X})}^T{\mathbf{t}}_{\left(s\right)},\end{array}$$

(27)

2.3.1. NIPALS Algorithm for Two Binary Matrices

Then, we will search for the parameters $\mathbf{T}$, $\mathbf{P}$ y ${\mathbf{p}}_0$ for the predictors and $\mathbf{T}$, $\mathbf{Q}$ y ${\mathbf{q}}_0$ for the responses. We will obtain the scores $\mathbf{T}$ for $\mathbf{X}$ that better predict $\mathbf{Y}$ with a logistic regression model. The scores are also related to predictors with a logistic response model.

The calculations can be organized in an alternating algorithm that calculates the parameters ${\mathbf{t}}_{\left(s\right)}=(t_{1s},\cdots ,t_{Is})$, ${\mathbf{p}}_{\left(s\right)}=(p_{1s},\cdots ,p_{Js})$ and ${\mathbf{q}}_{\left(s\right)}=(q_{1s},\cdots ,q_{Ks})$ for each dimension s, while fixing the parameters already obtained in the previous dimensions. First, we need to calculate the constant (dimension 0) ${\mathbf{p}}_{\left(0\right)}=(p_{10},\cdots ,p_{J0})$ and ${\mathbf{q}}_{\left(0\right)}=(p_{10},\cdots ,q_{K0})$ separately. This is so because data can not be centered before hand.

The algorithm is as follows.

Algorithm 3 (PLS for two binary matrices - BPLS)

1: procedureBPLS($\mathbf{X},\mathbf{Y},S,\alpha$, $ϵ$)   2:     Set $\mathbf{T}=\left[\right]$, $\mathbf{P}=\left[\right]$ and $\mathbf{Q}=\left[\right]$ (empty matrices)   3:     ${\mathbf{q}}_0={\displaystyle \frac{1}{I}}{\mathbf{Y}}^T{\mathbf{1}}_I$ and ${\mathbf{p}}_0={\displaystyle \frac{1}{I}}{\mathbf{X}}^T{\mathbf{1}}_I$ or any other good guesses.     ▹ Init: ${\mathbf{q}}_0$ and ${\mathbf{p}}_0$   4:     repeat   5:         ${\Pi }^y=e^{\mathbf{1}{\mathbf{q}}_{\mathbf{0}}^{\mathbf{T}}}÷(1+e^{\mathbf{1}{\mathbf{q}}_0^T})$   6:         ${\mathbf{q}}_0\Leftarrow {\mathbf{q}}_0-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{1}}_I$                ▹ Update: ${\mathbf{q}}_0$   7:     until ${\mathbf{q}}_0$ does not change   8:     repeat   9:         ${\Pi }^x=e^{\mathbf{1}{\mathbf{p}}_{\mathbf{0}}^{\mathbf{T}}}÷(1+e^{\mathbf{1}{\mathbf{p}}_0^T})$   10:         ${\mathbf{p}}_0\Leftarrow {\mathbf{p}}_0-\alpha {({\Pi }^x-\mathbf{X})}^T{\mathbf{1}}_I$                  ▹ Update: ${\mathbf{p}}_0$   11:     until ${\mathbf{p}}_0$ does not change   12:     for $s=1\to S$ do    13:         Init: ${\mathbf{t}}_{\left(s\right)}$, ${\mathbf{q}}_{\left(s\right)}$ and ${\mathbf{p}}_{\left(s\right)}$ with random values or any other good guesses.    14:         $\mathbf{Q}=[\mathbf{Q};{\mathbf{q}}_{\left(s\right)}]$, $\mathbf{P}=[\mathbf{P};{\mathbf{p}}_{\left(s\right)}]$ and $\mathbf{T}=[\mathbf{T};{\mathbf{t}}_{\left(s\right)}]$  ▹ Add column ${\mathbf{q}}_{\left(s\right)}$ to $\mathbf{Q}$, ${\mathbf{p}}_{\left(s\right)}$ to $\mathbf{P}$ and ${\mathbf{t}}_{\left(s\right)}$ to $\mathbf{T}$   15:         repeat   16:            ${\mathbf{M}}^x=\mathbf{1}{\mathbf{p}}_0^T+\mathbf{T}{\mathbf{P}}^T$   17:            ${\Pi }^x=e^{{\mathbf{M}}^x}÷(1+e^{{\mathbf{M}}^x})$   18:            repeat   19:                ${\mathbf{p}}_{\left(s\right)}\Leftarrow {\mathbf{p}}_{\left(s\right)}-\alpha {({\Pi }^x-\mathbf{X})}^T{\mathbf{t}}_{\left(s\right)}$   20:                ${\mathbf{p}}_s\leftarrow {\mathbf{p}}_{\left(s\right)}$    ▹ Change the $s-th$ column ${\mathbf{p}}_s$ of $\mathbf{P}$ for the new guess ${\mathbf{p}}_{\left(s\right)}$   21:                ${\mathbf{M}}^x=\mathbf{1}{\mathbf{p}}_0^T+\mathbf{T}{\mathbf{P}}^T$   22:                ${\Pi }^x=e^{{\mathbf{M}}^x}÷(1+e^{{\mathbf{M}}^x})$ 23:            until ${\mathbf{p}}_{\left(s\right)}$ does not change   24:            repeat   25:                ${\mathbf{t}}_{\left(s\right)}\Leftarrow {\mathbf{t}}_{\left(s\right)}-\alpha ({\Pi }^x-\mathbf{X}){\mathbf{p}}_{\left(s\right)}$   26:                ${\mathbf{t}}_s\leftarrow {\mathbf{t}}_{\left(s\right)}$    ▹ Change the $s-th$ column ${\mathbf{t}}_s$ of $\mathbf{T}$ for the new guess ${\mathbf{t}}_{\left(s\right)}$   27:                ${\mathbf{M}}^x=\mathbf{1}{\mathbf{p}}_0^T+\mathbf{T}{\mathbf{P}}^T$   28:                ${\Pi }^x=e^{{\mathbf{M}}^x}÷(1+e^{{\mathbf{M}}^x})$   29:            until ${\mathbf{t}}_{\left(s\right)}$ does not change   30:            ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   31:            ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   32:            $L_{old}^y\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{Y}\odot {\Pi }^y)+((1-\mathbf{Y})\odot (1-{\Pi }^y))){\mathbf{1}}_J$   33:            repeat   34:                ${\mathbf{q}}_{\left(s\right)}\Leftarrow {\mathbf{q}}_{\left(s\right)}-\alpha {({\Pi }^y-\mathbf{Y})}^T{\mathbf{t}}_{\left(s\right)}$   35:                ${\mathbf{q}}_s\leftarrow {\mathbf{q}}_{\left(s\right)}$    ▹ Change the $s-th$ column ${\mathbf{q}}_s$ of $\mathbf{Q}$ for the new guess ${\mathbf{q}}_{\left(s\right)}$   36:                ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   37:                ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   38:            until ${\mathbf{q}}_{\left(s\right)}$ does not change   39:            repeat   40:                ${\mathbf{t}}_{\left(s\right)}\Leftarrow {\mathbf{t}}_{\left(s\right)}-\alpha ({\Pi }^y-\mathbf{Y}){\mathbf{q}}_{\left(s\right)}$   41:                ${\mathbf{t}}_s\leftarrow {\mathbf{t}}_{\left(s\right)}$    ▹ Change the $s-th$ column ${\mathbf{t}}_s$ of $\mathbf{T}$ for the new guess ${\mathbf{t}}_{\left(s\right)}$   42:                ${\mathbf{M}}^y=\mathbf{1}{\mathbf{q}}_0^T+\mathbf{T}{\mathbf{Q}}^T$   43:                ${\Pi }^y=e^{{\mathbf{M}}^y}÷(1+e^{{\mathbf{M}}^y})$   44:            until ${\mathbf{t}}_{\left(s\right)}$ does not change   45:            $L_{new}^y\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{Y}\odot {\Pi }^y)+((1-\mathbf{Y})\odot (1-{\Pi }^y))){\mathbf{1}}_J$   46:         until $\mid L_{old}^y-L_{new}^y\mid <ϵ$   47:     end for   return ${\mathbf{p}}_0$, ${\mathbf{q}}_0$, $\mathbf{T}=\left[{\mathbf{t}}_{\left(1\right)},\cdots ,{\mathbf{t}}_{\left(S\right)}\right]$, $\mathbf{P}=\left[{\mathbf{p}}_{\left(1\right)},\cdots ,{\mathbf{p}}_{\left(S\right)}\right]$ and $\mathbf{Q}=\left[{\mathbf{q}}_{\left(1\right)},\cdots ,{\mathbf{q}}_{\left(S\right)}\right]$ 48: end procedure

In logistic regression, there is a potential problem known as separation. This occurs when the data are perfectly separable by a linear combination of the predictors. In this case, the maximum likelihood estimator does not exist and tends to infinity [29]. Even if the separation is not perfect (quasi-separation), the estimator can be highly unstable.

To avoid separation, the usual solution is to use a penalized likelihood [30]. This adds a penalty term to the likelihood function that encourages the coefficients to be smaller. In our case, we will use the Ridge penalty [31].

$$L^y=\sum _{i=1}^I\sum _{k=1}^K\left[-y_{ik}log\left({\pi }_{ik}^y\right)-(1-y_{ik})log\left(1-{\pi }_{ik}^y\right)\right]+\lambda \sum _{k=1}^K\sum _{s=1}^S{q}_{ks}^2.$$

Gradient adaptation for the penalized cost function is automatic:

$${\displaystyle \frac{\partial L^y}{\partial q_{ks}}}=\sum _{i=1}^I{t}_{is}({\pi }_{ik}^y-y_{ik})+2\lambda q_{ks}.$$

To address the separation problem, we can introduce a regularization parameter $\lambda$. This parameter penalizes large coefficients, making the solution less likely to be separated.

The BPLS algorithm (Algorithm 3) can be sensitive to the random start. To obtain the best possible solution, it is recommended to test several starting points. In practice, it has been observed that the solutions from different starting points produce similar predictions of the responses.

The BPLS algorithm can also suffer from the separation problem.

If we eliminate the part where ${\mathbf{p}}_0$ and $\mathbf{P}$ are updated in algorithms 3 and 2 we obtain a reduced dimension approximation of a single data matrix analogous to a multidimensional logistic two parameter model in Item Response Theory, see [32,33]. It is also similar to Principal Components Analysis for binary data as in [34,35,36,37], or to Factor analysis for binary data [9].

3. Biplots

A biplot [4] is a type of graph that simultaneously displays two sets of information: the scores of the observations and the loadings of the variables in a reduced-dimensional space. It is commonly used in techniques like Principal Component Analysis (PCA) and Partial Least Squares (PLS) regression, among many others, to visualize high-dimensional data in two or three dimensions.

Scores represent the projections of the original data points (observations) onto a reduced number of components or latent variables. Loadings indicate the contributions of the original variables to these components, showing how each variable relates to the reduced dimensions.

In a biplot, observations (samples) are often represented by points, while the loadings (variables) are shown as vectors (arrows). The angles and lengths of the arrows help interpret how variables are correlated and how they contribute to the variation in the data.

Formally, a biplot represents the approximate decomposition of a data matrix $\mathbf{X}$ into the product of two lower-rank matrices (typically two or three), as follows:

$$\mathbf{X}=E\left[\mathbf{X}\right]+\mathbf{E}=\widehat{\mathbf{X}}+\mathbf{E}=\mathbf{A}{\mathbf{B}}^T+\mathbf{E}$$

(28)

This decomposition approximates the original matrix as accurately as possible based on specific criteria, with $\mathbf{A}$ and $\mathbf{B}$ capturing the main structure and $\mathbf{E}$ representing the residuals or error. In fact, a biplot decomposes the expected approximate structure of the data matrix into the product of two other matrices. The matrix $\widehat{\mathbf{X}}$ captures the majority of the data’s characteristics based on a predefined criterion.

Matrices $\mathbf{A}$ and $\mathbf{B}$ can be represented as points into a two or three dimensional space in such a way that the inner product

$$x_{ij}\approx {\widehat{x}}_{ij}={\mathbf{a}}_i^T{\mathbf{b}}_j$$

(29)

approximates the element $x_{ij}$.

Biplots were initially proposed in connection to PCA [4] but have been used also accompanying other techniques as Multivariate Analysis of Variance (MANOVA) [5,6,38], STATIS-ACT [39], Correspondence Analysis [40] among many others. In the context of PLS for continuous data it is clear that Equations (1) and (2) define a biplot. The properties of those representations can be developed by Oyedele and Lubbe [18] and Vicente-Gonzalez and Vicente-Villardon [19].

Most biplots are based on the assumption of linear relationships between observed and latent variables and are primarily designed for continuous data. For binary data, however, these representations are not ideal, much like how linear regression struggles to accurately model relationships with a binary response.

Vicente-Villardon et al. [10] introduces a method where the relationship between the observed variables and the dimensions is modeled using a logistic approach. This method is linked to psychometric techniques, such as item response theory or latent trait analysis. The original paper proposes an algorithm based on alternating generalized regressions using the Newton-Raphson method to maximize the likelihood. However, the algorithm faces the same challenge as logistic regression when separation occurs. Specifically, for some binary variables, individuals displaying the presence of a characteristic are fully separated from those without it (by a hyperplane) in the final representation. This method is referred to as the Logistic Biplot.

Demey et al. [11] introduced a new algorithm that combines Principal Coordinates Analysis, Cluster Analysis, and Logistic Regression to create this type of biplot. The result is termed an "external logistic biplot" because it uses a two-step approach (Principal Coordinates for individuals and Logistic Regressions for variables) rather than simultaneously obtaining the row and column markers. This simplification helps to avoid the separation problem but comes at the cost of reduced goodness of fit.

More recently, Babativa-Márquez and Vicente-Villardón [41] proposed an algorithm that improves parameter fitting by employing the nonlinear conjugate gradient algorithm.

More recently, in the context of PLS methods, Vicente-Gonzalez and Vicente-Villardon [19] proposed a generalization of the NIPALS algorithm for binary responses, where the components are obtained iteratively using the gradient method. The paper also introduces a procedure to address the separation issue encountered in earlier algorithms for logistic biplots.

3.1. Logistic Biplots

Linear biplots decompose the expected values of the approximation into a product of two lower-rank matrices. In contrast, the logistic biplot decomposes the expected probabilities using the logit link function, similar to how generalized linear models operate.

For any binary data matrix ${\mathbf{X}}_{I\times J}$, if we call $E\left[\mathbf{X}\right]=\Pi$ the decomposition, using the logit link, is

$$\mathrm{logit}\left(\Pi \right)={\mathbf{1}}_I{\mathbf{b}}_0^T+{\mathbf{AB}}^T.$$

(30)

The constants ${\mathbf{b}}_0^T$ have to be included because data can not be centered beforehand as in the continuous case. With this decomposition, the inner product of the biplot markers (coordinates) is the logit of a expected probability, except for a constant,

$$logit\left({\pi }_{ij}\right)=b_{j0}+{\mathbf{a}}_i^T{\mathbf{b}}_j.$$

(31)

The logits are easily converted into probabilities

$${\pi }_{ij}={\displaystyle \frac{e^{\left(b_{j0}+{\mathbf{a}}_i^T{\mathbf{b}}_j\right)}}{1+e^{\left(b_{j0}+{\mathbf{a}}_i^T{\mathbf{b}}_j\right)}}}.$$

(32)

Then, by projecting a row marker ${\mathbf{a}}_i$ onto a column marker ${\mathbf{b}}_j$ we obtain, except for a constant, the expected logit and thus the expected probability for the entry $x_{ij}$ of $\mathbf{X}$. The constant $b_{j0}$ serves to determine the exact point where the logit is zero or any other value.

Due to the generalized nature of the model, the geometric interpretation closely resembles that of linear biplots. Computational procedures are analogous to those employed in previous cases, with the addition of the constant term. For example, in two dimensions, we can determine, on the direction of the vector ${\mathbf{b}}_j=\left(b_{j1},b_{j2}\right)$, what point predicts a given probability p. Let’s $\left(x,y\right)$ that point, satisfying the equation

$$y={\displaystyle \frac{b_{j2}}{b_{j1}}}x.$$

Prediction also verifies that

$$\mathrm{logit}\left(p\right)=b_{j0}+b_{j1}x+b_{j2}y.$$

Therefore, we obtain

$$x={\displaystyle \frac{\left(\mathrm{logit}\left(p\right)-b_{j0}\right)b_{j1}}{b_{j1}^2+b_{j2}^2}};y={\displaystyle \frac{\left(\mathrm{logit}\left(p\right)-b_{j0}\right)b_{j2}}{b_{j1}^2+b_{j2}^2}}.$$

(33)

The point in the direction of ${\mathbf{b}}_j$ that predicts $0.5$ ($logit\left(0,5\right)=0$), is

$$x={\displaystyle \frac{-b_{j0}{b}_{j1}}{b_{j1}^2+b_{j2}^2}};y={\displaystyle \frac{-b_{j0}{b}_{j2}}{b_{j1}^2+b_{j2}^2}}.$$

Using Equations (33), we can position markers for different probabilities along the direction of ${\mathbf{b}}_j$, creating a graded scale similar to a coordinate axis. The interpretation remains fundamentally similar to that of linear biplots; however, markers representing equidistant probabilities may not be spaced equally in the graph.

As an illustrative simple example with data taken from [42]. The data companies collect can include the expected, such as your name, date of birth, and email address, as well as more unusual details, such as your pets, hobbies, height, weight, and even personal preferences in the bedroom. They may also store your banking information and links to your social media accounts, along with any data you share there.

How companies use this data varies depending on their business, but it often leads to targeted advertising and optimizing website management.

The infographics taken from the site is shown in Figure 1. The picture shows information on your face, voice, and environment, that some internet companies collect.

The original data has been transformed into a binary matrix, where the rows represent internet companies, the columns represent the types of information collected, and each entry is marked as 1 if the company collects that information, and 0 if it does not. The data matrix is shown in Table 1. With this data we have performed a Logistic Biplot that summarizes the information in a two dimensional graph. We obtain a joint representation of rows and columns of the data matrix. Distances among companies will be interpreted as similarities, i. e., companies lying near on the biplot display have similar profiles in relation to the collected information. Angles between variables (kind of information) is interpreted as correlations. Acute small angles mean strong positive correlations, near straight angles, strong negative correlations and almost rights angles mean no correlation. Projecting the row markers onto column markers, we have the expected probability for each entry.

Figure 2 shows the graphical results as a typical logistic biplot with scales for probabilities 0.05, 0.25, 0.5, 0.75 and 0.95. We have used percents (5, 25, 50, 75, 95) rather than probabilities to simplify the plot. We could have used any other values for the probabilities.

The plot displays a two-dimensional scattergram, where points represent the rows (companies) of the data matrix, and vectors or directions correspond to the variables. Distances among points for companies are inversely related to their similarity, for example, points representing Facebook, Instagram and Tik-Tok coincide, indicating that the three companies have exactly the same profiles (the three collect all the available information). Companies with a single point have unique patterns, for example, Offerup, Uber or eBay, and groups with close points have several characteristics in common.

As previously mentioned, the angles between the variable directions provide an approximation of the correlations between them. For example, face recognition, environment recognition and languages have high positive correlations because the increasing probabilities point into the same direction, i. e. have small acute angles between them. Collecting the contacts is not correlated to this group because they form an straight angle with it. Contacts and image library are negatively correlations because the probabilities increase in opposite directions.

The angles are connected to tetrachoric correlation due to the close relationship between this representation and the factorization of the tetrachoric correlation matrix, as we will explain later.

Projecting companies into the direction of a variable, we have the expected probabilities given by the approximation. Figure 3 shows the projections onto the variable "Face Recognition". The point predicting 0.5 has been marked with a circle because is usually the cut point to predict presence or absence of the characteristic.

To simplify the plot we can place marks for just a few probabilities, for example, 0.5 and 0.75 indicating the point were expected logit is 0 (probability 0.5) and the direction of increasing probabilities. Figure 4 shows a simplified logistic biplot with probabilities 0.5 and 0.75.

We can use the expected probabilities from equation (32) to generate predicted binary values: $\left({\widehat{x}}_{ij}=1\right)$ if ${\pi }_{ij}\ge 0,5$ and 0 otherwise. This yields an expected binary matrix $\widehat{\mathbf{X}}=\left({\widehat{x}}_{ij}\right)$. In the plot, this indicates that for each variable, the plane is divided into two prediction regions by a line perpendicular to the variable’s arrow and passing through the point predicting 0.5. The region on the side of the arrow predicts the presence of the characteristic, while the opposite region predicts its absence. Figure 5 illustrates this statement.

A guide for the interpretation of the biplot can be found also in the supplementary material by Demey et al. [11].

A common overall goodness-of-fit measure is the percentage of correctly classified entries of the data matrix. To assess the quality of representation for individual rows (individuals) and columns (variables), we can calculate these percentages separately. We can also calculate the percentage of true positives (sensitivity) and true negatives (specificity).

To evaluate the quality of representation for each binary variable, we need a suitable measure that generalizes the measures of predictiveness for continuous data as in [43]. Equation (32) essentially defines a logistic regression model for each variable when the row coordinates are fixed, that allows to evaluate goodness of fit using pseudo-$R^2$ measures commonly associated to that model, for example, McFaden, Cox-Snell or Nagelkerke.

Table 2 contains fit measures for the Internet Companies Data. It contains fit measures for each column separately and for the complete table.

94% of all entries in the data table are correctly predicted by the previously described procedure, as shown in Figure 5. However, it is more insightful to examine the indices for each individual variable. This is because, in some cases, only a few variables are accurately represented. Such a situation can arise when dealing with a large set of variables, where only a small subset is actually relevant to the problem. For instance, in a genetic polymorphisms matrix, only certain variables may be of real interest.

In our example, all variables show a reasonable percentage of correct classifications, except for "Image Library," which correctly classifies only 80% of all values and 66.67% of the negatives. Additionally, we included the deviance for each variable in the model, using the latent dimensions as predictors, along with a p-value as an indicator of model fit. This measure, potentially adjusted for multiple comparisons, could be used to select the most significant variables, as demonstrated in [11]. We also show pseudo-$R^2$ measures. All have high values except for "Image Library" and "Voice Recog." that are worse predicted or represented on the graph.

By assuming that the observed categorical responses are discretized versions of continuous processes, we can leverage the close relationship between the presented model and the factorization of tetrachoric correlations to perform a factor analysis.

Consider a binary variable $X_j$. This variable is assumed to arise from an underlying continuous variable $X_j^{*}$ that follows a standard normal distribution. There is a threshold, ${\tau }_j$, which divide the continuous variable $X^{*}$ into two categories.

The relationship between $X_j$ and $X_j^{*}$ can be expressed as:

$$X_j=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{X}_j^{*}\le {\tau }_j\hfill \\ 1\hfill & \mathrm{if}{X}_j^{*}>{\tau }_j\hfill \end{array}\right.$$

The tetrachoric correlations are the correlations among the $X_j^{*},j=1,\cdots ,J$. Let $\mathbf{R}$ a $J\times J$ matrix containing the tetrachoric correlations among the J binary variables and let ${\tau }_j,j=1,\cdots ,J$ the thresholds.

We can factorize the matrix $\mathbf{R}$ as

$$\mathbf{R}\cong \mathrm{\Lambda }{\mathrm{\Lambda }}^T$$

(34)

where $\mathrm{\Lambda }$ contains the loadings of a linear factor model for the underlying variables.

It can be shown that there is a close relation between the factor model in Equation (34) and the model in Equation (32). You can find the details in [9]. The model in Equation (32) is actually equivalent to the multidimensional logistic two parameter model in Item Response Theory (IRT).

The factor model loadings $\mathrm{\Lambda }=({\lambda }_{js},j=1,\cdots ,J;s=1,\cdots ,S)$ and he thresholds ${\tau }_j,j=1,\cdots ,J)$ can be calculated from the parameters for the variables in our model as:

$${\tau }_j=b_{j0}{\left(1-\sum _{s=1}^S{b}_{js}^2\right)}^{-1/2}$$

(35)

$${\lambda }_{js}=b_{js}{\left(1-\sum _{s=1}^S{b}_{js}^2\right)}^{-1/2}$$

(36)

In this way we provide a classical factor interpretation to our model adding the loadings and communalities (sum of squares of all dimension loadings of each variable). Loadings measure the correlation among the dimensions or latent factors and the communalities the amount of variance of each variable explained by the factors. For our data, the loadings and communalities are in Table 3.

The two dimensional solution explains the 88% of the variance. All the comunalities are higher than 0.71 meaning that a good amount of the variance is explained by the dimensions. The communalities serve also to select important variables. Variables with low communalities may be explained by other dimension or have little importance for the description of the rows.

We can see also that Face Recognition, Environment Recognition and Language, are mostly related to the first dimension while, Product Recognition, Your contacts are related to the first dimension. Voice Recognition is related to both.

3.2. Triplot Representations for Binary Partial (Least Squares) Regression

Equations (1) and (2) define two biplots for the continuous case that share the scores for individuals in $\mathbf{T}$. Some of the properties of those biplots have been described by Oyedele and Lubbe [18].

In the same way, Equations (1) and (8) define two biplots that share scores in $\mathbf{T}$ but now the biplot for responses is a logistic biplot as described in subSection 3.1. The details of this representation can be found in [19]. For this case, Equation (19) defines a biplot for the regression coefficients in terms of the observed variables that can help in determining the more important predictors to explain the responses.

The Equations in (20) and (21) define two logistic biplots with the particularity that both share the row scores in $\mathbf{T}$. We can plot those together with the loadings for responses $\mathbf{Q}$ and predictors $\mathbf{P}$ to obtain a joint representation of individuals, responses and predictors, i. e. a triplot with three sets of markers. Both biplots and then, the triplot, are interpreted using the rules in subSection 3.1.

The measures of quality (Correct classifications, sensitivity, specificity, pseudo-$R^2$, factor loadings or communalities) for matrix $\mathbf{Y}$ can be used to identify which responses are correctly explained and what dimensions are needed for the explanation. Global qualities can be used to determine the number of dimensions needed for an adequate prediction, taking into account that some of the responses may not be well predicted by the model. In this case angles among directions for responses or projections of the scores on them are less important because are not optimized to explain the variability in data but to capture the relations with the predictors. The same occurs for the joint interpretation of rows and columns

4. Illustrative Example

Anderson et al. [44] estimates that approximately 30% of plant diseases are caused by phytopathogenic fungi. These fungi can have a significant impact on both ecology and agriculture, causing crop losses and environmental damage.

We focus on a specific fungus, Colletotrichum graminicola (or Glomerella graminícola). This fungus is a major pathogen of corn, causing a disease called anthracnose. Anthracnose produces spots on the leaves and stems of corn plants, which can lead to crop losses. There are some studies, such as that of O’Connell et al. [45], which collect the economic impact it had in the United States.

The data used in this study comes from a project carried out by the CIALE. The project collected samples of corn plants infected with Colletotrichum graminicola from 9 different countries: Argentina (AR), Brazil (BR), Canada (CA), Croatia (HR), Slovenia (SI), France (FR), Portugal (PT), Switzerland (CH) and the United States (US). The response matrix has the indicators (dummy variables) for the countries,

The predictors are a set of binary polymorphisms obtained from the RNA sequencing of the fungus Colletotrichum graminicola yielded a total of 13,183 variables. However, the sequencing method generates duplicate information, a variable for the presence and another for the absence, that are complementary; so we first removed those variables for the absence. This reduced the number of variables to 6,419. We then selected the 2,379 variables that were most significant based on their correlation with the response variable, i.e., the presence or absence of anthracnose.

The statistical analysis was carried out with the statistical software R [26]. We used the MultBiplotR package [25], which includes functions for performing the proposed method including biplot analysis.

The biplot representation of the responses is shown in Figure 6. Rather than the rows, the convex hulls containing all the points of each origin, have been represented.

The BPLS model achieved an overall classification accuracy of 90.94% considering the binary variables separately. The classification accuracy for each country, together with some additional measures of fit are shown in Table 4.

In our example, the strains from Brazil correctly classify all of their individuals. The rest have lower percentages of correct classification. Argentina, Croatia, France, Portugal and Slovenia have 0% Sensitivity, that is none of the presences is correctly classified. Switzerland has a sensitivity of 22.22 %. This kind of classification is more useful when we have several (possibly) correlated variables rather than a set of indicators for several classes.

In his particular case, in which we have the indicators of some disjoint classes, we can calculate another classification assigning each individual to the class with higher expected probability. The results of this classification are shown in Table 5.

With this classification only 49.51% of the strains are correctly classified into their origin. All the strains from Brazil, Canada, and France are correctly classified, but the strains from the USA are mistakenly grouped with those from Canada, while the remaining strains (Argentina, Croatia, France, Portugal, Slovenia, and Switzerland) are classified under France. This likely suggests that the characteristics of all European varieties (including Argentina) are quite similar, while Canada and the USA also share many similarities, with only Brazil showing distinct genetic characteristics. Experts have informed us that the Argentinian strains originated in Europe, which further supports this observation.

We have not included the predictors and its interpretation in this plot because of the big amount. We will show the interpretation in the next step of the analysis.

In conclusion we identify three groups with different genetic profiles according to the origin of the strains: Brazil, North America (USA, Canada) and Europe (Argentina, Croatia, France, Portugal, Slovenia, and Switzerland). We repeat the analysis recoding the whole set into those groups named: Brazil, North America and Europe (including Argentina).

The biplot for the new responses are in Figure 7. It is quite clear that all the groups are correctly classified.

Table 6 contains measures of fit for the responses. All measures are very high so no further comment is needed.

Table 7 contains the loadings and the comunalities that are very high for this data what means that the relations among the latent variables with the responses is very high.

The procedure is also helpful in identifying the variables that contribute to the prediction by examining the predictors most strongly associated with the latent variables used for prediction. We can add the predictors to the previous biplot to obtain a triplot. The resulting plot is shown in Figure 8.

Due to the high number of predictors the complete plot is too crowded and difficult to see and we may select just one fraction of the variables. We can use for example predictors that have high comunalities, or any other measure of fit.

We are interested in the predictors that point in the direction of each response and are highly related to the latent variables. For example, we can calculate the squared cosines between the directions of predictors and response to identify the most important polymorfisms for each group. Those cosines are easily calculated from the matrices in Equation (20). If the vector ${\mathbf{l}}_P={\left({\mathbf{P}}^2{\mathbf{1}}_K\right)}^{1/2}$ contains the lengths of the row vectors of $\mathbf{P}$ (directions of predictors on the plot) and ${\mathbf{l}}_Q={\left({\mathbf{Q}}^2{\mathbf{1}}_I\right)}^{1/2}$ contains the lengths of the row vectors of $\mathbf{Q}$ (directions of responses on the plot), where the exponents apply to each element of the matrices, then the cosines are in the matrix

$$\mathbf{C}={\mathbf{D}}_{{\mathbf{l}}_P}^{-1}\mathbf{P}{\mathbf{Q}}^T{\mathbf{D}}_{{\mathbf{l}}_Q}^{-1}$$

(37)

where ${\mathbf{D}}_{{\mathbf{l}}_P}^{-1}$ and ${\mathbf{D}}_{{\mathbf{l}}_Q}^{-1}$ are the diagonal matrices with the lengths in the diagonal. Each column of $\mathbf{C}$ contains the cosines of all the predictors with a single response, the near are cosines to 1 or -1, the closer are the direction of responses and predictors.

For each group we select the variables that have squared cosines higher than 0.9 and communalities higher than 0.85 to select also those related to the latent dimensions. The limits for the selection can be changed to select more or less predictors.

Figure 9 shows the variables most related to the direction that predicts Brazil and Table 8 the goodness of fit.

Figure 10 shows the variables most related to the direction that predicts Europe.

The number of variables highly associated with European tribes is significantly lower than those observed for Brazilian or American tribes. Nevertheless, the goodness of fit, as shown in Table 9, remains robust. Interestingly, all genes are directly related to the response variable, except for gene "409898", which shows an inverse relationship. Notably, this gene has a higher goodness of fit compared to the rest, highlighting its distinct contribution to the model. However, it is worth noting that the classification performance for the European strains does not reach the near-perfect results observed for Brazil. In particular, the percentage of cases correctly classified, as well as the sensitivity and specificity, do not reach 100%. This indicates a slightly lower overall predictive accuracy when analyzing the European dataset compared to the Brazilian one.

Figure 11 shows the variables most related to the direction that predicts North America.

The number of genes associated with the North American strain is significantly higher than those associated with the European and Brazilian strains. Only 10 genes are directly associated with the North American strain, while the rest have an inverse relationship with the response variable. The goodness of fit, as shown in Table 10, is high, indicating robust model performance. Nevertheless, the percentage of correctly classified cases does not reach 100%. Interestingly, in certain cases, either sensitivity or specificity achieves perfect accuracy, demonstrating the ability of the method to capture certain aspects of the data with exceptional precision.

5. Conclusions

• Partial Least Squares Regression (PLSR) has been assessed as a robust alternative to Multivariate Linear Regression for exploring relationships between two datasets, particularly in situations where predictor variables are either highly numerous or exhibit significant collinearity. This evaluation highlights the versatility and effectiveness of PLSR in addressing complex modeling challenges.
• We propose an extension of the PLS technique applicable to scenarios involving binary predictors and responses. This novel approach, termed Binary Partial Least Squares (BPLS), is based on a generalization of the NIPALS algorithm, adapted to accommodate binary data structures, thereby broadening the applicability of PLS techniques.
• Custom functions implementing the BPLS method for binary data have been developed in the R programming environment. These functions are included in the MultBiplotR package, which has been prepared for submission to the CRAN repository to ensure wider broader accessibility and usability for researchers and practitioners.
• The usefulness and performance of the proposed BPLS technique was validated through a real-world case study, focusing on the classification of Colletotrichum graminicola strains. The results showed promising classification accuracy, highlighting the potential of BPLS as a powerful tool for binary data analysis.