8 Análisis de atributos
Instalación/carga librerías/datos utilizados
Fecha última actualización: 2025-01-07
if (!require(highcharter)) install.packages('highcharter')
library(highcharter)
if (!require(openxlsx)) install.packages('openxlsx')
library(openxlsx)
if (!require(plotly)) install.packages('plotly')
library(plotly)
if (!require(ggplot2)) install.packages('ggplot2')
library(ggplot2)
if (!require(tidyverse)) install.packages('tidyverse')
library(tidyverse)
8.1 Introducción
De manera general, podemos considerar que el análisis de atributos consiste en estudiar, comparar y manipular las variables de nuestro análisis de datos, transformándolas o combinándolas para crear nuevas variables (o eliminar algunas) con objeto de conseguir una nueva colección de variables que representen mejor, y de manera más compacta, la información contenida globalmente en la colección de variables inicial. En lo que sigue, para desarrollar los conceptos de esta tema, consideraremos que nuestras variables son una colección de atributos, como población, superficie, PIB, etc.., asociados a los países del mundo, tal y como aparece en la tabla owid_country
. Cada atributo (o variable) la denotaremos por el vector xj=(x1j,x2j,..,xmj) donde m representa el número de países. Es decir, que si el primer atributo (j=1) es la población, entonces, xi1 representaría la población del país i.
La primera fase de nuestro análisis de atributos consiste en comparar cada par de atributos xj y xj′ entre sí. La manera más fácil de comparar variables es calcular su correlación, r, que viene dada por la fórmula:
r=∑mi=1(xij−ˉxj)(xij′−ˉxj′)√∑mi=1(xij−ˉxj)2∑mi=1(xij′−ˉxj′)2
Donde, ˉxj y ˉxj′ representan la media de las respectivas variables. El valor de la correlación r se mueve en el rango entre -1 y 1 y se puede interpretar de la siguiente manera:
- Si r=1 la variable xj′ se puede obtener de forma exacta y lineal a partir de xj, si el valor de una variable aumenta también aumenta la otra, y si dibujamos la gráfica de dispersión de ambas variables, todos los puntos caen de forma exacta sobre la recta de regresión
- Si r=0 no existe relación lineal entre ambas variables y la gráfica de dispersión, normalmente, nos dará una nube de puntos dispersa
- Si r=−1, el caso es idéntico a r=1 con la única diferencia de que, en este caso, cuando una variable aumenta, la otra disminuye.
En general, cuanto más cerca de 1 o -1 está la correlación, mejor es la relación lineal entre las variables y mejor alineados estarán sus valores sobre la recta de regresión. Una propiedad importante de la correlación es que es independiente de la unidad de medida de las variables (siempre que la transformación entre las medidas sea lineal), es decir, da igual medir una población en millones o en miles, o una temperatura en grados Celsius o Fahrenheit, la correlación de estas variables con otras variables da el mismo valor independientemente de la unidad de medida utilizada.
Veamos un ejemplo comparando el total de delitos cometidos y el total de consumo de pan y cereales obtenidos, por comunidades autónomas, a través de consultas al INE. Para calcular la correlación usamos la función cor
.
INE.CONSUMO.PAN.DELITOS <- read_csv("https://ctim.es/AEDV/data/INE.CONSUMO.PAN.DELITOS.csv") %>%
as_tibble()
INE.CONSUMO.PAN.DELITOS
## # A tibble: 19 × 3
## `Comunidades y Ciudades Autónomas` total.consumo.pan.cereales total.delitos
## <chr> <dbl> <dbl>
## 1 01 Andalucía 6957. 68853.
## 2 02 Aragón 1117. 8458.
## 3 03 Asturias, Principado de 927. 7825.
## 4 04 Balears, Illes 944. 10590.
## 5 05 Canarias 1805. 19146.
## 6 06 Cantabria 500. 4392.
## 7 07 Castilla y León 2130. 14720.
## 8 08 Castilla - La Mancha 1718. 11670.
## 9 09 Cataluña 6140. 58416.
## 10 10 Comunitat Valenciana 4167. 44809.
## 11 11 Extremadura 926. 7206.
## 12 12 Galicia 2389. 17904.
## 13 13 Madrid, Comunidad de 5289. 45146.
## 14 14 Murcia, Región de 1198. 12126.
## 15 15 Navarra, Comunidad Foral de 528. 4244.
## 16 16 País Vasco 1834. 15788.
## 17 17 Rioja, La 263 2372.
## 18 18 Ceuta 65.6 1759.
## 19 19 Melilla 62.8 1458.
paste("correlación =",
cor(
INE.CONSUMO.PAN.DELITOS$total.consumo.pan.cereales,
INE.CONSUMO.PAN.DELITOS$total.delitos)
)
## [1] "correlación = 0.98987954033385"
Observamos una correlación alta entre ambas variables. Vamos a ilustrar la relación entre las variables con un diagrama de dispersión:
p <- INE.CONSUMO.PAN.DELITOS %>%
ggplot(aes(total.consumo.pan.cereales,total.delitos,label=`Comunidades y Ciudades Autónomas`)) +
geom_point() +
geom_smooth(method = lm, se = FALSE)
ggplotly(p)
Observamos que como la correlación es alta, los puntos se encuentran próximos a la recta de regresión. Este ejemplo nos ilustra además que una alta regresión entre variables no implica que exista una causalidad entre la variables observadas, es decir, el consumo de pan y cereales no es un factor que provoque delincuencia. La correlación alta es debida a que ambas variables crecen con la población y la relación entre ellas es indirecta dado que las comunidades con más población consumen más pan y a la vez se producen más delitos. Establecer relaciones de causalidad, exclusivamente a partir de la correlación entre variables, es un error habitual que el buen científico de datos debe evitar.
Puede ser que la correlación entre las variables mejore si le aplicamos previamente a alguna de ellas (o a las dos) una transformación no -lineal (como un logaritmo) antes de calcular la correlación. Esto es lo que hacemos en el cuadro de mandos de ejemplo Dashboard3.Rmd
del tema anterior, en él, se estudia como transformar los atributos, usando la función logaritmo o de Box-Cox para estudiar la relación entre pares de atributos. Por ejemplo observamos que si aplicamos a PIB per capita un logaritmo, el resultado de su correlación con la esperanza de vida mejora mucho. Esto nos indica que la relación entre las dos variables es no-lineal.
Diremos que una variable xj no aporta información nueva a nuestro análisis de datos si se puede obtener a partir de las otras variables. Es decir si
xij=F({xik}k≠j)∀i
para una cierta función F. Por tanto, podemos quitar xj de nuestro análisis de datos sin perder información. En particular, si la correlación entre dos variables está cerca de 1 o -1, podemos quitar una de ellas porque una se puede obtener a partir de la otra. Por ejemplo, en la tabla owid_country
se guardan como indicadores el porcentaje de personas mayores de 65 años y el de personas mayores de 70 años. Estos dos indicadores están muy relacionados entre sí. Si usamos el tablero de mando Dashboard3.Rmd
, observamos que el factor de correlación es de 0.99, y por tanto, se puede calcular con precisión una en función de la otra a partir de la recta de regresión, que en este caso viene dada por
aged_70_older=0.68⋅aged_65_older−0.39 Por tanto, podríamos eliminar el indicador de porcentaje de mayores de 70 años de nuestro análisis de datos sin perder información relevante. Otro ejemplo de variable que podríamos desechar son las variables que poseen un valor constante (en este caso la función F sería una constante). Este tipo de variables se consideran desechables porque no aportan nada a la hora de discriminar entre unos países y otros.
En este tema vamos a estudiar como combinar las variables y reducir su número sin perder información relevante. A este problema se le denomina reducción de la dimensionalidad.
Una limitación de las técnicas que vamos a analizar es que solo usaremos combinaciones lineales de las variables, es decir sumas de variables multiplicadas por unos coeficientes, eso deja fuera dependencias más complejas como las que vienen dadas por la función logaritmo o la transformación de Box-Cox utilizadas en el cuadro de mandos Dashboard3.Rmd
. Una forma de añadir transformadas no-lineales sería, antes de empezar con nuestro análisis lineal, añadir/sustituir en la tabla, variables que sean el resultado de aplicar transformaciones no-lineales a las variables iniciales. Por ejemplo, a la vista de los resultados del cuadro de mandos Dashboard3.Rmd
, podríamos sustituir la variable gdp_per_capit
por log(gdp_per_capit)
que posee una mejor relación lineal con las otras que la variable inicial.
8.2 La matriz de correlación
Tal y como vimos en la sección anterior, una manera sencilla de eliminar variables es estudiar la correlación entre pares de ellas y si la correlación es muy cercana a 1 o -1, podemos eliminar una de ellas. Para ello se utiliza la matriz de correlación donde cada valor de la matriz corresponde al valor de la correlación, r, entre 2 variables. Por ejemplo en la siguiente gráfica interactiva visualizamos la matriz de correlación de la tabla owid_country
. Previamente hay que quitar las variables no numéricas, puesto que la correlación solo se calcula para variables numéricas y además hay que eliminar los valores NA
antes de hacer el cálculo.
owid_country %>%
select(-location,-continent,-iso_code) %>% # eliminamos de la tabla las variables no-numéricas
cor(use='complete.obs') %>% # cálculo matriz correlación eliminando NA previamente
hchart() # dibujo interactivo matriz de correlación
Figure 8.1: Matriz de correlación usando hchart
De la exploración de los valores de correlación, observamos, que aged_65_older
y aged_70_older
están muy relacionadas, que population
y cardiovasc_death_rat
están muy poco relacionadas con el resto de variables y que extreme_poverty
, tiene, en general, correlación negativa con el resto de variables, es decir que cuando esta variable crece, las otras decrecen.
8.3 Análisis de componentes principales (ACP)
El análisis de componentes principales, ACP, (PCA en inglés) transforma las variables xj en otras nuevas, yj, haciendo combinaciones lineales de las iniciales. El objetivo del ACP es buscar las combinaciones lineales de xj que maximizan la varianza de yj en orden descendente, es decir, y1 es la combinación lineal con mayor varianza posible e yn, es la combinación lineal con menor varianza posible. Cuando la varianza de yj es muy pequeña consideramos que no aporta información relevante y la desechamos, de esta forma reducimos la dimensionalidad quedándonos solo con las variables yj con varianza significativa, que son las que denominamos componentes principales.
Fundamento teórico
Conceptos algebraícos necesarios
Definición: El producto escalar de dos vectores u y v de tamaño n se define como (u,v)=u1v1+⋯+unvn. Los vectores siempre los consideraremos en vertical, es decir podemos escribir (u,v)=uTv, donde uT es el vector traspuesto.
Definición: Una matriz cuadrada U es ortonormal (o unitaria) si se cumple UTU=Id, es decir que la traspuesta por ella misma es la matriz identidad. Por tanto, para una matriz ortonormal la inversa es la traspuesta.
Teorema: Los vectores columna de una matriz ortonormal U=(u1,⋯,un) (con uk=(u1k,…,unk)) cumplen que (uk,uk)=1 y si k≠k′ entonces (uk,uk′)=0 es decir los vectores son perpendiculares.
Definición: el número real λk es un autovalor de la matriz cuadrada A, si existe un vector uk=(u1k,…,unk), no nulo, denominado autovector, tal que Auk=λkuk. Si uk es autovector, entonces, al multiplicarlo por cualquier número diferente de cero, sigue siendo autovector. Por ello siempre podemos ajustar el autovector para que cumpla (uk,uk)=1.
Teorema: Si la matriz cuadrada A de dimensión n es simétrica, entonces posee n autovalores que podemos ordenar en forma decreciente, λ1≥λ2≥⋯≥λn y además posee una matriz ortonormal de autovectores U=(u1,u2,..,un) formada por los autovectores uk en columnas.
Teorema: Si X es una matriz cualquiera de tamaño mxn, entonces A=XTX es cuadrada de tamaño nxn, simétrica y todos sus autovalores son mayores o iguales que cero.
Definición: Dados dos variables xj y xj′, la covarianza Cov(xj,xj′) entre ellas se define como:
Cov(xj,xj′)=∑mi=1(xij−ˉxj)(xij′−ˉxj′)m−1 donde m es el número de observaciones.
Definición: Si X=(x1,x2,...,xn) es una matriz que representa n variables x1,x2,...,xn, la matriz de covarianza de X viene dada por la matriz Σ de dimensión n×n cuyos elementos son Σjj′=Cov(xj,xj′)
Definición: Tipificar una variable, xj, consiste en restarle su media y dividirla por su desviación típica. Las variables tipificadas tienen media cero y varianza y desviación típica igual a 1.
Teorema: Si X=(x1,x2,...,xn) es una matriz que representa n variables x1,x2,...,xn tipificadas, entonces la matriz de covarianza Σ cumple:
Σ=XTXm−1 además
λ1+λ2+..+λn=n donde λj son los autovalores de la matriz de covarianza Σ. Además, como la matriz de covarianza Σ es proporcional a XTX, sus autovectores son los mismos y los autovalores de XTX son los de Σ multiplicados por m−1.
Fundamento matemático del ACP
Consideremos n variables xj (los vectores columna de nuestra tabla), y las observaciones en forma de matriz X=xij de tamaño mxn. Vamos a sustituir X por una nueva colección de variables Y=XU donde U es una matriz ortonormal de tamaño n. Es decir, las nuevas variables yj=Xuj son combinaciones lineales de las anteriores obtenidas a partir de una matriz ortonormal que se pueden expresar como :
y1=u11x1+u21x2+⋯+un1xny2=u12x1+u22x2+⋯+un2xn⋯yn=u1nx1+u2nx2+⋯+unnxn
El análisis por componentes principales es una técnica para elegir la matriz ortonormal U de tal manera que las variables yj=Xuj (las componentes principales) estén ordenadas de forma descendente en función de su varianza, es decir, la primera componente principal y1=Xu1 corresponde a la combinación lineal de las variables originales que tiene la mayor varianza posible, y2=Xu2 sería la combinación ortogonal a la anterior ((u1,u2)=0) con la mayor varianza posible, y3=Xu3 sería la combinación ortogonal a las anteriores ((u1,u3)=0 y (u2,u3)=0) con la mayor varianza posible, y así sucesivamente. Cuanto mayor es la varianza de yj más información contiene la variable en términos, por ejemplo, de distinguir entre las diferentes observaciones. Nótese que una variable con valor constante (es decir varianza cero) no aporta ninguna información para discernir entre unas observaciones y otras. Dado que las variables xj pueden tener magnitudes y varianzas muy distintas, para equilibrar la aportación de cada variable, primero se tipifican, es decir se le resta su media y se divide por la desviación típica. Es decir, consideraremos que las variables xj tienen media cero y desviación típica 1. Si la media de cada variable xj es cero, la media de yj=Xuj que es una combinación lineal de xj también tendrá media cero y por tanto su varianza será
Var(yj)=∑mi=1y2ijm−1=yTjyjm−1=uTjXTXujm−1 la clave para maximizar la varianza son los autovectores de la matriz A=XTX. Si λj son los autovalores ordenados de A en forma descendente, es decir λ1≥λ2≥⋯≥λn y uj sus correspondientes autovectores, entonces
Var(yj)=uTjXTXujm−1=uTjλjujm−1=λjuTjujm−1=λjm−1 por tanto, el procedimiento para calcular las componentes principales de X se puede dividir en los siguientes pasos:
Tipificar la variables originales xj para que tenga media cero y varianza 1.
Calcular los autovalores λj (ordenados de forma descendente) y autovectores uj de la matriz XTX
Las componentes principales son yj=Xuj y su varianza es λjm−1
Ejemplo de ACP
Para ilustrar el proceso del cálculo de las componentes principales vamos a tomar un caso sencillo extraído de la tabla owid_country
donde usaremos 3 variables (x1=human_development_index
, x2=aged_65_older
y x3=aged_70_older
) y las observaciones de 5 países (China, Canada, Nigeria, Brazil y Germany). Para calcular el ACP la tabla solo puede tener variables numéricas. Para conservar el nombre de los países se usa la función column_to_rownames
que asigna los países como nombres de cada fila, de esta manera se gestionan como un elemento externo a la tabla y no interfieren con el cálculo del ACP.
data <- owid_country %>%
filter( # filtramos los 5 países que se usarán
location == "China" |
location == "Canada" |
location == "Nigeria" |
location == "Brazil" |
location == "Germany"
) %>%
column_to_rownames(var="location") %>%
select(human_development_index,aged_65_older,aged_70_older)
data
## human_development_index aged_65_older aged_70_older
## Brazil 0.765 8.552 5.060
## Canada 0.929 16.984 10.797
## China 0.761 10.641 5.929
## Germany 0.947 21.453 15.957
## Nigeria 0.539 2.751 1.447
A continuación imprimimos las variables tipificadas usando la función scale
y su matriz de correlación que muestra una alta correlación entre las variables .
## human_development_index aged_65_older aged_70_older
## Brazil -0.1409164 -0.4824390 -0.4932478
## Canada 0.8552170 0.6718445 0.5253853
## China -0.1652124 -0.1964691 -0.3389525
## Germany 0.9645488 1.2836202 1.4415691
## Nigeria -1.5136370 -1.2765565 -1.1347540
## attr(,"scaled:center")
## human_development_index aged_65_older aged_70_older
## 0.7882 12.0762 7.8380
## attr(,"scaled:scale")
## human_development_index aged_65_older aged_70_older
## 0.1646366 7.3049644 5.6320575
Por tanto la matriz X está definida por:
x1x2x3Brazil−0.14−0.48−0.49Canada0.860.670.53China−0.17−0.20−0.34Germany0.961.281.44Nigeria−1.51−1.28−1.13 Las nuevas variables yj que vamos a construir como combinación lineal de x1,x2 y x3 tienen la forma:
yj=u1j(−0.140.86−0.170.96−1.51)+u2j(−0.480.67−0.201.28−1.28)+u3j(−0.490.53−0.341.44−1.13)=Xuj
Como yj tiene siempre media cero, su varianza es yTjyj4=uTjXTXuj4. Si uj es un autovector normalizado de XTX para el autovalor λj, tendremos que la varianza de yj es uTjXTXuj4=λjuTjuj4=λj4
Calculamos ahora la matriz XTX: XTX=(4.003.853.683.854.003.963.683.964.00) A continuación, a mano, o usando cualquier software de cálculo matemático, calculamos la matriz ortonormal U con los autovectores de XTX
U=(0.570.780.260.58−0.16−0.790.58−0.610.55)
Simultáneamente se calculan los autovalores de XTX que son: λ1=11.66,λ2=0.33, y λ3=0.02 (observamos que la suma de los autovalores es aproximadamente 12=n⋅(m−1)=3⋅4) y por tanto las desviaciones típicas de yj son σj=√λj/4, es decir: σ1=1.71, σ2=0.29, σ3=0.07. Las componentes principales quedan entonces como:
y1=0.57⋅x1+0.58⋅x2+0.58⋅x3y2=0.78⋅x1−0.16⋅x2−0.60⋅x3y3=0.26⋅x1−0.79⋅x2+0.55⋅x3
y, por tanto, la nueva tabla de valores, Y=XU, da como resultado:
y1y2y3Brazil−0.650.270.08Canada1.180.24−0.02China−0.400.11−0.073Germany2.13−0.330.02Nigeria−2.27−0.28−0.00
Cálculo ACP con prcomp
La función prcomp
realiza de forma automática todo el proceso para el cálculo de las componentes principales. Vamos a utilizarla para el ejemplo anterior:
La función prcomp
devuelve una estructura (en este ejemplo la hemos llamado pca1
) que almacena, en el campo sdev
, la desviación típica de las componentes principales yj, en el campo rotation
la matriz de autovectores U y en el campo x
las componentes principales Y=XU. A continuación se imprimen estos resultados. Se puede observar que los resultados son los mismos que hemos calculado para este ejemplo en la sección anterior.
## [1] 1.70716366 0.28734312 0.05501071
## PC1 PC2 PC3
## human_development_index 0.5708408 0.7790411 0.2592986
## aged_65_older 0.5845568 -0.1638432 -0.7946375
## aged_70_older 0.5765709 -0.6051863 0.5489222
## PC1 PC2 PC3
## Brazil -0.6468462 0.2677714 0.0760700585
## Canada 1.1838460 0.2382161 -0.0237206224
## China -0.4045875 0.1086123 -0.0727761422
## Germany 2.1321196 -0.3313071 0.0214026069
## Nigeria -2.2645318 -0.2832927 -0.0009759008
Visualización ACP
El objetivo principal del ACP es reducir la dimensionalidad quedándonos con las primeras componentes principales que acumulan mayor varianza y desechando el resto. El criterio para decidir con cuantas componentes principales nos quedamos se basa justamente en el análisis de su varianza. Normalmente se analiza la varianza de las componentes principales yj en términos relativos, es decir el porcentaje que representa la varianza de yj respecto a la varianza global dada por la suma de todas las varianzas (∑jVar(yj)). Para visualizar esto, utilizamos un gráfico de barras interactivo, creado con plotly
, con el porcentaje de varianza explicada por cada componente principal.
p <- tibble(
label=paste("PC",1:length(pca1$sdev)), # creación etiquetas para el eje horizontal
varPercent = pca1$sdev^2/sum(pca1$sdev^2) * 100 # cálculo porcentaje de varianza explicada
) %>%
ggplot(aes(x=label,y=varPercent)) + # creación gráfico de barras interactivo
geom_bar(stat = "identity") +
labs(x= "Componentes Principales",
y= "Porcentaje varianza explicada")+
geom_text(aes(label=paste0(round(varPercent,digits = 1),"%")),
size=3,
colour = "blue",
nudge_y = 1.
)
ggplotly(p)
Figure 8.2: Porcentaje varianza explicada por cada componente principal usando ggplotly
Para decidir con cuantas componentes principales nos quedamos podemos, por ejemplo, quedarnos con las componentes principales que acumulen más de un 90% de la varianza explicada. Otro criterio, sería, por ejemplo, eliminar las componentes principales con la varianza explicada inferior al 2%. Para este ejemplo, con el primer criterio, solo nos quedaríamos con la primera componente principal, y con el segundo criterio con las dos primeras.
Un gráfico habitual para visualizar los resultados del ACP es hacer un gráfico de dispersión donde se pone el resultado de los valores de las dos primeras componentes principales para cada observación, es decir se representan los puntos (yi1,yi2) para i=1,⋯,m, de esta manera se aprecia como están distribuidos los valores de las dos primeras componentes. Además, para visualizar la contribución de cada variable xj a las dos primeras componentes, se dibuja, para cada j un segmento que va desde el (0,0) hasta el (uj1,uj2) para j=1,⋯,n. uj1 es el coeficiente de xj en la combinación lineal de la primera componente y1 y
uj2 es el coeficiente de xj en la combinación lineal de la segunda componente y2. Si este segmento se alinea con el eje x, significa que xj aporta poco a la segunda componente, y si además el segmento es de gran magnitud, significa que xj aporta a la primera componente más que el resto de variables. Para realizar este gráfico utilizaremos la función hchart
:
Figure 8.3: Gráfica ilustrativa de las dos primeras componentes principales usando hchart
A continuación vamos a reproducir el experimento, pero tomando todos los indicadores y para todos los países. Primero preparamos los datos y calculamos el ACP. Como el cálculo del ACP solo permite variables numéricas, tenemos que eliminar las variables no numéricas de la tabla. Para conservar los nombres de los países, en lugar de gestionarlo como variable, lo gestionamos como nombre de cada registro, de esta manera no interfiere en el cálculo del ACP
data <- owid_country %>%
na.omit() %>% # quitamos los registros con algún NA
column_to_rownames(var="location") %>% # asignamos el campo location como nombre de las filas
select(-continent,-iso_code) # eliminamos las variables no-numéricas
pca2 <- prcomp(data,scale = TRUE) # calculamos el ACP
Ahora realizamos la gráfica de barras con el porcentaje de varianza explicada por cada componente principal. Por defecto, la gráfica de barras ordena por orden alfabético las etiquetas del eje x, para que use el orden usado al crear las etiquetas usamos la función fct_inorder
que impide que se cambie el orden en el vector de etiquetas.
p <- tibble(
label=fct_inorder(paste("PC",1:length(pca2$sdev))),
varPercent = pca2$sdev^2/sum(pca2$sdev^2) * 100
) %>%
ggplot(aes(x=label,y=varPercent)) +
geom_col() +
labs(x= "Componentes Principales",y= "Porcentaje varianza explicada")+
geom_text(aes(label=paste0(round(varPercent,digits = 1),"%")), size=3, colour = "blue",nudge_y = 1.)
ggplotly(p)
Figure 8.4: Porcentaje varianza explicada por cada componente principal usando ggplotly
Por último presentamos la gráfica que ilustra el resultado para las dos primeras componentes (con el parámetro choices
seleccionamos las dos componentes que se van a representar).
Figure 8.5: Gráfica ilustrativa de las dos primeras componentes principales usando hchart
Referencias
[Ke23] Zoumana Keita. Principal Component Analysis in R Tutorial, 2023.
[Ku22] Joshua Kunst. Hchart Function, 2022.
[UGA] Data Analysis in the Geosciences (Principal Components Analysis), University of Georgia.