Biplot - CoDA

Análisis de Datos Composicionales (CoDA)
Maestría en Ciencias - Estadística, UNAL

Authors

Camilo Gómez

Sergio A. Prieto

1 Biplot usual

Un biplot es una herramienta gráfica que permite representar las filas y columnas de una matriz \(\mathbf{X}\) de tamaño \(n \times p\) como puntos en un espacio Euclidiano de menor dimensión donde \(n\) representa el número de individuos y \(p\) el número de variables.

En la práctica, la matriz \(\mathbf{X}\) se suele transformar ya sea centrando los datos, estandarizándolos, etc., obteniendo una matriz \(\mathbf{Z}\) de rango \(r\). Toda matriz rectangular se puede descomponer de infinitas formas como el producto de dos matrices \(\mathbf{F}\) y \(\mathbf{G}\), de tamaños \(n\times r\) y \(p \times r\) respectivamente, tal que \[\mathbf{Z} = \mathbf{F} \mathbf{G}^\top,\]

donde las filas de \(\mathbf{F}\) y las filas de \(\mathbf{G}\) corresponden a las coordenadas de los \(n\) puntos para las filas y \(p\) puntos para las columnas en el espacio Euclidiano de dimensión \(r\).

Una selección en particular de dichas matrices, se obtiene de la descomposición en valores singulares (SVD), que además cumple con ser la mejor aproximación de rango \(q^*\) en el sentido en el que \(\mathbf{Z}^*_q\) minimiza \[\lVert \mathbf{Z} - \mathbf{Y} \rVert^2 = \sum_i \sum_j (z_{ij} - y_{ij})^2,\] sobre todas las matrices \(\mathbf{Y}\) de rango \(q^*\) (Teorema de Eckart-Young).

Para obtener la descomposición SVD de \(\mathbf{Z}\), es necesario formar las matrices de vectores propios de \(\mathbf{Z}\mathbf{Z}^\top\) (que se denota como \(\mathbf{U}\)) y de \(\mathbf{Z}^\top \mathbf{Z}\) (que se denota como \(\mathbf{V}\)) y sus \(r\) valores propios \(\lambda_1\geq\lambda_2 \geq\cdots \geq \lambda_r > 0\), de forma que \(\mathbf{U}\) y \(\mathbf{V}\) tienen \(r\) columnas ortonormales. Luego, \(\mathbf{Z}\) se puede escribir como \[\mathbf{Z} = \mathbf{U} \mathbf{\Gamma} \mathbf{V}^\top,\] donde \(\mathbf{\Gamma} = \operatorname{diag}\{{\gamma_1}, \ldots, \gamma_r\}\), con \(\gamma_i=\sqrt{\lambda_i}\), \(i=1,\ldots,r\).

Con el objetivo de reducir la dimensionalidad, se consideran únicamente las componentes que explican la mayor variabilidad de los datos. Si se retienen \(q\) valores singulares, entonces la proporción de variabilidad explicada por estas \(q\) componentes se calcula como \[\frac{\sum_{i=1}^q \lambda_i}{\sum_{j=1}^r \lambda_j}.\] Normalmente, se toma \(q=2\) o \(q=3\). Si la aproximación de \(\mathbf{Z}\) mediante \(\mathbf{Z}_q\) es lo suficientemente buena, luego el biplot representará correctamente las características presentes en los datos.

De esta manera, las matrices \(\mathbf{F}\) y \(\mathbf{G}\) tienen la forma: \[\mathbf{F} = \mathbf{U} \mathbf{\Gamma}^\alpha, \qquad \mathbf{G} = \mathbf{\Gamma}^{1-\alpha} \mathbf{V}^\top ,\] tomando \(q=2\) están dadas por \[\mathbf{F} = \begin{pmatrix} \gamma_1^\alpha \mathbf{u}_1 & \gamma_2^\alpha \mathbf{u}_2 \end{pmatrix} , \qquad \mathbf{G} = \begin{pmatrix} \gamma_1^{1-\alpha} \mathbf{v}_1 & \gamma_2^{1-\alpha} \mathbf{v}_2 \end{pmatrix} ,\] donde \(\alpha\in[0,1]\) puede tomar distintos valores dependiendo del aspecto que se quiera resaltar de la matriz de datos.

