Clase AID: Clusters

cat("\014")



rm(list = ls())
library(ggplot2) 

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     

library(ggrepel)
library(devtools)

Loading required package: usethis

#install_github("vqv/ggbiplot")
library(ggbiplot)

Loading required package: plyr
Loading required package: scales
Loading required package: grid

library(readxl)
library(nortest)

situacion_analisis= c(1,2,3)      #grupo de alumnos a analizar según situacion academica
cuat_ini_analisis = 20090         #cuatrimestre inicial de analisis (siempre hasta el presente)
columnas = c(1,8,9,12:17)         #variables a tener en cuenta (menos que las de egresados)

datos_crudos=read_xlsx("base_de_alumnos.xlsx",
                       sheet = 1,col_types = "numeric")

datos_numericos_1 = datos_crudos[(datos_crudos$acuat_ing>=cuat_ini_analisis)&
                                  datos_crudos$sit_ac %in% situacion_analisis,]

datos_numericos_2 = data.frame(datos_numericos_1[columnas])
#se queda con casos completos
datos_completos = datos_numericos_2[complete.cases(datos_numericos_2),]
datos = datos_completos[c(3:ncol(datos_completos))]

#analsis exploratorio
par(mfcol = c(1,length(datos)))
    for (k in 1:length(datos)){
      boxplot(datos[k],main = names(datos[k]))
      grid()
    }

#analsis exploratorio
par(mfcol = c(2,4))
#par(mfcol = c(1,length(datos)))
    for (k in 1:length(datos)){
      hist(datos[,k],proba = T,main = names(datos[,k]),10)
      x0 <- seq(min(datos[, k]), max(datos[, k]), le = 50) 
      lines(x0, dnorm(x0, mean(datos[,k]), sd(datos[,k])), col = "red", lwd = 2) 
      grid()
    }

#analsis exploratorio
pval = list() 
par(mfcol = c(1,length(datos)))
    for (k in 1:length(datos)){
      qqnorm(datos[,k],main = names(datos[k]))
      qqline(datos[,k],col="red") 
      pval[k] = ad.test(datos[,k])$p.value
      grid()
    }

pval

[[1]]
[1] 3.7e-24

[[2]]
[1] 3.7e-24

[[3]]
[1] 3.7e-24

[[4]]
[1] 3.7e-24

[[5]]
[1] 3.7e-24

[[6]]
[1] 3.7e-24

[[7]]
[1] 3.7e-24

#estandarizo datos
datos_estandarizados = data.frame(scale(datos))
boxplot(datos_estandarizados)

#dispersograma - Ojo que tarda
pairs(datos_estandarizados,pch=19,cex=0.8)

#matriz de correlaciones aplicada directamente a "datos"
matriz_de_correlaciones = cor(datos)
round(matriz_de_correlaciones,2)

              edad_1c n_mat_1c n_fallas_1c prom_1c n_mat_1y2c n_fallas_1y2c prom_1y2c
edad_1c          1.00    -0.20       -0.13    0.03      -0.25         -0.12      0.03
n_mat_1c        -0.20     1.00        0.08    0.08       0.76          0.02      0.12
n_fallas_1c     -0.13     0.08        1.00   -0.06       0.13          0.79     -0.06
prom_1c          0.03     0.08       -0.06    1.00       0.05         -0.06      0.80
n_mat_1y2c      -0.25     0.76        0.13    0.05       1.00          0.07      0.11
n_fallas_1y2c   -0.12     0.02        0.79   -0.06       0.07          1.00     -0.08
prom_1y2c        0.03     0.12       -0.06    0.80       0.11         -0.08      1.00

desc_mat_cor = eigen(matriz_de_correlaciones)
autovalores_cor = desc_mat_cor$values
round(autovalores_cor,2)

[1] 2.05 1.94 1.50 0.87 0.24 0.21 0.19

#cuanta variabilidad concentra cada autovalor?
#es lo mismo calcular esa variabilidad con los autovalores de una matriz u otra?
variabilidad_cor = autovalores_cor/sum(autovalores_cor)
round(variabilidad_cor,2)

[1] 0.29 0.28 0.21 0.12 0.03 0.03 0.03

• Componentes principales

#comando que ejecuta el metodo de componentes principales
datos.pc = prcomp(datos,scale = TRUE)
#datos.pc$sdev #raiz cuadrada de los autovalores
round(datos.pc$sdev^2,2)

[1] 2.05 1.94 1.50 0.87 0.24 0.21 0.19

round(datos.pc$rotation,2) #autovectores (en columna)

                PC1   PC2   PC3   PC4   PC5   PC6   PC7
