1 Librerías

library(readr)    #Lectura de archivos .csv
library(dplyr)    #Filtros
library(sm)       #Comparar densidades
library(ggplot2)  #Graficos
library(GGally)   #Diagramas de dispersion
library(yacca)    #Correlacion canonica
library(CCP)      #Pruebas de hipotesis (correlacion canonica)
library(reshape2) #Mapas de calor
library(CCA)      #Biplot

2 Lectura de datos

Los datos fueron extraídos de https://www.kaggle.com/datasets/sorelyss/icfes-colombia-20182021 y contienen información acerca de los puntajes en los distintos componentes que se consideran en la prueba de estado Icfes Saber 11 que se realiza en Colombia a los estudiantes que están en último año de bachillerato.

Para el análisis solo se considerarán los colegios de carácter oficial ubicados en la ciudad de Bogotá D.C. y los resultados del periodo 2018-II.

datos <- read.csv("datos_icfes.csv")[,2:6]
names(datos) <- c("Matemáticas","C_Naturales","L_Crítica","Inglés","C_Ciudadanas")
str(datos)
## 'data.frame':    48383 obs. of  5 variables:
##  $ Matemáticas : int  28 52 58 54 80 51 46 30 58 54 ...
##  $ C_Naturales : int  29 49 50 50 76 49 46 27 60 56 ...
##  $ L_Crítica   : int  39 57 63 60 70 54 43 35 61 61 ...
##  $ Inglés      : int  37 59 49 53 71 55 50 35 45 48 ...
##  $ C_Ciudadanas: int  29 46 55 47 76 49 47 26 58 47 ...
head(datos)
##   Matemáticas C_Naturales L_Crítica Inglés C_Ciudadanas
## 1          28          29        39     37           29
## 2          52          49        57     59           46
## 3          58          50        63     49           55
## 4          54          50        60     53           47
## 5          80          76        70     71           76
## 6          51          49        54     55           49

3 Grupos de variables

3.1 Grupo X: Ciencias

#Grupo X: Ciencias
X <- datos[,c(1,2)]
head(X)
##   Matemáticas C_Naturales
## 1          28          29
## 2          52          49
## 3          58          50
## 4          54          50
## 5          80          76
## 6          51          49

3.2 Grupo Y: Humanidades

#Grupo Y: Humanidades
Y <- datos[,c(3,5,4)]
head(Y)
##   L_Crítica C_Ciudadanas Inglés
## 1        39           29     37
## 2        57           46     59
## 3        63           55     49
## 4        60           47     53
## 5        70           76     71
## 6        54           49     55

4 Análisis descriptivo

Media Mediana Mínimo Máximo Desv. Est.
Matemáticas 51.985 52 16 100 10.013
C_Naturales 51.198 51 0 100 8.881
L_Crítica 54.005 54 0 100 8.910
Inglés 52.134 52 0 100 9.483
C_Ciudadanas 49.857 50 18 100 10.555

5 Matriz de correlación

round(cor(datos),3) #Matriz de correlacion completa
##              Matemáticas C_Naturales L_Crítica Inglés C_Ciudadanas
## Matemáticas        1.000       0.755     0.673  0.575        0.678
## C_Naturales        0.755       1.000     0.705  0.626        0.759
## L_Crítica          0.673       0.705     1.000  0.584        0.763
## Inglés             0.575       0.626     0.584  1.000        0.604
## C_Ciudadanas       0.678       0.759     0.763  0.604        1.000

6 Variables canónicas

6.1 Objeto yacca

# Objeto del analisis de correlacion canonico - Libreria "yacca"
acc <- cca(X, Y, standardize.scores = T)
# Matrices de covarianzas
Sx  <- cov(X)
Sy  <- cov(Y)
Sxy <- cov(X,Y)
Syx <- t(Sxy)

## U ---------------------------------------------------------------------------

E1 <- solve(Sx)%*%Sxy%*%solve(Sy)%*%Syx

#Vectores propios
vect_prop_u <- eigen(E1)$vectors

#Variables canonicas U* (coeficientes sin estandarizar)
u1 <- vect_prop_u[1,1]*X[,1] + vect_prop_u[2,1]*X[,2]
u2 <- vect_prop_u[1,2]*X[,1] + vect_prop_u[2,2]*X[,2]
  
# Estandarizar coeficientes
u11 <- vect_prop_u[1,1]/sd(u1)
u12 <- vect_prop_u[2,1]/sd(u1)