2 Biplot Composicional

Sea \(\mathbf{X}\) la matriz de datos composicionales con \(n\) filas y \(D\) columnas, donde cada columna corresponde a una parte \(X_1,\ldots, X_D\). Se hallan las coordenadas clr de la matriz transformada con filas \[\mathbf{z}_i = \mathrm{clr}\left( \mathbf{x}_i \ominus \mathbf{g} \right)\] Véase que la transformación clr centra sobre las filas. Sea \[z_{ij}^*=\ln{\frac{x_{ij}}{g(X)}} = \ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})}\] fijando \(i\) y calculando la media de los \(z_{ij}^*\) \[\begin{aligned} \frac{1}{D}\sum_{j=1}^D z_{ij}^* &= \frac{1}{D}\sum_{j=1}^D \left(\ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})} \right)\\ &= \frac{1}{D} \left[\sum_{j=1}^D \ln{(x_{ij})} - \sum_{j=1}^D \left(\frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})}\right)\right]\\ &=\frac{1}{D} \left[\sum_{j=1}^D \ln{(x_{ij})} - D \left(\frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})}\right)\right]\\ &=\frac{1}{D} \left[\sum_{j=1}^D \ln{(x_{ij})} - \sum_{j=1}^D \ln{(x_{ij})}\right] = 0. \end{aligned}\]

Con lo que se tiene que la media de los elementos de \(Z\) por fila es cero.

Solo resta centrar la matriz de composiciones sobre las columnas, tal que \[\begin{aligned} z_{ij} &= \ln{\frac{x_{ij}}{g(X)}} - \frac{1}{n}\sum_{i=1}^n clr{(x_{ij})} \\ &= \ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})} - \frac{1}{n} \sum_{i=1}^n \left[ \ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})} \right] \\ &= \ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})} - \frac{1}{n} \sum_{i=1}^n \ln{(x_{ij})} +\frac{1}{nD} \sum_{i=1}^n \sum_{j=1}^D \ln{(x_{ij})} \end{aligned}\]

Por tanto se tiene que la matriz \(\mathbf{Z}\) está doblemente centrada, por filas y columnas.
Véase que en efecto se encuentra centrada sobre las columnas, para ello fijando la columna \(j\) \[\begin{aligned} \frac{1}{n} \sum_{i=1}^n z_{ij} &= \frac{1}{n} \sum_{i=1}^n \left[ \ln{(x_{ij})} - \frac{1}{D}\sum_{j=1}^D \ln{(x_{ij})} - \frac{1}{n} \sum_{i=1}^n \ln{(x_{ij})} +\frac{1}{nD} \sum_{i=1}^n \sum_{j=1}^D \ln{(x_{ij})} \right]\\ &= \frac{1}{n} \left[ \sum_{i=1}^n \ln{(x_{ij})} - \frac{1}{D}\sum_{i=1}^n\sum_{j=1}^D \ln{(x_{ij})} - \frac{1}{\cancel{n}} \cancelto{\cancel{n}}{\sum_{i=1}^n}\sum_{i=1}^n \ln{(x_{ij})} +\frac{1}{\cancel{n}D} \cancelto{\cancel{n}}{\sum_{i=1}^n}\sum_{i=1}^n \sum_{j=1}^D \ln{(x_{ij})} \right]\\ &= \frac{1}{n} \left[\cancel{\sum_{i=1}^n \ln{(x_{ij})}} - \bcancel{\frac{1}{D}\sum_{i=1}^n\sum_{j=1}^D \ln{(x_{ij})}} - \cancel{\sum_{i=1}^n \ln{(x_{ij})}} + \bcancel{\frac{1}{D} \sum_{i=1}^n \sum_{j=1}^D \ln{(x_{ij})}} \right]\\ &= \frac{1}{n}[0] = 0 \end{aligned}\]

