Logistic Biplots for Ordinal Variables Based on Alternated Gradient Descent on the Cumulative Probabilities, with an Application to Survey Data

1. Introduction

The Biplot method [9,19], allows for simultaneous representation of the I rows and J columns of a data matrix $\mathbf{X}$, in reduced dimensions, where rows usually correspond to individuals, objects or samples and columns to a set of variables measured on them. Classical Biplot methods are a graphical representation of a Principal Components Analysis (PCA) or Factor Analysis (FA) that are used to obtain linear combinations that successively maximize the extracted variability. From another point of view, Classical Biplots can be obtained from alternated regressions and calibrations in [14]. This approach is essentially an alternated least squares algorithm.

In the context multivariate analysis, biplots are not that popular, although principal components analysis still are. The reason for this is that biplots are probably not well understood and its properties not well explained and also because some lack of connection between the statistical and applied frameworks. We think that every time a principal components analysis is used to analyse data, a biplot could improve significantly such analysis.

For data with distributions from the exponential family, [11] describes bilinear regression as a method to estimate biplot parameters, but the procedure have never been implemented and the geometrical properties of the resulting representations have never been studied in detail. [6] proposes Principal Components Analysis for binary data based on an alternate procedure in which each iteration is performed using iterative majorization and [25] extends the procedure for sparse data matrices, none of those describe the associated biplot. [33] propose a biplot based on logistic responses called “Logistic Biplot” that is linear, and study the geometry of this kind of biplots and uses a estimation procedure that is slightly different from Gabriel’s method. A heuristic version of the procedure for large data matrices in which scores for individuals are calculated with an external procedure as Principal Coordinates Analysis is described in [8]. Method is called “External Logistic Biplot”. Binary Logistic Biplots have been successfully applied to different data sets, see for example, [5,8,16,31] or [17]. More recently, [30] have proposed some additional algorithms to calculate the parameters of the logistic principal components and [1] an algorithm for the biplot.

For nominal data, [22] proposed a biplot representation based on convex prediction regions for each category of a nominal variable. An alternated algorithm is used for the parameter estimation. In Section 2, biplots for continuous, binary and, in lesser extent, nominal variables are described.

Linear, binary or nominal logistic biplots are not adequate when data are ordinal and techniques as CATPCA or IRT for ordinal items can be used instead. The first represents rows and columns of the data matrix in a biplot and the second does not usually have a graphical representation. In Section 3 we extend the logistic biplot to ordinal data. A draft of this extension was already publish in a repository [35]. The resulting method is termed “ordinal logistic biplot” (OLB). Row scores are computed to have ordinal logistic responses along the dimensions and column parameters produce logistic response surfaces that, projected onto the space spanned by the row scores, define a linear biplot. A proportional odds model is used, obtaining a multidimensional model similar to a graded response model in the IRT framework. We study the geometry of such a representation and construct computational algorithms for estimating the model parameters and calculating the prediction directions (or axes) used for visualization on the biplot. Ordinal logistic biplots extend both CATPCA and IRT in the sense that gives a graphical representation for IRT similar in some way to the biplot for CATPCA. The algorithm is based on an alternated gradient descent method.

In Section 4 the main results are applied to the study of job satisfaction of doctorate (PhD) holders in Spain. Holders of doctorate degrees or other research qualifications are crucial to the creation, commercialization and dissemination of knowledge and to innovation. The proposed methods are used to extract useful information from the Spanish data from the international ’Survey on the careers of doctorate holders (CDH)’, jointly carried out by Eurostat, the Organization for Economic Co-operation and Development (OECD) and UNESCO’s Institute for Statistics (UIS).

2. Biplot for Continuous, Binary or Nominal Data

In this section we describe the biplot for continuous, binary and nominal data, being the first two treated in a greater extent because of the closer relation to the proposal of this paper.

2.1. Linear biplot for continuous data

Let ${\mathbf{X}}_{I\times J}$ be a data matrix containing the measures of J variables on I individuals. The S-dimensional principal components (biplot) model is formulated as:

$$x_{ij}=b_{j0}+\sum _{s=1}^S{b}_{js}{a}_{is}+e_{ij}=b_{j0}+{\mathbf{a}}_i^T{\mathbf{b}}_j+e_{ij},$$

(1)

or in matrix form

$$\mathbf{X}={\mathbf{1}}_I{\mathbf{b}}_0^T+\mathbf{A}{\mathbf{B}}^T+\mathbf{E}$$

(2)

where ${\mathbf{b}}_0^T$ is a vector of constants. ${\mathbf{b}}_0^T$ can be estimated beforehand and contains the column means (${\mathbf{b}}_0=\overline{\mathbf{x}}$), the scores $\mathbf{A}$ and the loadings $\mathbf{B}$ are matrices of rank S with I and J rows respectively, and $\mathbf{E}$ is a $I\times J$ matrix of errors or residuals. Then

$$E\left[\mathbf{X}\right]={\mathbf{1}}_I{\mathbf{b}}_0^T+\mathbf{A}{\mathbf{B}}^T,$$

(3)

In practice we avoid the constants by centring the data matrix. The reduced rank approximation of the centred data matrix is

$$\tilde{\mathbf{X}}=E\left[\mathbf{X}-{\mathbf{1}}_I{\mathbf{b}}_0^T\right]=\mathbf{A}{\mathbf{B}}^T$$

(4)

is usually obtained from its singular value decomposition (SVD),

$$\tilde{\mathbf{X}}=\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^T,$$

(5)

where $\mathbf{U}$ contains the eigenvectors of $\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T$, $\mathbf{V}$ the eigenvectors of ${\tilde{\mathbf{X}}}^T\tilde{\mathbf{X}}$ and $\mathbf{\Lambda }$ is a diagonal matrix containing the non-zero eigenvalues of both matrices that are the same. The SVD is closely related to the principal components.

We can choose $\mathbf{A}$ and $\mathbf{B}$ in different ways, using the first S columns of the matrices in the SVD, for example,

$$\mathbf{A}=\mathbf{U}\mathbf{\Lambda }\mathbf{B}=\mathbf{V}.$$

(6)

$\mathbf{A}$ contains the scores on the principal components and are used for visualization on a scattergram, $\mathbf{B}$ contains the coefficients of the principal components. This is called a JK-Biplot or RMP-Biplot by [9].

Another possible choice is

$$\mathbf{A}=\mathbf{U}\mathbf{B}=\mathbf{V}\mathbf{\Lambda }.$$

(7)

$\mathbf{B}$ contains the factor loadings and $\mathbf{A}$ the standardized factor scores. If the columns of $\tilde{\mathbf{X}}$ are standardized, $\mathbf{B}$a contains the correlations of the observed variables and the principal components. This is called a GH-Biplot or CMP-Biplot by [9].

Normally the scores $\mathbf{A}$ in equation 6 are used for visualization od the individuals and the loadings or correlations $\mathbf{B}$ in 7 to interpret the indirect relation to the components in a process called reification.

The joint representation is called a biplot because simultaneously plots the individuals and variables using the rows of $\mathbf{A}={({\mathbf{a}}_1,\dots ,{\mathbf{a}}_I)}^{\prime }$ and $\mathbf{B}={({\mathbf{b}}_1,\dots ,{\mathbf{b}}_J)}^{\prime }$ as markers, in such a way that the inner product ${\mathbf{a}}_i^T{\mathbf{b}}_j$ approximates the element ${\tilde{x}}_{ij}$ as close as possible.

If we consider the row markers $\mathbf{A}$ as fixed and the data matrix previously centred, the column markers can be computed by regression trough the origin:

$${\mathbf{B}}^{\prime }={\left({\mathbf{A}}^{\prime }\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\prime }(\mathbf{X}-{\mathbf{1}}_I{\overline{\mathbf{x}}}^{\prime }).$$

(8)

In the same way, fixing $\mathbf{B}$, $\mathbf{A}$ can be obtained as:

$${\mathbf{A}}^{\prime }={\left({\mathbf{B}}^{\prime }\mathbf{B}\right)}^{-1}{\mathbf{B}}^{\prime }{(\mathbf{X}-{\mathbf{1}}_I{\overline{\mathbf{x}}}^{\prime })}^{\prime }.$$

(9)

Alternating the steps (8) and (9) the product converges to the same subspace generated by the SVD of the centred data matrix. Regression step in equation (8) adjusts a separate linear regression for each column (variable) and interpolation step in equation (9), interpolates an individual using the column markers as reference. The procedure is in some way a kind of EM-algorithm in which the regression step is the maximization part and the interpolation step is the expectation part. An extension for frequency matrices can be found in [13]. In summary, the expected values on the original data matrix are obtained on the biplot using a simple scalar product, that is, projecting the point ${\mathbf{a}}_i$ onto the direction defined by ${\mathbf{b}}_j$. This is why row markers are usually represented as points and column markers as vectors (also called biplot axis by [19]).

The biplot axis can be completed with scales to predict individual values of the data matrix. To find the point on the biplot direction, that predicts a fixed value $\mu$ of the observed variable when an individual point is projected, we look for the point $(x,y)$ lies on the biplot axis, i. e. that verifies

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

and

$$\mu =b_{j0}+b_{j1}x+b_{j2}y$$

Solving for x and y, we obtain

$$x=(\mu -b_{j0}){\displaystyle \frac{b_{j1}}{b_{j1}^2+b_{j2}^2}}$$

and

$$y=(\mu -b_{j0}){\displaystyle \frac{b_{j2}}{b_{j1}^2+b_{j2}^2}}$$

or

$$(x,y)=(\mu -b_{j0}){\displaystyle \frac{{\mathbf{b}}_j}{{\mathbf{b}}_j^T{\mathbf{b}}_j}}$$

Therefore, the unit marker for the j-th variable is computed by dividing the coordinates of its corresponding marker by its squared length, several points for specific values of $\mu$ can be labelled to obtain a reference scale. If data are centred then $b_{j0}=0$ and labels can be calculated by adding the average to the value of $\mu$ ($\mu +{\overline{x}}_j$). The resulting representation will be like the one shown in Figure 1. The data for the representation was used to illustrate a biplot by [10].

Data contains the consumption of different types of proteins at different European countries. The aim is to detect groups of countries with similar behaviours and the variables responsible for them.

The biplot allows for the direct interpretation of the relations among rows and columns rather than the indirect relation trough the principal components.

If $\widehat{\mathbf{X}}$ the expected values of the biplot in reduced dimension $\widehat{\mathbf{X}}\approx {\mathbf{A}\mathbf{B}}^T$, the global goodness of fit is the amount of variability accounted by the prediction, that is,

$${\rho }^2=tr\left({\widehat{\mathbf{X}}}^T\widehat{\mathbf{X}}\right)/tr\left({\mathbf{X}}^T\mathbf{X}\right)$$

(10)

Even for matrices in which we obtain a good global fit, there may be some columns whose information is not accounted in the plot. Goodness of fit for each column is

$${\rho }_j^2=diag\left({\widehat{\mathbf{X}}}^T\widehat{\mathbf{X}}\right)÷diag\left({\mathbf{X}}^T\mathbf{X}\right)$$

(11)

where ÷ means the element by element operation. The quantity ${\rho }_j^2$ is also the R-Squared of the regression of each column of $\mathbf{X}$ on $\mathbf{A}$ as in equation 8. We call that quantity quality of representation of the variable. [18] call this quantity predictiveness of the column. The measures are used to identify which variables are most related to the representation or the variables whose information is preserved on the biplot. If two dimensions are not enough to account for the desired information, additional components can be used and qualities can be calculated adapting equations 10 and 11 to the desired dimensions.

2.2. Logistic Biplot for Binary data

2.2.1. Formulation and Geometry

Let $p_{ij}(i=1,\dots ,I;j=1,\dots ,J)$ the observed probability, either 0 or 1, that an individual i have a characteristic j, resulting in a binary data matrix ${\mathbf{P}}_{I\times J}$. The S-dimensional logistic principal components (biplot) is

$${\pi }_{ij}={\displaystyle \frac{e^{b_{j0}+{\sum }_k{b}_{jk}{a}_{ik}}}{1+e^{b_{j0}+{\sum }_k{b}_{jk}{a}_{ik}}}}={\displaystyle \frac{1}{1+e^{-(b_{j0}+{\sum }_k{b}_{jk}{a}_{ik})}}}.$$

(12)

Or in $logit$ scale,

$$logit\left({\pi }_{ij}\right)=log\left({\displaystyle \frac{{\pi }_{ij}}{1-{\pi }_{ij}}}\right)=b_{j0}+\sum _{s=1}^S{b}_{js}{a}_{is}=b_{j0}+{\mathbf{a}}_i^T{\mathbf{b}}_j,$$

(13)

where ${\pi }_{ij}=E\left(p_{ij}\right)$ are the expected probabilities that the character j be present at individual i and $a_{is}$ and $b_{js}$, $(i=1,\dots ,I;j=1,\dots ,J;s=1,...,S)$, are the model parameters used as row and column markers respectively. The model is a generalized (bi)linear model having the $logit$ as a link function. In matrix form:

$$logit\left(\mathbf{\Pi }\right)={\mathbf{1}}_I{{\mathbf{b}}^{\prime }}_0+{\mathbf{AB}}^{\prime },$$

(14)

where $\mathbf{\Pi }$ is the matrix of expected probabilities, ${\mathbf{1}}_I$ is a vector of ones and ${\mathbf{b}}_0=(b_{10},\dots ,b_{J0})$ is the vector containing intercepts that have been added because it is not possible to center the data matrix in the same way as in linear biplots. The points predicting different probabilities are on parallel straight lines on the biplot; this means that predictions on the logistic biplot are made in the same way as on the linear biplots, i. e., projecting a row marker ${\mathbf{a}}_i=(a_{i1},a_{i2})$ onto a column marker ${\mathbf{b}}_j=(b_{j1},b_{j2})$. (See [33] or [8]). The calculations for obtaining the scale markers are simple. To find the marker for a fixed probability $\pi$, we look for the point $(x,y)$ that predicts $\pi$ and is on the ${\mathbf{b}}_j$ direction, i. e., on the line joining the points (0, 0) and $(b_{j1},b_{j2})$, that is

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

and

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

We obtain

$$x={\displaystyle \frac{(logit\left(\pi \right)-b_{j0})b_{j1}}{b_{j1}^2+b_{j2}^2}}$$

and

$$y={\displaystyle \frac{(logit\left(\pi \right)-b_{j0})b_{j2}}{b_{j1}^2+b_{j2}^2}}$$

Several points for specific values of $\pi$ can be labelled to obtain a reference scale. From a practical point of view the most interesting value is $\pi =0.5$ because the line passing trough that point and perpendicular to the direction, divides the representation into two regions, one that predicts presence and other absence. Plotting that point and an arrow pointing to the direction of increasing probabilities should be enough for most practical applications.

A typical representation of a binary logistic biplot is given in Figure 2.

2.2.2. Parameter Estimation

The model in (12) is also a latent trait or item response theory model, in that ordination axes are considered as latent variables that explain the association between the observed variables. In this framework we suppose that individuals respond independently to variables, and that the variables are independent for given values of the latent traits. With these assumptions the likelihood function is:

$$\mathrm{Prob}(p_{ij}/({\mathbf{b}}_0,\mathbf{A},\mathbf{B}))=\prod _{i=1}^I\prod _{j=1}^J{\pi }_{ij}^{p_{ij}}{(1-{\pi }_{ij})}^{1-p_{ij}}.$$

(15)

Taking the logarithm of the likelihood function yields:

$$L=log(Prob(p_{ij}/({\mathbf{b}}_0,\mathbf{A},\mathbf{B}))=\sum _{i=1}^I\sum _{j=1}^J\left[p_{ij}log\left({\pi }_{ij}\right)+(1-p_{ij})log(1-{\pi }_{ij})\right].$$

(16)

For $\mathbf{A}$ fixed, (16) can be separated into J parts, one for each variable:

$$L=\sum _{J=1}^J{L}_j=\sum _{J=1}^J\left(\sum _{i=1}^I\left[p_{ij}log\left({\pi }_{ij}\right)+(1-p_{ij})log(1-{\pi }_{ij})\right]\right).$$

(17)

Maximizing each $L_j$ is equivalent to performing a standard logistic regression using the j-$th$ column of $\mathbf{X}$ as a response and the columns of $\mathbf{A}$ as regressors.

In the same way the probability function can be separated into several parts, one for each row of the data matrix, $L={\sum }_{i=1}^I{L}_i$. The details of the procedure to calculate the row or individual scores can be found in [33]. The procedure may be affected by the separation problem, that is, if in the reduced dimension subspace there is an hyperplane that separates the presences and the absences for any of the observed variables, the Maximum Likelihood method does not converge.

Other useful methods for the estimation of the parameters as Marginal Maximum Likelihood can be found in [2] in the context of the Item Response Theory (IRT).

[32] have used gradient descent methods, the cost function is the negative maximum likelihood

$$L=\sum _{i=1}^I\sum _{k=1}^K\left[-p_{ij}log\left({\pi }_{ij}\right)-(1-p_{ij})log\left(1-{\pi }_{ij}\right)\right],$$

(18)

We interpret here the function as a cost to minimize rather than a likelihood to maximize. We look for the parameters $\mathbf{A}$, $\mathbf{B}$ and ${\mathbf{b}}_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 b_{j0}\Leftarrow b_{j0}-\alpha {\displaystyle \frac{\partial L}{\partial b_{j0}}}=b_{j0}-\alpha \sum _{i=1}^I({\pi }_{ij}-p_{ij})\end{array}$$

(19)

$$\begin{array}{c}\hfill a_{is}\Leftarrow a_{is}-\alpha {\displaystyle \frac{\partial L}{\partial a_{is}}}=a_{is}-\alpha \sum _{j=1}^J{b}_{js}({\pi }_{ij}-p_{ij})\end{array}$$

(20)

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

(21)

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

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

(22)

$$\begin{array}{c}\hfill {\mathbf{a}}_{\left(s\right)}\Leftarrow {\mathbf{a}}_{\left(s\right)}-\alpha (\mathbf{\Pi }-\mathbf{P}){\mathbf{b}}_{\left(s\right)}\end{array}$$

(23)

$$\begin{array}{c}\hfill {\mathbf{b}}_{\left(s\right)}\Leftarrow {\mathbf{b}}_{\left(s\right)}-\alpha {(\mathbf{\Pi }-\mathbf{P})}^T{\mathbf{a}}_{\left(s\right)},\end{array}$$

(24)

where ${\mathbf{a}}_{\left(s\right)}=(a_{1s},\dots ,a_{Is})$ and ${\mathbf{b}}_{\left(s\right)}=(b_{1s},\dots ,b_{Js})$. For large data matrices it may be convenient to summarize the data into a matrix $\mathbf{P}$ containing unique patterns in the data and a vector with the frequencies of each pattern $\mathbf{f}$. The equations of the gradient descent method can be written now as:

$$\begin{array}{c}\hfill {\mathbf{b}}_0\Leftarrow {\mathbf{b}}_0-\alpha {(\mathbf{\Pi }-\mathbf{P})}^T\mathbf{f},\end{array}$$

(25)

$$\begin{array}{c}\hfill {\mathbf{a}}_{\left(s\right)}\Leftarrow {\mathbf{a}}_{\left(s\right)}-\alpha (\mathbf{\Pi }-\mathbf{P}){\mathbf{b}}_{\left(s\right)},\end{array}$$

(26)

$$\begin{array}{c}\hfill {\mathbf{b}}_{\left(s\right)}\Leftarrow {\mathbf{b}}_{\left(s\right)}-\alpha {(\mathbf{\Pi }-\mathbf{P})}^T(\mathbf{f}\odot {\mathbf{a}}_{\left(s\right)}),\end{array}$$

(27)

where ⊙ is the Hadamard product.

We can organize the calculations in an alternated algorithm that alternatively calculates parameters for rows and columns for each dimension s fixing the parameters already obtained for each previous dimensions. Before that, the constants ${\mathbf{b}}_0=(b_{10},\dots ,b_{J0})$ have to be fitted separately. This is equivalent to subtracting means to obtain the centred data matrix in the continuous biplot. The procedure would be as follows:

Algorithm 1 Algorithm to calculate the Components for Binary Data

1: procedureP-Binary-Components($\mathbf{P},S$) 2:     Choose $\alpha$ 3:     Init: ${\mathbf{b}}_0=random$ 4:     repeat 5:         Update: ${\mathbf{b}}_0$ with Equations 22 or 25 6:         Update : $\mathbf{\Pi }=\left({\pi }_{ij}\right)$ with Equation 12 7:     until ${\mathbf{b}}_0$ do not change 8:     for $s=1\to S$ do 9:         Init: ${\mathbf{a}}_{\left(s\right)}\leftarrow random$ (or any other choice) 10:         repeat 11:            repeat 12:                Update: ${\mathbf{b}}_{\left(s\right)}$ with Equation 24 or 27 13:                Update : $\mathbf{\Pi }=\left({\pi }_{ij}\right)$ with Equation 12 14:            until ${\mathbf{b}}_{\left(s\right)}$ do not change 15:            repeat 16:                Update: ${\mathbf{a}}_{\left(s\right)}$ with Equation 23 or 26 17:                Update : $\mathbf{\Pi }=\left({\pi }_{ij}\right)$ with Equation 12 18:            until ${\mathbf{a}}_{\left(s\right)}$ do not change 19:            $L\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{P}\odot \mathbf{\Pi })+((1-\mathbf{P})\odot (1-\mathbf{\Pi }))){\mathbf{1}}_J$ 20:         until L does not change 21:     end forreturn ${\mathbf{b}}_0$, $\mathbf{B}=\left[{\mathbf{b}}_{\left(1\right)},\dots ,{\mathbf{b}}_{\left(S\right)}\right]$, $\mathbf{A}=\left[{\mathbf{a}}_{\left(1\right)},\dots ,{\mathbf{a}}_{\left(S\right)}\right]$ 22: 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 ([34]) developed in the R language ([27]).