u21 <- vect_prop_u[1,2]/sd(u2)
u22 <- vect_prop_u[2,2]/sd(u2)

#Variables canonicas U (coeficientes estandarizados)
u1 <- u11*X[,1] + u12*X[,2]
u2 <- u21*X[,1] + u22*X[,2]

#Restriccion de varianza 1
var(u1)
## [1] 1
var(u2)
## [1] 1
#Restriccion cor(U1,U2) = 0
round(cor(u1,u2),4)
## [1] 0

6.2 Comparación con librería yacca

#Comparacion con libreria "yacca"
tab_u <- cbind(c(u11,u12,u21,u22), c(acc$xcoef))
colnames(tab_u) <- c("Manual", "yacca")
rownames(tab_u) <- c("u11", "u12", "u21", "u22")
knitr::kable(round(tab_u,3), align = "c")
Manual yacca
u11 -0.036 0.036
u12 -0.078 0.078
u21 -0.148 0.148
u22 0.153 -0.153 
#Variables canonicas U (coeficientes estandarizados)
u1 <- acc$xcoef[1,1]*X[,1] + acc$xcoef[2,1]*X[,2]
u2 <- acc$xcoef[1,2]*X[,1] + acc$xcoef[2,2]*X[,2]
## V --------------------------------------------------------------------------

#Variables canonicas V (libreria yacca)
v1 <- acc$ycoef[1,1]*Y[,1] + acc$ycoef[2,1]*Y[,2] + acc$ycoef[3,1]*Y[,3]
v2 <- acc$ycoef[1,2]*Y[,1] + acc$ycoef[2,2]*Y[,2] + acc$ycoef[3,2]*Y[,3]

#Restriccion de varianza 1, cor(V1,V2) = 0, cor(V1,U2) = 0 y cor(V2,U1) = 0
var(v1)
## [1] 1
var(v2)
## [1] 1
round(cor(v1,v2),4)
## [1] 0
round(cor(v1,u2),4)
## [1] 0
round(cor(v2,u1),4)
## [1] 0

6.2.1 Correlaciones canónicas

tab_rho <- cbind(sqrt(eigen(E1)$values), acc$corr)
colnames(tab_rho) <- c("Manual", "yacca")
rownames(tab_rho) <- c("Corr_canonica_1", "Corr_canonica_2")
knitr::kable(round(tab_rho,3), align = "c")
Manual yacca
Corr_canonica_1 0.825 0.825
Corr_canonica_2 0.096 0.096 
# cor(U,V)
matrix(c(round(c(abs(cor(u1,v1)),acc$corr[1]),3),
         round(c(abs(cor(u2,v2)),acc$corr[2]),3)), ncol = 2, byrow = T)
##       [,1]  [,2]
## [1,] 0.825 0.825
## [2,] 0.096 0.096

6.3 Coeficientes canónicos

# Coeficientes canonicos estandarizados

#U
round(acc$xcoef,3)
##              CV 1   CV 2
## Matemáticas 0.036  0.148
## C_Naturales 0.078 -0.153

\[U_1 = 0.036*\textit{Matemáticas} + 0.078*\textit{C_Naturales}\]

\[U_2 = 0.148*\textit{Matemáticas} - 0.153*\textit{C_Naturales}\]

#V
round(acc$ycoef,3)
##               CV 1   CV 2
## L_Crítica    0.038  0.167
## C_Ciudadanas 0.049 -0.134
## Inglés       0.029  0.001

\[V_1 = 0.038*\textit{L_Crítica} + 0.049*\textit{C_Ciudadanas} + 0.029*\textit{Inglés}\]

\[V_2 = 0.167*\textit{L_Crítica} - 0.134*\textit{C_Ciudadanas} + 0.001*\textit{Inglés}\]

7 Prueba de hipótesis

#Test asintotico
p.asym(rho = acc$corr, N = dim(datos)[1], p = dim(X)[2], q = dim(Y)[2], 
       tstat = "Wilks")
## Wilks' Lambda, using F-approximation (Rao's F):
##               stat     approx df1   df2 p.value
## 1 to 2:  0.3161451 12554.2828   6 96756       0
## 2 to 2:  0.9907945   224.7446   2 48379       0
Wilks’ Lambda, using F-approximation (Rao’s F)
Statistic Stat Approx df1 df2 p.value
1 to 2 0.3161451 12554.2828 6 96756 0
2 to 2 0.9907945 224.7446 2 48379 0