Con lo que queda probado que la matriz \(Z\) está centrada con respecto a filas y columnas.
A continuación, se sigue el procedimiento descrito anteriormente, donde se construyen las matrices \(\mathbf{U}\) con los vectores propios de \(\mathbf{Z}\mathbf{Z}^\top\), \(\mathbf{V}\) con los vectores propios de \(\mathbf{Z}\top\mathbf{Z}\) y \(\mathbf{\Gamma}\) con los \(r\leq \min\{ D-1,n \}\) valores singulares positivos de \(\mathbf{Z}\mathbf{Z}^\top\) o \(\mathbf{Z}\top\mathbf{Z}\).

  • Las matrices \(\mathbf{U}\) y \(\mathbf{V}\) son ortogonales, de forma que \(\mathbf{U}^\top \mathbf{U} = \mathbf{I}_r\) y \(\mathbf{V}^\top\mathbf{V} = \mathbf{I}_r\).

  • Cuando \(\mathbf{Z}\) se forma con las coordenadas clr centradas, las filas suman 0 y entonces \(r\leq D-1\).

  • Cada fila de la matriz \(\mathbf{V^\top}\) es la coordenada clr de algún elemento de una base ortonormal en el simplex (cargas o loadings).

  • Las filas de la matriz producto \(\mathbf{U} \operatorname{diag}\{{\gamma_1}, \ldots, \gamma_r\}\) contienen las coordenadas de cada dato composicional con respecto a la base ortonormal descrita por \(\mathbf{V}^\top\) (scores en el PCA). Es decir, contienen las coordenadas ilr de los datos composicionales centrados.

Al igual que antes, tomando \(q=2\) en los productos \(\mathbf{Z} = \mathbf{U} \mathbf{\Gamma}^\alpha, \mathbf{G} = \mathbf{\Gamma}^{1-\alpha} \mathbf{V}^\top\), si se elige \(\alpha=1\), se obtiene el biplot de covarianza, donde la magnitud de cada rayo es proporcional a la desviación estándar de las coordenadas clr. Si, en cambio, se toma \(\alpha=0\), entonces se obtiene el biplot de forma, que ayuda a identificar qué tan bien representadas se encuentran las componentes clr.
Véase

\[\begin{aligned} \hat {\mathbf{Z}}_2 = & \begin{pmatrix} \gamma_1^\alpha \mathbf{u}_1 & \gamma_2^\alpha \mathbf{u}_2 \end{pmatrix} \begin{pmatrix} \gamma_1^{1-\alpha} \mathbf{v}^\top_1 \\ \gamma_2^{1-\alpha} \mathbf{v}^\top_2 \end{pmatrix} \\ =& \begin{pmatrix} \sqrt{n}\gamma_1^{1-\alpha}u_{11} & \sqrt{n}\gamma_2^{1-\alpha}u_{21} \\ \sqrt{n}\gamma_1^{1-\alpha}u_{12} & \sqrt{n}\gamma_2^{1-\alpha}u_{22} \\ \vdots & \vdots \\ \sqrt{n}\gamma_1^{1-\alpha}u_{1n} & \sqrt{n}\gamma_2^{1-\alpha}u_{2n} \end{pmatrix} \begin{pmatrix} \dfrac{\gamma_1^\alpha v_{11}}{\sqrt{n}} & \dfrac{\gamma_1^\alpha v_{21}}{\sqrt{n}} & \cdots & \dfrac{\gamma_1^\alpha v_{D1}}{\sqrt{n}} \\ \dfrac{\gamma_2^\alpha v_{12}}{\sqrt{n}} & \dfrac{\gamma_2^\alpha v_{22}}{\sqrt{n}} & \cdots & \dfrac{\gamma_2^\alpha v_{D2}}{\sqrt{n}} \end{pmatrix} \\ =& \begin{pmatrix} \mathbf{a}_1 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{a}_n \end{pmatrix} \begin{pmatrix} \mathbf{b}_1 & \mathbf{b}_2 & \cdots & \mathbf{b}_D \end{pmatrix} \end{aligned}\]