An direct extension of this algorithm can be found in [1].

2.3. Logistic Biplot for Nominal data

Let ${\mathbf{X}}_{I\times J}$ be a data matrix containing the values of J nominal variables, each with $K_j$$(j=1,\dots ,J)$ categories, for I individuals, and let ${\mathbf{P}}_{I\times L}$ be the corresponding indicator matrix with $L={\sum }_j{K}_j$ columns. The last (or the first) category of each variable will be used as a baseline. Let ${\pi }_{ij\left(k\right)}$ denote the expected probability that the category k of variable j be present at individual i. A multinomial logistic latent trait model with S latent traits, states that the probabilities are obtained as:

$${\pi }_{ij\left(k\right)}={\displaystyle \frac{e^{b_{j\left(k\right)0}+\sum _{s=1}^S{b}_{j\left(k\right)s}{a}_{is}}}{\sum _{l=1}^{K_j}{e}^{b_{j\left(l\right)0}+\sum _{s=1}^S{b}_{j\left(l\right)s}{a}_{is}}}},(k=1,\dots ,K_j).$$

(28)

Using the last category as a baseline in order to make the model identifiable, the parameter for that category are restricted to be 0, i.e., $b_{j\left(K_j\right)0}=b_{j\left(K_j\right)s}=0$, $(j=1,\dots ,J;s=1,\dots ,S)$.The model can be rewritten as:

$${\pi }_{ij\left(k\right)}={\displaystyle \frac{e^{b_{j\left(k\right)0}+\sum _{s=1}^S{b}_{j\left(k\right)s}{a}_{is}}}{1+\sum _{l=1}^{K_j-1}{e}^{b_{j\left(l\right)0}+\sum _{s=1}^S{b}_{j\left(l\right)s}{a}_{is}}}},(k=1,\dots ,K_j-1).$$

(29)

With this restriction we assume that the log-odds of each response (relative to the last category) follows a linear model:

$$log\left({\displaystyle \frac{{\pi }_{ij\left(k\right)}}{{\pi }_{ij\left(K_j\right)}}}\right)=b_{j\left(k\right)0}+\sum _{s=1}^S{b}_{j\left(k\right)s}{a}_{is}=b_{j\left(k\right)0}+{\mathbf{a}}_i^T{\mathbf{b}}_{j\left(k\right)},$$

where $a_{is}$ and $b_{j\left(k\right)s}(i=1,\dots ,I;j=1,\dots ,J;k=1,\dots ,K_j-1;s=1,\dots ,S)$ are the model parameters. In matrix form:

$$\mathbf{O}={\mathbf{1}}_I{\mathbf{b}}_0^T+{\mathbf{AB}}^{\prime },$$

(30)

where ${\mathbf{O}}_{I\times (L-J)}$ is the matrix containing the expected log-odds, defines a biplot for the odds. Although the biplot for the odds may be useful, it would be more interpretable in terms of predicted probabilities and categories. This biplot will be called “nominal logistic biplot”, and it is related to the latent nominal models in the same way as classical linear biplots are related to factor or principal components analysis or binary logistic biplots are related to the item response theory or latent trait analysis for binary data.

The points predicting different probabilities are no longer on parallel straight lines; this means that predictions on the logistic biplot are not made in the same way as in the linear biplots, the response surfaces define now prediction regions for each category as shown in [22]. The nominal logistic biplot is described here in less detail because its geometry is less related to our proposal than linear or binary logistic biplots.

3. Logistic Biplot for Ordinal Data

3.1. Formulation and Geometry

Let ${\mathbf{X}}_{I\times J}$ be a data matrix containing the measures of I individuals on J ordinal variables with $K_j,(j=1,\dots ,J)$ ordered categories each, and let ${\mathbf{P}}_{I\times L}$ the indicator matrix with $L={\sum }_j\left(K_j\right)$ columns. The indicator ${\mathbf{P}}_j$ matrix of size $I\times K_j$ for each categorical variable $X_j$ contains binary indicators for each category and $\mathbf{P}=\left({\mathbf{P}}_1,\dots ,{\mathbf{P}}_J\right)$. Each row of ${\mathbf{P}}_j$ sums 1 and each row of $\mathbf{P}$ sums J. Then $\mathbf{P}$ is the matrix of observed probabilities for each category of each variable. We can also define the cumulative observed probabilities as $p_{ij\left(k\right)}^{*}=1$ if $x_{ij}\le k$ and $p_{ij\left(k\right)}^{*}=0$ otherwise, $(k=1,\dots ,K_j-1)$, that is, the indicators of the cumulative categories. Observe that $p_{ij\left(K_j\right)}^{*}=1$, so we can eliminate the last category for each variable. We organize the observed cumulative probabilities for each variable into a $I\times (K_j-1)$ matrix ${\mathbf{P}}_{\mathbf{j}}^{*}=\left(p_{ij\left(k\right)}^{*}\right)$. Then ${\mathbf{P}}^{*}=(P_1^{*},\dots ,P_J^{*})$ is the $I\times (L-J)$ matrix of observed cumulative probabilities.

Let ${\pi }_{ij\left(k\right)}^{*}=P(x_{ij}\le k)$ be the (expected) cumulative probability that individual i has a value lower than k on the $j-th$ ordinal variable, and let ${\pi }_{ij\left(k\right)}=P(x_{ij}=k)$ the (expected) probability that individual i takes the $k-th$ value on the $j-th$ ordinal variable. Then ${\pi }_{ij\left(K_j\right)}^{*}=P(x_{ij}\le K_j)=1$ and ${\pi }_{ij\left(k\right)}={\pi }_{ij\left(k\right)}^{*}-{\pi }_{ij(k-1)}^{*}$ (with ${\pi }_{ij\left(0\right)}^{*}=0$). A multidimensional (S-dimensional) logistic latent trait model for the cumulative probabilities can be written for ($1\le k\le K_j-1$) as