edad_1c        0.31 -0.07  0.16 -0.93  0.04 -0.01  0.01
n_mat_1c      -0.52 -0.17 -0.38 -0.26 -0.68  0.07 -0.15
n_fallas_1c   -0.41  0.37  0.43 -0.10  0.05  0.69  0.16
prom_1c       -0.12 -0.57  0.41  0.07 -0.18 -0.15  0.66
n_mat_1y2c    -0.54 -0.14 -0.37 -0.20  0.67 -0.14  0.19
n_fallas_1y2c -0.37  0.39  0.46 -0.07 -0.08 -0.68 -0.18
prom_1y2c     -0.15 -0.58  0.37  0.05  0.20  0.15 -0.67

round(datos.pc$center,2) #vector de medias (de casualidad coinciden)

      edad_1c      n_mat_1c   n_fallas_1c       prom_1c    n_mat_1y2c n_fallas_1y2c     prom_1y2c 
        20.09          1.80          0.44          6.63          3.45          0.94          6.66 

round(datos.pc$scale,2) #vector de desvios

      edad_1c      n_mat_1c   n_fallas_1c       prom_1c    n_mat_1y2c n_fallas_1y2c     prom_1y2c 
         2.45          0.99          0.71          1.35          2.02          1.18          1.11 

#loadings
carga1 = data.frame(cbind(X=1:length(datos),
                          primeracarga=data.frame(datos.pc$rotation)[,1]))
carga2 = data.frame(cbind(X=1:length(datos),
                          segundacarga=data.frame(datos.pc$rotation)[,2]))
round(cbind(carga1,carga2),2)

ggplot(carga1, aes(X,primeracarga) ,
       fill=tramo ) + geom_bar ( stat="identity" ,
       position="dodge" ,
       fill ="royalblue" ,
       width =0.5 ) + xlab( 'Tramo' ) + ylab('Primeracarga ' )

ggplot( carga2 , aes ( X , segundacarga ) ,
        fill =X ) + geom_bar ( stat="identity" , position="dodge" ,
           fill ="royalblue" ,
           width =0.5 ) +
xlab('Tramo') + ylab('Segundacarga')

ggbiplot(datos.pc, obs.scale=1 ,var.scale=1,alpha=0.0) #cambiando el alfa?

#calculo de scores de cada individuo
CP1 = as.matrix(datos_estandarizados)%*%as.matrix(carga1[2])
CP2 = as.matrix(datos_estandarizados)%*%as.matrix(carga2[2])
datos_estandarizados$CP1 = CP1
datos_estandarizados$CP2 = CP2

ggbiplot(datos.pc, obs.scale=0.5 ,var.scale=1,
         alpha=0.5,groups=factor(datos_completos$sit_ac)) +
  scale_color_manual(name="situacion academica", values=c("black", "red","green"),labels=c("cursantes", "desertores","egresados")) +  
theme(legend.direction ="horizontal", legend.position = "top")

library(MASS)
datos_train = datos_completos[datos_completos$sit_ac  %in% c(2,3),]    #se realiza un analisis discriminante con las clases 2 y 3
modelo_lda <- lda(formula = as.factor(sit_ac) ~ 
                  edad_1c    + n_mat_1c + n_fallas_1c   + prom_1c 
                + n_mat_1y2c + n_fallas_1y2c + prom_1y2c , data = datos_train)
modelo_lda

Call:
lda(as.factor(sit_ac) ~ edad_1c + n_mat_1c + n_fallas_1c + prom_1c + 
    n_mat_1y2c + n_fallas_1y2c + prom_1y2c, data = datos_train)

Prior probabilities of groups:
        2         3 
0.7139551 0.2860449 

Group means:
   edad_1c n_mat_1c n_fallas_1c  prom_1c n_mat_1y2c n_fallas_1y2c prom_1y2c
2 20.55818 1.464506   0.2722853 6.553522   2.615559      0.660940  6.573485
3 19.36003 2.360032   0.1270227 6.484709   5.292071      0.263754  6.667348

Coefficients of linear discriminants:
                      LD1
edad_1c       -0.05140860
n_mat_1c      -0.18549629
n_fallas_1c   -0.26415157
prom_1c       -0.11196861
n_mat_1y2c     0.63107883
n_fallas_1y2c -0.25074001
prom_1y2c      0.03414477