Donde los vectores \(\mathbf{a}_i \in A\), con \(i= 1,\dots, n\) y los vectores \(\mathbf{b}_j \in B\), con \(j=1,\dots, D\) , donde \(A\) es una matriz \((n,2)\) y \(B\) es una matriz \((D,2)\).

  • Los vectores \(\mathbf{a}_i\) se denominan marcadores de fila de \(\hat {\mathbf{Z}}_2\) que corresponden a las proyecciones de las \(n\) muestras en el plano definido por los dos primeros autovectores de \(\mathbf{Z}\mathbf{Z}^\top\).

  • Los vectores \(\mathbf{b}_j\) son los marcadores de columna de \(\hat {\mathbf{Z}}_2\) que corresponden a las proyecciones del los \(D\) coeficientes-clr en el plano definido por los dos primeros autovectores de \(\mathbf{Z}^\top \mathbf{Z}\).

Luego se sobreponen los dos planos para la visualización de la relación entre los puntos y los coeficientes-clr.

2.1 Interpretación del bidimensional composicional biplot

Las características del biplot composicional son:

  1. Origen: Representa el centro de los datos composicionales, \(O\).

  2. Vértice: Un vértice en cada \(\mathbf{b}_j\), para cada una de las \(D\) variables clr.

  3. Marcador de Observación: punto ubicado en \(\mathbf{a}_i\), asociado a cada una de las \(n\) observaciones.

  4. Rayo: La unión entre el origen \(O\) y el vértice \(\mathbf{b}_j\), \(\overline{O\mathbf{b}_j}\). \[\lvert \overline{O \mathbf{b}_j} \rvert^2 \approx \mathrm{var}\left(\ln{\left( \frac{X_j}{g_m({X})} \right)}\right)\]

  5. Enlace: La unión de dos vértices \(\mathbf{b}_j\) y \(\mathbf{b}_k\), \(\overline{\mathbf{b}_j\mathbf{b}_k}\). \[\lvert \overline{\mathbf{b}_j\mathbf{b}_k} \rvert^2 \approx \mathrm{var}\left(\ln{\left( \frac{X_j}{X_k} \right)}\right)\]

A partir de la geometría de los vértices, rayos y enlaces pueden estudiarse las principales relaciones de variabilidad entre las partes:

  1. Longitud de los enlaces: Un enlace corto entre \(\mathbf{b}_j\) y \(\mathbf{b}_k\) indica que \[\mathrm{var}\left(\ln{\left( \frac{X_j}{X_k} \right)}\right) \approx 0\] y, entonces, \(X_j\) y \(X_k\) son proporcionales. Por el contrario, un enlace largo indica una alta variabilidad de dicho log-ratio y, por tanto, una mayor contribución de dicha relación a la variabilidad composicional.

  2. Ángulo entre enlaces: Proporcionan información sobre la correlación entre subcomposiciones. Si dos enlaces, \(\overline{\mathbf{b}_j\mathbf{b}_k}\) y \(\overline{\mathbf{b}_i\mathbf{b}_l}\), se intersectan en \(M\), \[\cos{ \left(\overline{\mathbf{b}_j\mathbf{b}_k} M \overline{\mathbf{b}_i\mathbf{b}_l} \right)} \approx \mathrm{Corr}\left( \ln{\left( \frac{X_j}{X_k} \right)}, \ln{\left( \frac{X_i}{X_l} \right)} \right)\]

    1. Enlaces aproximadamente ortogonales sugieren una correlación cercana a cero entre los log-ratios correspondientes.

    2. Si varios vértices se encuentran aproximadamente sobre una misma línea, esto puede indicar que la subcomposición formada por esas partes presenta una estructura de variación unidimensional (subcomposicones colineales).

    3. Dos conjuntos de subcomposicones colineales cuyos rayos forman un ángulo de 90°, posiblemente son ortogonales.

  3. Análisis Subcomposicional: Los cocientes entre componentes se conservan al formar subcomposiciones. Por eso, el biplot de una subcomposición puede obtenerse seleccionando los vértices correspondientes a sus partes y tomando como centro el centroide de esos vértices.

  4. Marcadores de Observación: Estos se pueden interpretar proyectando cada observación sobre el eje formado por los rayos o los enlaces:

    1. La dirección del eje corresponde al aumento de la cantidad representada, mientras que la dirección opuesta corresponde a su disminución. Si el marcador se proyecta en la dirección de la parte situada en el numerador, el log-ratio correspondiente es mayor que su promedio. Si se proyecta en la dirección contraria, el log-ratio es menor que su promedio.

    2. Si se proyecta una observación y el centro sobre el eje formado por el enlace \(\overline{\mathbf{b}_i\mathbf{b}_j}\), y sus proyecciones tienen una distancia pequeña, entonces esa observación tiene un log-ratio \(\ln{(X_i/X_j)}\) muy similar al log-ratio de todos los datos.

    3. Si la distancia entre la proyección de una observación y la proyección del centro sobre el enlace \(\overline{\mathbf{b}_i\mathbf{b}_j}\) equivale a una longitud del enlace, entonces \(\ln\left(\frac{X_i}{X_j}\right)\) se encuentra aproximadamente una desviación estándar por encima de su valor promedio.