$${\pi }_{ij\left(k\right)}^{*}={\displaystyle \frac{1}{1+e^{-\left(d_{j\left(k\right)}+{\sum }_{s=1}^S{a}_{is}{b}_{js}\right)}}}={\displaystyle \frac{1}{1+e^{-(d_{j\left(k\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}$$

(31)

where ${\mathbf{a}}_i={(a_{i1},\dots ,a_{iS})}^{\prime }$ is the vector of latent trait scores for the $i-th$ individual and $d_{j\left(k\right)}$ and ${\mathbf{b}}_j={(b_{j1},\dots ,b_{jS})}^{\prime }$ the parameters for each item or variable. Observe that we have defined a set of binary logistic models, one for each category, where there is a different intercept for each but a common set of slopes for all. In the context of IRT, this is known as the “Graded Response Model” or Samejima’s model[28]. The main difference with IRT models is that we don’t have the restriction that the probability of obtaining a higher category must increase along the dimensions. Our variables are not necessarily items of a test, but the models are formally the same for both cases. For the unidimensional case that corresponds to a model with an unique discrimination $b_j$ for all categories and different threshold, boundaries, difficulties or location parameters $d_{j\left(k\right)}$. The two-dimensional cumulative model is shown in Figure 3.

The ${\mathbf{a}}_i$ scores can be represented in a scatter diagram and used to establish similarities and differences among individuals or searching for clusters with homogeneous characteristics, i. e., the representation is like the one obtained from any multidimensional scaling method. In the following we will see that the ${\mathbf{b}}_j$ parameters can also be represented on the graph as directions on the scores space that best predict probabilities and are used to help in searching for the variables or items responsible for the differences among individuals.

In logit scale, the model is

$$logit\left({\pi }_{ij\left(k\right)}^{*}\right)=d_{j\left(k\right)}+\sum _{s=1}^S{a}_{is}{b}_{js}=d_{j\left(k\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j,k=1,\dots ,K_j-1$$

(32)

That defines a binary logistic biplot for the cumulative categories.

In matrix form:

$$logit\left({\mathbf{\Pi }}^{*}\right)=\mathbf{1}{\mathbf{d}}^{\prime }+{\mathbf{AB}}^{\prime }$$

(33)

where ${\mathbf{\Pi }}^{*}=({\mathbf{\Pi }}_1^{*},\dots ,{\mathbf{\Pi }}_J^{*})$ is the $I\times (L-J)$ matrix of expected cumulative probabilities, ${\mathbf{1}}_I$ is a vector of ones and $\mathbf{d}={({\mathbf{d}}_1^T,\dots ,{\mathbf{d}}_J^T)}^T$, with ${\mathbf{d}}_j^T=(d_{j\left(1\right)},\dots ,d_{j(K_j-1)})$, is the vector containing thresholds, $\mathbf{A}={({\mathbf{a}}_1^T,\dots ,{\mathbf{a}}_I^T)}^{\prime }$ with ${\mathbf{a}}_i^T=(a_{i1},\dots ,a_{iS})$ is the $I\times S$ matrix containing the individual scores matrix and $\mathbf{B}={({\mathbf{B}}_1^T,\dots ,{\mathbf{B}}_J^T)}^{\prime }$ with ${\mathbf{B}}_j={\mathbf{1}}_{K_j-1}\otimes {\mathbf{b}}_j^T$ and ${\mathbf{b}}_j^T=(b_{j1},\dots ,b_{jS})$, is the $(L-J)\times S$ matrix containing the slopes for all the variables. This expression defines a biplot for the odds that will be called ordinal logistic biplot. Each equation of the cumulative biplot shares the geometry described for the binary case [33], moreover, all curves share the same direction when projected on the biplot. The set of parameters $\left\{d_{j\left(k\right)}\right\}$ provide a different threshold for each cumulative category, the second part of (32) does not depend on the particular category, meaning that all the $K_j-1$ curves share the same slopes. In the following paragraphs we will obtain the geometry for the general case an algorithm to perform the calculations.

The expected probability of individual i responding in category k to item j, with $(k=1,\dots ,K_j)$, that we denote by ${\pi }_{ij\left(k\right)}=P(x_{ij}=k)$ must be obtained by subtracting cumulative probabilities:

$${\pi }_{ij\left(k\right)}={\pi }_{ij\left(k\right)}^{*}-{\pi }_{ij(k-1)}^{*}$$

then using the equations in (31):

$$\begin{array}{cc}\hfill {\pi }_{ij\left(1\right)}& =P(x_{ij}=1)={\displaystyle \frac{1}{1+e^{-(d_{j\left(1\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}\hfill \\ \hfill {\pi }_{ij\left(k\right)}& =P(x_{ij}=k)=P(x_{ij}\le k)-P(x_{ij}\le (k-1))\hfill \\ \hfill & ={\displaystyle \frac{1}{1+e^{-(d_{j\left(k\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}-{\displaystyle \frac{1}{1+e^{-(d_{j(k-1)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}\hfill \\ \hfill & ={\displaystyle \frac{e^{-\left({\mathbf{a}}_i^T{\mathbf{b}}_j\right)}(e^{-d_{j(k-1)}}-e^{-d_{j\left(k\right)}})}{(1+e^{-(d_{j\left(k\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)})(1+e^{-(d_{j(k-1)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)})}},1<k<K_j\hfill \\ \hfill {\pi }_{ij\left(K_j\right)}& =P(x_{ij}=K_j)=1-{\displaystyle \frac{1}{1+e^{-(d_{j(K_j-1)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}\hfill \end{array}$$

(34)

If the row scores where known, obtaining the parameters of the model in 34 is equivalent to fitting a proportional odds model using each item as a response and the row scores as regressors. The response surfaces for such a model are shown in Figure 4. Although the response surfaces are no longer sigmoidal, the level curves are still straight lines, so the set of points on the representation (generated by the columns of $\mathbf{A}$) predicting a particular value for the probability of a category lie on a straight line, and different probabilities for all the categories of a particular variable or item lie on parallel straight lines. A perpendicular to all those lines can be used as "biplot axis" in the sense of [19] and is the direction that better predicts the probabilities of all the categories, that is, projecting any individual point onto that direction, we should obtain an optimal prediction of the category probabilities. As all the categories share the same biplot direction, it would be very difficult to place a different graded scales for each and we will represent just the line segments in which the probability of a category is higher then the probability of the rest. That will result, except for some pathological cases, in as many segments as categories ($K_j$), separated by $K_j-1$ points in which the probabilities of two (contiguous) categories are equal. See Figure 5 in which we show the parallel lines representing the points that predict equal probabilities for two contiguous categories and a line, perpendicular to all, that is the "biplot axis". The three parallel lines divide the space spanned by the columns of $\mathbf{A}$ into four regions, each predicting a particular category of the variable. For a biplot representation we don’t need the whole set of lines but just the "axis" and the points on it, intersecting the boundaries of the prediction regions.

3.2. Obtaining the Biplot Representation

So, if we denote $(x,y)$ one of those intersection points, it must be on the biplot direction, that is,

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

(35)

and the probability of two, possibly contiguous, categories (for example l and m) at this point, must be equal,

$${\pi }_{j\left(l\right)}={\pi }_{j\left(m\right)}({\pi }_{j\left(l\right)}^{*}-{\pi }_{j(l-1)}^{*}={\pi }_{j\left(m\right)}^{*}-{\pi }_{j(m-1)}^{*}).$$

(36)

We have omitted index i because probabilities are for a general point and not for a particular individual. Using the condition in 35 we can rewrite the cumulative probabilities (or its $logit$) as

$$logit\left({\pi }_{j\left(k\right)}^{*}\right)=d_{j\left(k\right)}+x{b}_{j1}+y{b}_{j2}=d_{j\left(k\right)}+z$$

(37)

with

$$z=x\left({\displaystyle \frac{b_{j1}^2+b_{j2}^2}{b_j}}\right)$$

(38)

Changing the values of z we can obtain the probabilities of each category along the biplot axis. So, finding the point $(x,y)$ is equivalent to find the values of z in which 36 holds. From those values the original point is obtained solving for x in 38 and then calculating y from 35.

There are some pathological cases in which the probability of one or several categories are never higher than the probability of the rest, in such cases we say that the category is “hidden” or “never predicted” and the number of separating points will be lower than $K_j-1$. Those pathological cases have to be taken into account when calculating the intersection points.

Figure 6. Probability curves for a variable with 6 categories in which two (4 and 5) are hidden or never predicted.

Figure 6. Probability curves for a variable with 6 categories in which two (4 and 5) are hidden or never predicted.

The existence of abnormal cases means that, not just contiguous, but any pair of categories may have to be compared. Then, many comparisons are possible because the equations are different for each case

• 1-2
• 1-$l(l<K_j)$
• 1-$K_j$
• l-$K_j(l>1)$
• l-$(l+1)$ with $l>1$
• l-j with $j>(l+1),l>1$
• $(K_j-1)$-$K_j$

For example, in case (3) 1-$K_j$, is simple to deduce that

$$z={\displaystyle \frac{-(d_{j(K_j-1)}+d_{j\left(1\right)})}{2}}.$$

Cases (1), (3), (5) and (7) are simple. In the other 3 combinations we have to solve a quadratic equation to obtain the intersection points. For example, in case (2), the first with the $l-th$ categories, we have to solve ${\pi }_{j\left(1\right)}={\pi }_{j\left(l\right)}$, that is

$${\displaystyle \frac{1}{1+e^{-(d_{j\left(1\right)}+z)}}}={\displaystyle \frac{e^{-z}(e^{-d_{j(l-1)}}-e^{-d_{j\left(l\right)}})}{(1+e^{-(d_{j\left(l\right)}+z)})(1+e^{-(d_{j(l-1)}+z)})}}$$

Taking

$$w=e^{-z}$$

we have to solve the quadratic equation

$$\alpha w^2-\beta w-1=0$$

with $\alpha =(e^{-(d_{j\left(1\right)}+d_{j(l-1)})}-e^{-(d_{j\left(1\right)}+d_{j\left(l\right)})}-e^{-(d_{j(l-1)}+d_{j\left(l\right)})})$ and $\beta =2e^{-d_i}$.

If the roots of the equation are both negative, curves don’t intersect. If it has a positive root, we can calculate the intersection points solving for w and then reversing the transformations to obtain $(x,y)$. In a similar way we can calculate the intersection points for cases (4) i-$K_j(i>1)$ and (6) i-j with $j>(i+1)$.

A procedure to calculate the representation of an ordinal variable on the biplot would be as follows:

• Calculate the biplot axis with equation $y={\displaystyle \frac{b_{j2}}{b_{j1}}}x$.
• Calculate the intersection points z and then $(x,y)$ of the biplot axis with the parallel lines used as boundaries of the prediction regions for each pair of categories, in the following order:

  $$\begin{array}{cc}\hfill {\pi }_{j\left(1\right)}& ={\pi }_{j\left(2\right)}\hfill \\ \hfill {\pi }_{j(l-1)}& ={\pi }_{j\left(l\right)},1<l<(K_j-1)\hfill \\ \hfill {\pi }_{j(K_j-1)}& ={\pi }_{j\left(K_j\right)}\hfill \end{array}$$
• If the values of z are ordered, there are not hidden categories and the calculations are finished.
• If the values of z are not ordered we can do the following:

  (a)

  Calculate the z values for all the pairs of curves, and the probabilities for the two categories involved.

  (b)

  Compare each category with the following, the next to represent is the one with the highest probability at the intersection.

  (c)

  If the next category to represent is $K_j$ the process is finished. It not go back to the previous step, starting with the new category.

A simpler numeric procedure could also be developed to avoid the explicit solution of the equations.

• Calculate the predicted category for a set of values for z. For example a sequence from -6 to 6 with small steps (0,01 or 0,01 for example). (The precision of the procedure can be changed with the step)
• Search for the z values in which the prediction changes from one category to another.
• Calculate the mean of the two z values obtained in the previous step and then the $(x,y)$ values. Those are the points we are searching for.

Hidden categories are the ones with zero frequencies in the predictions obtained by the algorithm.

3.3. Parameter Estimation Based on an Alternated Gradient Descent Algorithm on the Cumulative Probabilities

The alternated algorithm described in [33], can be easily extended replacing binary logistic regressions by ordinal logistic regressions. The problem with this approach is that the parameters for the individuals can not be estimated when the individual has 0 or 1 in all the variables for the binary case, or all the responses are at the baseline category for the ordinal case. There is also some problems when perfect probability predictions are present and then there is a perfect separation of some of the categories. The procedure would be similar to the joint maximum likelihood estimation of Item Response Theory (IRT) models. An algorithm for the cumulative probabilities has been developed in [7]. We propose an alternative algorithm based on an alternated gradient descent procedure and the biplot representation described before.

We use cumulative probabilities together with a recursive procedure similar to the one used for binary (data) based on gradient descent methods.

The cost function based on cumulative probabilities is

$$L=-\sum _{i=1}^I\sum _{j=1}^J\sum _{k=1}^{K_j}{p}_{ij\left(k\right)}^{*}log\left({\pi }_{ij\left(k\right)}^{*}\right).$$

(39)

where the expected cumulative probabilities are

$${\pi }_{ij\left(k\right)}^{*}={\displaystyle \frac{1}{1+e^{-\left(d_{j\left(k\right)}+{\sum }_{s=1}^S{a}_{is}{b}_{js}\right)}}}={\displaystyle \frac{1}{1+e^{-(d_{j\left(k\right)}+{\mathbf{a}}_i^T{\mathbf{b}}_j)}}}$$

(40)

We will use a recursive procedure as in the binary case, estimating first the thresholds $d_{jk}$ and fixing them the parameters for the dimensions one at a time. That would be equivalent to subtract the mean for continuous data.

The gradients for the binary case are easily generalized to the ordinal case by using the cumulative probabilities. The updates are then

$$\begin{array}{c}\hfill d_{j\left(k\right)}\Leftarrow d_{j\left(k\right)}-\alpha {\displaystyle \frac{\partial L}{\partial d_{j\left(k\right)}}}=d_{j\left(k\right)}-\alpha \sum _{i=1}^I({\pi }_{ij\left(k\right)}^{*}-p_{ij\left(k\right)}^{*}),\end{array}$$

(41)

$$\begin{array}{c}\hfill a_{is}\Leftarrow a_{is}-\alpha {\displaystyle \frac{\partial L}{\partial a_{is}}}=a_{is}-\alpha \sum _{j=1}^J{b}_{js}\sum _{k=1}^{K_j}({\pi }_{ij\left(k\right)}^{*}-p_{ij\left(k\right)}^{*}),\end{array}$$

(42)

$$\begin{array}{c}\hfill b_{js}\Leftarrow b_{js}-\alpha {\displaystyle \frac{\partial L}{\partial b_{js}}}=b_{js}-\alpha \sum _{i=1}^I{a}_{is}\sum _{k=1}^{K_j}({\pi }_{ij\left(k\right)}^{*}-p_{ik\left(k\right)}^{*}).\end{array}$$

(43)

Let ${\mathbf{a}}_s={(a_{1s},\dots ,a_{Is})}^T$ and ${\mathbf{b}}_s={(b_{1s},\dots ,b_{Js})}^T$ the vectors containing the row and column parameters for dimension s. Let $\mathbf{Z}$ an $(L-J)\times J$ indicator matrix in which, the $j-th$ column, takes the value 1 for the rows corresponding to categories of the $j-th$ variable, and 0 elsewhere. Then the columns of the matrix

$$\mathbf{R}=({\mathbf{\Pi }}^{*}-{\mathbf{P}}^{*})\mathbf{Z}$$

(44)

contain $r_{ij}=\sum _{k=1}^{K_j}({\pi }_{ij\left(k\right)}^{*}-p_{ij\left(k\right)}^{*})$ for each variable. The updates are written in matrix form

$$\begin{array}{c}\hfill \mathbf{d}\Leftarrow \mathbf{d}-\alpha {({\mathbf{\Pi }}^{*}-{\mathbf{P}}^{*})}^T{\mathbf{1}}_I,\end{array}$$

(45)

$$\begin{array}{c}\hfill {\mathbf{a}}_s\Leftarrow {\mathbf{a}}_s-\alpha \mathbf{R}{\mathbf{b}}_s,\end{array}$$

(46)

$$\begin{array}{c}\hfill {\mathbf{b}}_s\Leftarrow {\mathbf{b}}_s-\alpha {\mathbf{R}}^T{\mathbf{a}}_s.\end{array}$$

(47)

For large data matrices it may be convenient to summarize the data into a matrix $\mathbf{X}$ (and the corresponding $\mathbf{P}$) containing unique patterns in the data and a vector with the frequencies of each pattern $\mathbf{f}$. The equations of the gradient descent method can be written now as:

$$\begin{array}{c}\hfill \mathbf{d}\Leftarrow \mathbf{d}-\alpha {({\mathbf{\Pi }}^{*}-{\mathbf{P}}^{*})}^T\mathbf{f},\end{array}$$

(48)

$$\begin{array}{c}\hfill {\mathbf{a}}_s\Leftarrow {\mathbf{a}}_s-\alpha \mathbf{R}{\mathbf{b}}_s,\end{array}$$

(49)

$$\begin{array}{c}\hfill {\mathbf{b}}_s\Leftarrow {\mathbf{b}}_s-\alpha {\mathbf{R}}^T(\mathbf{f}\odot {\mathbf{a}}_s).\end{array}$$

(50)

We can organize the calculations in an alternated algorithm that alternatively calculates parameters for rows and columns for each dimension s fixing the parameters already obtained for each previous dimensions. Before that, the constants $d_{j\left(k\right)}$ have to be fitted separately. The algorithm will produce decreasing values of the cost function and eventually will converge at least to a local minima. A good strategy would be starting at different points and select the solution with lower cost.

The procedure would be as follows:

Algorithm 2 Algorithm to calculate the Components for Ordinal Data

1: procedureP-Ordinal-Components($\mathbf{P},S$) 2:     Choose $\alpha$ 3:     Init: ${\mathbf{b}}_0=random$ 4:     repeat 5:         Update: $\mathbf{d}$ with Equation 45 6:         Update : ${\mathbf{\Pi }}^{*}=\left({\pi }_{ij}^{*}\right)$ with Equation 40 7:     until ${\mathbf{b}}_0$ do not change 8:     for $s=1\to S$ do 9:         Init: ${\mathbf{a}}_{\left(s\right)}\leftarrow random$ (or any other choice) 10:         repeat 11:            repeat 12:                Update: ${\mathbf{b}}_{\left(s\right)}$ with Equation 46 13:                Update : ${\mathbf{\Pi }}^{*}=\left({\pi }_{ij}^{*}\right)$ with Equation 40 14:            until ${\mathbf{b}}_{\left(s\right)}$ do not change 15:            repeat 16:                Update: ${\mathbf{a}}_{\left(s\right)}$ with Equation 47 17:                Update : ${\mathbf{\Pi }}^{*}=\left({\pi }_{ij}^{*}\right)$ with Equation 40 18:            until ${\mathbf{a}}_{\left(s\right)}$ do not change 19:            $L\leftarrow (-1){\mathbf{1}}_I^T((\mathbf{P}\odot \mathbf{\Pi })+((1-\mathbf{P})\odot (1-\mathbf{\Pi }))){\mathbf{1}}_J$ 20:         until L does not change 21:     end forreturn ${\mathbf{b}}_0$, $\mathbf{B}=\left[{\mathbf{b}}_{\left(1\right)},\dots ,{\mathbf{b}}_{\left(S\right)}\right]$, $\mathbf{A}=\left[{\mathbf{a}}_{\left(1\right)},\dots ,{\mathbf{a}}_{\left(S\right)}\right]$ 22: end procedure

In practice, the choice of $\alpha$ can be avoided by using a pre-programmed optimization routine. This algorithm will be implemented in the new version of package MultBiplotR ([34]) developed for the R language ([27]).

3.4. Factorization of the Polychoric Correlation Matrix

There is another way of obtaining the parameters: The factorization of the polychoric correlations matrix.

The idea is that the observed categorical responses are discretized versions of a continuous process.

Consider an ordinal variable $X_j$ with $K_j$ categories. This variable is assumed to arise from an underlying continuous variable $X_j^{*}$ that follows a standard normal distribution. There are $K_j-1$ thresholds, ${\tau }_{j\left(1\right)},{\tau }_{j\left(2\right)},\dots ,{\tau }_{j(K_j-1)}$, which divide the continuous variable $X^{*}$ into $K_j$ ordinal categories.

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

$$X_j=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{X}_j^{*}\le {\tau }_{j\left(1\right)}\hfill \\ 2\hfill & \mathrm{if}{\tau }_{j\left(1\right)}<X_j^{*}\le {\tau }_{j\left(2\right)}\hfill \\ 3\hfill & \mathrm{if}{\tau }_{j\left(2\right)}<X_j^{*}\le {\tau }_{j\left(3\right)}\hfill \\ ⋮\hfill & ⋮\hfill \\ k\hfill & \mathrm{if}{X}_j^{*}>{\tau }_{j,(K_j-1)}\hfill \end{array}\right.$$

The polychoric correlations are the correlations among the $X_j^{*},j=1,\dots ,J$. Let $\mathbf{R}$ a $J\times J$ matrix containing the polychoric correlations among the J ordinal variables and let ${\tau }_{j\left(k\right)},j=1,\dots ,J;K=1,\dots ,K_j$ the thresholds.

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

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

(51)

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

It can be shown the there is a close relation between the factor model in 51 and the model in 31. You can find the details in [23].

If the factor model loadings $\mathbf{\Lambda }=({\lambda }_{js},j=1,\dots ,J;s=1,\dots ,S)$ and he thresholds ${\tau }_{j\left(k\right)},j=1,\dots ,J;K=1,\dots ,K_j$ are known, the parameters for the variables in our model (31) can be calculated as:

$$d_{j\left(k\right)}={\tau }_{j\left(k\right)}{\left(1-\sum _{s=1}^S{\lambda }_{js}^2\right)}^{-1/2}$$

(52)

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

(53)

The remaining parameters for te individuals $a_{is}$ can then be estimated using the gradient descent method with the cost function in (39) with fixed column parameters and equations (), () or ().

On the other hand, if we have the biplot parameters as in 31, we can obtain the factor model in 51 as:

$${\tau }_{j\left(k\right)}=d_{j\left(k\right)}{\left(1+\sum _{s=1}^S{b}_{js}^2\right)}^{-1/2}$$

(54)

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

(55)

Then, if we the biplot parameters, we can also interpret the dimensions as in a traditional factor analysis sing factor loadings and communalities.

3.5. Goodness of Fit

The log-likelihood can be used as a measure of the overall goodness of fit, specially for comparison of different models containing, for example, a different number of dimensions. Statistical tests similar to those proposed in the context of IRT models could be used here, in particular, tests based on ordinal logistic regressions. [26] propose tests based on logistic regressions for binary data that could easily be extended to ordinal data. On one hand, the statistical tests for this situations have many problems and on the other hand we see the procedure more as a descriptive exploratory model, so we are less interested in global statistical tests and more in goodness of fit indices or tests for each separate variable. For IRT models, the items are closely related to one or several latent dimensions and all are useful to describe the latent factors, so there must be an adequate overall goodness of fit. In a more general exploratory situation some of the variables may not be useful for describing the problem, inducing a lower overall fit that can lead to miss-interpretations. [8], conduct a simulation study in which some irrelevant variables and noise are added to a known structure, showing that the known structure is recovered even when the overall goodness of fit is not high. Using fit indices for each variable, they are able to identify the relevant and eliminate the irrelevant variables. A similar result can be found in [12].

Several tests and goodness of fit indices can be defined for each variable considering that has been fitted using an ordinal logistic regression with proportional odds. Here we use the likelihood ratio test to compare the model with constant probability (no latent dimensions) with the complete model in exactly the same way as in the standard logistic regression. The test should be interpreted here with care because the latent variables are also estimated inside the procedure and its variation is not considered explicitly in the test; nevertheless is an indication of the significance of the variable to describe the data. For classical linear biplots no such tests are usually performed but goodness of fit indices are calculated. In [18] predictivities as the percentage of the variance of each variable explained by the dimensions, i. e., is a measure of the prediction accuracy on the biplot, is defined. The predictivity is used as a measure of goodness of fit measure by the package BiplotGUI ([24]) In the context of Correspondence Analysis those quantities were called relative contributions of the axis to the elements (variables or individuals) (see [3,20] or [3]) extended to biplots by [15] or [20]. Contributions of the dimensions to the variables can be interpreted as squared correlations between observed variables and dimensions.

From another point of view, the predictivity, is the coefficient of determination $R_j^2$ for the regressions in 8. For ordinal responses, a pseudo $R_j^2$ as the Cox-Senell or Nagelkerke can be used.

Another possible measure of goodness of fit is the percentage of correct classifications. We calculate two percentages, one for the cumulative probabilities and another one for the original categories. The observed and predicted items are both ordinal so, for measuring the association between them we use the weighted kappa coefficient.

4. An Empirical Study

4.1. Dataset

In 2008, for a set of 26 countries worldwide, among which was Spain, setting 2006 as reference year and following the guidelines set by the Organization for Cooperation and Development(OECD), the department of statistics of UNESCO and Eurostat (Statistical Office of the European Union), began to do surveys to people that had obtained a PhD degree, and that therefore are doctorates, with the objective of having a clearer information about their characteristics. Most of the pioneer countries belonged the European Union, although members of the OECD, as USA or Australia, also participated. In the Spanish case, it was the National Institute of Statistics (Instituto Nacional de Estadistica - INE) which focused all efforts to carry out this new operation with the objective that the availability of information in this field had continuity in time. Thus, the so-called “Survey on Human Resources in Science and Technology” was established as part of the general plan of science and technology statistics carried out by the European Union Statistics Office (Eurostat). This need for information is evident in the European Regulation 753/2004 on Science and Technology, which specifies the production of statistics on human resources in science and technology.

Surveys on doctorates (CDH: Careers of doctorate holders) try to measure specific demographic aspects related to employment, so that the level of investigation of this group, the professional activity carried out, the satisfaction with their main job, the international mobility and the income of this group can be quantified in Spain.

The study focused on all the doctorate holders resident in Spain, younger than 70, that obtained their degree in a Spanish university, both public or private, between 1990 and 2006. The frame of this statistical operation was a directory of doctorate holders provided to the National Statistic Institute by the University Council, which includes all the persons who have defended a doctoral thesis in any Spanish university, according with their electronic databases, which was comprised of approximately 80000 people. Doctors belong to the level 6 of the international classification of education ISCED-97. This level is reserved for tertiary programs which lead to the award of an advanced research qualification and devoted to advanced study and original research and not based on course-work only.

As for the sampling design, a representative sample was designed for each region at NUTS-2 level. The NUTS classification (Nomenclature of Territorial Units for Statistics) is a hierarchical system for dividing up the economic territory of the EU for different purposes. (see Eurostat) , using a sampling with equal probabilities. The doctors were grouped according to their place of residence, and the selection has been done independently in each region by equal probability with random start systematic sampling. A sample of 17000 doctors was selected. The sample has been distributed between the regions assigning the 50% in an uniform way and the rest proportional to the size of them, measured in number of doctorate holders that have their residence in those regions.

The INE used a questionnaire harmonized at European level, structured in several modules, that can be found at the website of the Institute (INE). As a result of the collection process it was obtained a response rate at the national level of 72%.

We have 12193 doctorate’s answers to the questionnaire in Spain to develop this study. We will focus our attention on the module C (Employment situation) and specifically, in the subsection C.6.4, that tries to find out the level of satisfaction of the doctorate holders in aspects related to their main job. This question has several points of interest coded on a likert scale from 1 to 4 (see Figure 7).

Each item will be considered as ordinal, then we have 11 items or variables in total. For an easier interpretation codification the answers to each question have been inverted in such a way that higher numbers are associated to better satisfaction. Then the final codification is: 1-Very dissatisfied, 2-Somewhat dissatisfied, 3-Somewhat satisfied and 4-Very satisfied.

4.2. Results

Using the alternated algorithm to estimate the parameters of the two-dimensional model, we have obtained the fit indicators in Table 1, the factor loadings and communalities in Table 2 and the explained variances in Table 3.

All the variables fit reasonably well to the model as indicated for the different fit measures. All the $pseudo-R^2$ are reasonably high. For example, Nagelkerke pseudo $R^2$ values range from 0.36 to 0.65 that would be considered by most authors as reasonably high. The same would apply for the other pseudo-$R^2$ measures.

The percents of correct classifications for he cumulative probabilities are quite high , from 85.40% to 91.65%, with a global percent of 88.03%. For the original ordinal values is smaller, from 56.06% to 71.53%, with a global rate of 64.95%. We have to take into account that we have optimized the cumulative distributions rather than he original ordinal values.

The Kappa coefficients are not very high, probably due to the fact that we have optimized the cumulative distribution rather than the original values. In any case, all the indices are just indicators of what variables fit better to the model, that is, what variables have a higher relation to the latent factors.

Analyzing the interpretation of two dimensions or factors, using their loadings, the first has higher weights for the variables, opportunities for advancement, degree of independence, intellectual challenge, level of responsibility and contribution to society, which are features associated to intellectual satisfaction. The second factor presents higher values for salary, benefits, job security and working conditions, all related with economical and work satisfaction. Then we have two main almost independent factors, the first related to the intellectual aspects and the second to the job conditions. The job location and social status variables have similar loadings in both factors meaning that are related to both. The angle between the variables defining each factor is close to 90°, which means that satisfaction with income and work seems to be nod correlated with intellectual, something that also appears in other European countries, as is evident in the Austrian case [29].

The variance explained by two factors is 65.94% as shown in Table 3.

The biplot of ordinal data can be seen in Figure 8, in which points represent to individuals and the directions to variables. Together with the measures presented before, graphical representations are useful tools for interpreting patterns in data; biplots also can simultaneously find important features in individuals and variables and relations among them.

In general angles between directions are interpreted as correlations like in classical biplots, actually polychoric correlations among ordinal variables in this case. Acute small angles are interpreted as high positive correlation, almost straight angles as high negative correlation and right angles as no correlation.

Some of the variables have a similar behavior, such as level of responsibility, challenge or contribution to society, because angles between its directions are small. Although there are groups of variables whose directions are very similar, the position of the categories are quite different from each other within those groups. This happens with opportunities for advancement and degree of independence.

Markers for individuals are represented by small dots and not labeled because the plot would be too crowded and, in this case, a particular individual is not studied. In general, distances among points representing individuals are related to their similarity. The closer two points are, the more similar. Normally this allows to find clusters of similar individuals and the variables or items responsible for groupings. Cluster may also be defined by external nominal variables to investigate if the factors are related to it.

The projection of an individual into a direction for a variable predicts categories of the variables. The points separating the prediction of the categories ar marked on each direction (variable, item) with numbers. For example, points marked with number 1 are the threshold to change the prediction from the first to the second category, 2 is the threshold from the second to the third an so on.

The figure can be improved using different colors for each variable as in Figure 9. The colors may help identifying marks for each variable.

All the variables have three marks except job security that has 2 (1 and 2) meaning that the category 3 is never predicted or hidden, that is, the probability oh having the second category is never higher than probabilities for other categories.

The plot reflects the factor structure of the data set in a graphical way. To illustrate how the representation gets the structure of a particular variable we can add to the plot clusters with points with different observed categories. Figure 10 such a representation for variable Chalenge

We notice that the cluster centers are arranged nearly along the direction of Challenge, indicating that this axis effectively reflects the behavior of that variable. A similar analysis can be carried out for the other variables, confirming that the plots capture information from all of them, though to varying degrees. The fit measures provide an indication of how well each variable is represented.

The final plot can be also used to check the behavior of different groups of individuals defined by external nominal variables, for example, males and females to investigate the relation between satisfaction and sex. Figure 11.

The positions of men and women, measured as centroids of coordinates on the biplot for each group, are virtually indistinguishable (Figure 11), meaning that there is a very small difference between men and women in the way they perceive satisfaction with doctorate degrees.

5. Software Note

The procedures of his paper will be added to the package MultBiplotR ([34]) developed in the R language ([27]).