#predecimos la clase en base al modelo entrenado (o sea, como se comportan los mismos datos usados para el modelo)
pred_lda_train <- predict(modelo_lda,datos_train)
table(datos_train$sit_ac, pred_lda_train$class, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    2    3
         2 2729  356
         3  430  806

pred_lda_train_2 = ifelse(pred_lda_train$posterior[,1]>0.1,2,3)    #que significa esta linea? que es el 0.5??
table(datos_train$sit_ac, pred_lda_train_2, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    2    3
         2 3030   55
         3 1008  228

#posible intento de entender como se comporta la clase frente a las variables originales.
#Se puede hacer algo mejor?
library(klaR) 

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

partimat(as.factor(sit_ac) ~ edad_1c + n_mat_1c + n_fallas_1c + prom_1c, 
         data = datos_train, method = "lda", prec = 200,
         image.colors = c("darkgoldenrod1", "snow2"),col.mean = "firebrick")

pred_lda_completos <- predict(modelo_lda,datos_completos)

ggbiplot(datos.pc, obs.scale=0.1 ,var.scale=1,
         alpha=1,groups=factor(pred_lda_completos$class)) +
  scale_color_manual(name="prediccion lda", values=c("red","green"),labels=c("desertores","egresados")) +  
theme(legend.direction ="horizontal", legend.position = "top")

modelo_qda <- qda(as.factor(sit_ac) ~ 
                  edad_1c    + n_mat_1c + n_fallas_1c   + prom_1c 
                + n_mat_1y2c + n_fallas_1y2c + prom_1y2c , data = datos_train)
modelo_qda

Call:
qda(as.factor(sit_ac) ~ edad_1c + n_mat_1c + n_fallas_1c + prom_1c + 
    n_mat_1y2c + n_fallas_1y2c + prom_1y2c, data = datos_train)

Prior probabilities of groups:
        2         3 
0.7139551 0.2860449 

Group means:
   edad_1c n_mat_1c n_fallas_1c  prom_1c n_mat_1y2c n_fallas_1y2c prom_1y2c
2 20.55818 1.464506   0.2722853 6.553522   2.615559      0.660940  6.573485
3 19.36003 2.360032   0.1270227 6.484709   5.292071      0.263754  6.667348

pred_qda_train <- predict(modelo_qda,datos_train)
table(datos_train$sit_ac, pred_qda_train$class, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    2    3
         2 2450  635
         3  211 1025

# Y cuanto daba con el lda??
table(datos_train$sit_ac, pred_lda_train$class, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    2    3
         2 2729  356
         3  430  806

library(caret)

Loading required package: lattice
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

confusion_lda = confusionMatrix(factor(datos_train$sit_ac), pred_lda_train$class)
confusion_lda

Confusion Matrix and Statistics

          Reference
Prediction    2    3
         2 2729  356
         3  430  806
                                          
               Accuracy : 0.8181          
                 95% CI : (0.8063, 0.8295)
    No Information Rate : 0.7311          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.5465          
                                          
 Mcnemar's Test P-Value : 0.009219        
                                          
            Sensitivity : 0.8639          
            Specificity : 0.6936          
         Pos Pred Value : 0.8846          
         Neg Pred Value : 0.6521          
             Prevalence : 0.7311          
         Detection Rate : 0.6316          
   Detection Prevalence : 0.7140          
      Balanced Accuracy : 0.7788          
                                          
       'Positive' Class : 2               
                                          

confusion_qda = confusionMatrix(factor(datos_train$sit_ac), pred_qda_train$class)
confusion_qda

Confusion Matrix and Statistics

          Reference
Prediction    2    3
         2 2450  635
         3  211 1025
                                         
               Accuracy : 0.8042         
                 95% CI : (0.7921, 0.816)
    No Information Rate : 0.6158         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.5653         
                                         
 Mcnemar's Test P-Value : < 2.2e-16      
                                         
            Sensitivity : 0.9207         
            Specificity : 0.6175         
         Pos Pred Value : 0.7942         
         Neg Pred Value : 0.8293         
             Prevalence : 0.6158         
         Detection Rate : 0.5670         
   Detection Prevalence : 0.7140         
      Balanced Accuracy : 0.7691         
                                         
       'Positive' Class : 2              
                                         

pred_qda_completos <- predict(modelo_qda,datos_completos)

ggbiplot(datos.pc, obs.scale=0.5 ,var.scale=1,
         alpha=1,groups=factor(pred_qda_completos$class)) +
  scale_color_manual(name="prediccion qda", values=c("red","green"),labels=c("desertores","egresados")) +  
theme(legend.direction ="horizontal", legend.position = "top")

#regresión logistica
modelo_lg <- glm(as.factor(sit_ac) ~ 
                  edad_1c    + n_mat_1c + n_fallas_1c   + prom_1c 
                + n_mat_1y2c + n_fallas_1y2c + prom_1y2c , data = datos_train,family=binomial)
modelo_lg

Call:  glm(formula = as.factor(sit_ac) ~ edad_1c + n_mat_1c + n_fallas_1c + 
    prom_1c + n_mat_1y2c + n_fallas_1y2c + prom_1y2c, family = binomial, 
    data = datos_train)

Coefficients:
  (Intercept)        edad_1c       n_mat_1c    n_fallas_1c        prom_1c     n_mat_1y2c  n_fallas_1y2c      prom_1y2c  
      2.89307       -0.24816       -0.30848       -0.46255       -0.16191        0.88365       -0.52519       -0.06329  

Degrees of Freedom: 4320 Total (i.e. Null);  4313 Residual
Null Deviance:      5173 
Residual Deviance: 3287     AIC: 3303

pred_lg_train  <- predict(modelo_lg,type = "response")
clase_lg_train  = ifelse(pred_lg_train>0.5,1,0)  #ojo que el modelo genera la clase con 0 y 1 (no con 2 y 3)  
table(datos_train$sit_ac, clase_lg_train, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    0    1
         2 2764  321
         3  464  772

#regresión logistica 2. Modelo mas simple sacando las que parecían tener correlación.
#no es desable tener multicolinealidad de variables regresoras.
modelo2_lg <- glm(as.factor(sit_ac) ~ 
                  edad_1c    + n_mat_1c + n_fallas_1c   + prom_1c  , data = datos_train,family=binomial)

glm.fit: fitted probabilities numerically 0 or 1 occurred

modelo2_lg

Call:  glm(formula = as.factor(sit_ac) ~ edad_1c + n_mat_1c + n_fallas_1c + 
    prom_1c, family = binomial, data = datos_train)

Coefficients:
(Intercept)      edad_1c     n_mat_1c  n_fallas_1c      prom_1c  
    10.7459      -0.6190       0.8745      -0.9513      -0.1391  

Degrees of Freedom: 4320 Total (i.e. Null);  4316 Residual
Null Deviance:      5173 
Residual Deviance: 4050     AIC: 4060

pred2_lg_train  <- predict(modelo2_lg,type = "response")
clase2_lg_train  = ifelse(pred2_lg_train>0.5,1,0)    #que significa esta linea? que es el 0.5??

table(datos_train$sit_ac, clase2_lg_train, dnn = c("Clase real","Clase predicha"))

          Clase predicha
Clase real    0    1
         2 2840  245
         3  686  550

#Modelo support vector machine svm
library(e1071)
modelo_svm=svm(as.factor(sit_ac)~edad_1c    + n_mat_1c + n_fallas_1c   + prom_1c 
                + n_mat_1y2c+n_fallas_1y2c+prom_1y2c,
               data=datos_train,method="C-classification",kernel="radial",cost=10,gamma=.1)
pred_svm=predict(modelo_svm, datos_train)
table(datos_train$sit_ac, pred_svm, dnn = c("Clase real", "Clase predicha"))

          Clase predicha
Clase real    2    3
         2 2713  372
         3  340  896

error_svm<- mean(datos_train$sit_ac!= pred_svm) * 100
error_svm

[1] 16.47767

#comparacion de predicciones de diferentes modelos:
predicciones_varias = cbind(pred_lda_train$posterior[,2],
                            pred_qda_train$posterior[,2],
                            pred_lg_train,
                            pred2_lg_train)
cor(predicciones_varias)

                                   pred_lg_train pred2_lg_train
               1.0000000 0.9059432     0.9905854      0.7607068
               0.9059432 1.0000000     0.9295371      0.7554916
pred_lg_train  0.9905854 0.9295371     1.0000000      0.7833219
pred2_lg_train 0.7607068 0.7554916     0.7833219      1.0000000

#analisis de cluster
#este primer bloque es solo para definir funciones.
library(cluster)
library(pracma)

Attaching package: ‘pracma’

The following object is masked from ‘package:e1071’:

    sigmoid

# se define funcion de escalamiento disferente de la tipica normal.
esc01 <- function(x) { (x - min(x)) / (max(x) - min(x))} 
# se define una funcion para calcular metricas que orientan sobre el numero de clusters a elegir para el problema.
metrica = function(datA_esc,kmax,f) {
  
  sil = array()
  #sil_2 = array()
  sse = array()
  
  datA_dist= dist(datA_esc,method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
  for ( i in  2:kmax) {
    if (strcmp(f,"kmeans")==TRUE) {   #centroide: tipico kmeans
      CL  = kmeans(datA_esc,centers=i,nstart=50,iter.max = kmax)
      sse[i]  = CL$tot.withinss 
      CL_sil = silhouette(CL$cluster, datA_dist)
      sil[i]  = summary(CL_sil)$avg.width
        }
    if (strcmp(f,"pam")==TRUE){       #medoide: ojo porque este metodo tarda muchisimo 
      CL = pam(x=datA_esc, k=i, diss = F, metric = "euclidean")
      sse[i]  = CL$objective[1] 
      sil[i]  = CL$silinfo$avg.width
      }
  }
  sse
  sil
  return(data.frame(sse,sil))
}

#en este bloque se estudia cuantos clusters convendría generar respecto a indicadores tipicos de Clustering --> por ejemplo el "Silhouette"
kmax = 7
#2 opciones de escalamiento
#m1   = metrica(apply(datos,2,esc01),kmax,"kmeans")      #definida en la funcion esc01
m1   = metrica(scale(datos),kmax,"kmeans")               #tipica de la normal

did not converge in 7 iterationsdid not converge in 7 iterations

#graficos de los indicadores de clustering
par(mfrow=c(2,1))
plot(2:kmax, m1$sil[2:kmax],col=1,type="b", pch = 19, frame = FALSE, 
     xlab="Number of clusters K",
     ylab="sil") 

plot(2:kmax, m1$sse[2:kmax],type="b", pch = 19, frame = FALSE, 
     xlab="Number of clusters K",
     ylab="sse") 

par(mfrow=c(1,1))

#elegimos realizar 3 grupos
CL  = kmeans(scale(datos),5,nstart=50,iter.max = 10)
#CL  = kmeans(apply(datos,2,esc01),3,nstart=50,iter.max = 10)
datos$kmeans = CL$cluster

#en cuales 2 variables me conviene visualizar el cluster?
plot(datos$edad_1c,datos$prom_1c,col=datos$kmeans)+
grid()

integer(0)

#conviene en un biplot ya que tengo las flechas de las variables originales
ggbiplot(datos.pc, obs.scale=1 ,var.scale=1, alpha=0.5,groups = as.factor(datos$kmeans) )+
theme(legend.direction ="horizontal", legend.position = "top")

#lo hacemos finalmente para 3 grupos y lo visualizamos en un biplot
CL  = kmeans(scale(datos),3,nstart=50,iter.max = 10)
datos$kmeans = CL$cluster

ggbiplot(datos.pc, obs.scale=1 ,var.scale=1, alpha=0.5,groups = as.factor(datos$kmeans) )+
  scale_color_manual(name="Cluster kmeans", values=c("orange", "cyan","grey"),labels=c("grupo 1", "grupo 2","grupo 3")) +  
theme(legend.direction ="horizontal", legend.position = "top")

#cluster jerárquico
datos2=datos[,-8]#quito columna "kmeans"
datos2=scale(datos2)

# Matriz de distancias euclídeas 
mat_dist <- dist(x = datos2, method = "euclidean") 

# Dendrogramas (según el tipo de segmentación jerárquica aplicada)  
hc_complete <- hclust(d = mat_dist, method = "complete") 
hc_average  <- hclust(d = mat_dist, method = "average")
hc_single   <- hclust(d = mat_dist, method = "single")
hc_ward     <- hclust(d = mat_dist, method = "ward.D2")
#calculo del coeficiente de correlacion cofenetico
cor(x = mat_dist, cophenetic(hc_complete))

[1] 0.5846528

cor(x = mat_dist, cophenetic(hc_average))

[1] 0.7259612

cor(x = mat_dist, cophenetic(hc_single))

[1] 0.552105

cor(x = mat_dist, cophenetic(hc_ward))

[1] 0.4236285

# construccion de un dendograma usando los resultados de la técnica de Ward
plot(hc_ward )#no se ve bien si hay muchos datos
rect.hclust(hc_ward, k=3, border="red")#con 3 grupos

grupos<-cutree(hc_ward,k=3)#con 3 grupos
#split(rownames(datos),grupos)#devuelve una lista con las observaciones separadas por grupo

#visualizamos con 3 grupos el cluster jerarquico de la misma forma que kmeans
ggbiplot(datos.pc, obs.scale=1 ,var.scale=1, alpha=0.5,groups = as.factor(grupos) )+
  scale_color_manual(name="Cluster jerárquico Ward", values=c("orange", "cyan","grey"),labels=c("grupo 1", "grupo 2","grupo 3")) +  
theme(legend.direction ="horizontal", legend.position = "top")