Figure 1: Reglas de interpretación de los marcadores en un biplot composicional
Figure 2: Descripción de los elementos del biplot composicional

3 Ejemplo

Code 
library(compositions)

# Generacion de la muestra de 100 atletas
set.seed(123)
n_atletas = 100

prot_raw = rlnorm(n_atletas, meanlog = log(0.25), sdlog = 0.3)
carb_raw = rlnorm(n_atletas, meanlog = log(0.55), sdlog = 0.3)
gras_raw = rlnorm(n_atletas, meanlog = log(0.20), sdlog = 0.3)

# Construccion de la matriz
X_mat = data.frame(Proteinas = prot_raw,
                   Carbohidratos = carb_raw,
                   Grasas = gras_raw)

X_acomp = acomp(X_mat)

head(X_acomp)
     Proteinas   Carbohidratos Grasas     
[1,] "0.2026821" "0.4262887"   "0.3710293"
[2,] "0.2076017" "0.5285815"   "0.2638168"
[3,] "0.3645849" "0.4666581"   "0.1687570"
[4,] "0.2588950" "0.5024342"   "0.2386707"
[5,] "0.3057795" "0.4864089"   "0.2078117"
[6,] "0.3687251" "0.4784161"   "0.1528587"
attr(,"class")
[1] "acomp"
Code 
# Diagrama ternario de las composiciones
plot(
  X_acomp,
  cex = 0.7,
  col = "blue",
  main = "Composición nutricional de los atletas"
)

Diagrama ternario de la composición nutricional de los atletas
Code 
# Calculo del centro composicional (g(X))
g = mean(X_acomp)

print(g)
    Proteinas Carbohidratos        Grasas 
  "0.2577053"   "0.5342630"   "0.2080318" 
attr(,"class")
[1] "acomp"
Code 
# Centrado de los datos (X - g(X))
X_center = X_acomp - g

head(X_center)
     Proteinas   Carbohidratos Grasas     
[1,] "0.2335240" "0.2369126"   "0.5295634"
[2,] "0.2629943" "0.3229949"   "0.4140107"
[3,] "0.4564540" "0.2818158"   "0.2617302"
[4,] "0.3248746" "0.3041161"   "0.3710093"
[5,] "0.3832617" "0.2940742"   "0.3226641"
[6,] "0.4674208" "0.2925361"   "0.2400431"
attr(,"class")
[1] "acomp"
Code 
# Calculo de clr(X)
clr_X = clr(X_acomp)
clr_X_center = clr(X_center)

head(clr_X)
       Proteinas Carbohidratos     Grasas