8 Gráfico de variables canónicas

9 Índice de redundancia (según artículo [1])

#Y
#Porcentaje de la variabilidad de Y explicada por V1
#mean(cor(v1,Y)^2)

#Porcentaje de la variabilidad de Y explicada por V2
#mean(cor(v2,Y)^2)

#X
#Porcentaje de la variabilidad de X explicada por U1
#mean(cor(u1,X)^2)

#Porcentaje de la variabilidad de X explicada por U2
#mean(cor(u2,X)^2)

#Tabla de indice de redundancia
tab_IR <- matrix(c(mean(cor(u1,X)^2), mean(cor(v1,Y)^2),
                   mean(cor(u2,X)^2), mean(cor(v2,Y)^2)), ncol = 2, byrow = T)

colnames(tab_IR) <- c("IR_X", "IR_Y")
rownames(tab_IR) <- c("V_Canónica_1","V_Canónica_2")
knitr::kable(tab_IR, digits = 3, align = "c")
IR_X IR_Y
V_Canónica_1 0.867 0.763
V_Canónica_2 0.133 0.082

10 Cargas canónicas

# Cargas canonicas para el grupo X (Ciencias)
# Correlaciones de las X con cada varaible canonica U
knitr::kable(round(acc["xstructcorr"]$xstructcorr,3), align = "c")
CV 1 CV 2
Matemáticas 0.890 0.456
C_Naturales 0.971 -0.239 
# Cargas canonicas para el grupo Y (Humanidades)
# Correlaciones de las Y con cada varaible canonica V
knitr::kable(round(acc["ystructcorr"]$ystructcorr,3), align = "c")
CV 1 CV 2
L_Crítica 0.891 0.410
C_Ciudadanas 0.940 -0.276
Inglés 0.782 0.020

11 Cargas cruzadas

# Cargas cruzadas para el grupo X (Ciencias)
#(Correlaciones entre las X's y las variables canonicas del grupo Y)

knitr::kable(round(acc["xcrosscorr"]$xcrosscorr,3), align = "c")
CV 1 CV 2
Matemáticas 0.734 0.044
C_Naturales 0.801 -0.023 
# Cargas cruzadas para el grupo Y (Humanidades)
#(Correlaciones entre las Y's y las variables canonicas del grupo X)

knitr::kable(round(acc["ycrosscorr"]$ycrosscorr,3), align = "c")
CV 1 CV 2
L_Crítica 0.735 0.039
C_Ciudadanas 0.776 -0.027
Inglés 0.645 0.002

12 Proporción de varianza explicada

#Fraccion de la varianza del grupo X asociada con cada variable canonica U
p1 <- acc["xcanvad"]$xcanvad

#Fraccion de la varianza del grupo X asociada con cada variable canonica V
p2 <- acc["xvrd"]$xvrd

#Fraccion de la varianza del grupo Y asociada con cada variable canonica V
p3 <- acc["ycanvad"]$ycanvad

#Fraccion de la varianza del grupo Y asociada con cada variable canonica U
p4 <- acc["yvrd"]$yvrd

#Tabla de proporciones de varianza
tab_prop <- round(cbind(p1,p2,p3,p4), 3)
colnames(tab_prop) <- c("G1 por G1", "G1 por G2", "G2 por G2", "G2 por G1")
rownames(tab_prop) <- c("Variable canónica 1", "Variable canónica 2")
knitr::kable(tab_prop, align = "c")
G1 por G1 G1 por G2 G2 por G2 G2 por G1
Variable canónica 1 0.867 0.591 0.763 0.520
Variable canónica 2 0.133 0.001 0.082 0.001

13 Biplot (plano canónico)

cc <- cc(X,Y) 
plt.cc(cc, var.label = TRUE, d1 = 1, d2 = 2, type = "b")

14 Referencias

[1] Akour, I.; AL Rahamneh, A.; Al Kurdi, B.; and Alhamad, A. (2023).
“Using the Canonical Correlation Analysis Method to Study Students Levels in Face-to-Face and Online Education in Jordan”.
Information Sciences Letters: Vol. 12 : Iss. 2, PP -.
Disponible en: https://digitalcommons.aaru.edu.jo/isl/vol12/iss2/53

[2] Díaz Monroy, L y Morales Rivera, M. (2012). Análisis estadístico de datos multivariados. Universidad Nacional de Colombia.

[3] Rencher, A. C. (1998). Multivariate Statistical Inference and Applications, John Wiley and Sons, New York.