[1,] -0.44937344     0.2941046  0.1552688
[2,] -0.39140310     0.5431726 -0.1517695
[3,]  0.17448742     0.4213247 -0.5958121
[4,] -0.19390140     0.4691407 -0.2752393
[5,] -0.02598454     0.4382009 -0.4122164
[6,]  0.20670252     0.4671321 -0.6738346
attr(,"class")
[1] "rmult"
Code 
head(clr_X_center)
       Proteinas Carbohidratos      Grasas
[1,] -0.27772486   -0.26331832  0.54104318
[2,] -0.21975452   -0.01425035  0.23400486
[3,]  0.34613601   -0.13609824 -0.21003777
[4,] -0.02225282   -0.08828230  0.11053512
[5,]  0.14566405   -0.11922203 -0.02644202
[6,]  0.37835111   -0.09029089 -0.28806022
attr(,"class")
[1] "rmult"
Code 
# Factorización SVD
SVD_clr_X = svd(clr_X)
SVD_clr_X_center = svd(clr_X_center)


# PCA - sin centrar
pca_X = prcomp(clr_X, center = FALSE, scale. = FALSE)

dots = rep(".",times=nrow(X_acomp))

biplot(
  pca_X,
  scale = 0,
  main = "Biplot X sin centrar",
  xlabs = dots
)


# PCA - centrado
pca_X_center = prcomp(clr_X_center, center = FALSE, scale. = FALSE)

dots = rep(".",times=nrow(clr_X_center))

biplot(
  pca_X_center,
  scale = 0,
  main = "Biplot X centrado",
  xlabs = dots
)

4 Aplicación

Code 
GeoChemSed=read.csv("geochemsed.csv",header=TRUE,skip=1)

x=acomp(GeoChemSed,4:13)

pcx=princomp(x)

summary(pcx)
Importance of components:
                          Comp.1    Comp.2     Comp.3     Comp.4     Comp.5
Standard deviation     1.0474538 0.5416805 0.25687273 0.20221390 0.16660853
Proportion of Variance 0.7121545 0.1904544 0.04282926 0.02654157 0.01801769
Cumulative Proportion  0.7121545 0.9026089 0.94543820 0.97197976 0.98999745
                            Comp.6      Comp.7     Comp.8       Comp.9
Standard deviation     0.101967053 0.051561497 0.03956718 0.0280837332
Proportion of Variance 0.006748764 0.001725661 0.00101619 0.0005119343
Cumulative Proportion  0.996746215 0.998471876 0.99948807 1.0000000000

Para visualizar el porcentaje de varianza que capturará el biplot:

Code 
sum(pcx$sdev[1:2]^2)/mvar(x)
[1] 0.9026089

Figure 3: Tres diagramas ternarios donde el biplot sugiere un patrón aproximadamente unidimensional

5 Referencias

Aitchison, John, and Michael Greenacre. 2002. “Biplots of Compositional Data.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 51 (4): 375–92. https://doi.org/10.1111/1467-9876.00275.
Boogaart, K. Gerald van den, and Raimon Tolosana-Delgado. 2013. Analyzing Compositional Data with r. Springer. https://doi.org/10.1007/978-3-642-36809-7.
Cortés-Rodriguéz, María. 2023. “BIPLOT DE DATOS COMPOSICIONALES, UNA NUEVA HERRAMIENTA ESTADÍSTICA PARA EL ESTUDIO DE TEST PSICOLÓGICOS: APLICACIÓN AL CUESTIONARIO DE BIENESTAR DE CAROL RYFF.” PhD thesis, Universidad de Salamanca.
Filzmoser, Peter, Karel Hron, and Matthias Templ. 2018. Applied Compositional Data Analysis: With Worked Examples in r. Springer. https://doi.org/10.1007/978-3-319-96422-5.
Greenacre, Michael. 2010. “The Standard Biplot.” Statistics and Computing 20: 341–54. https://doi.org/10.1007/s11222-009-9143-x.
Pawlowsky-Glahn, Vera, Juan Egozcue, and Raimon Tolosana-Delgado. 2015. Modeling and Analysis of Compositional Data. Wiley.