PCA

PCA, análisis de componentes principales en inglés, es una técnica del análisis multivariante para identificar patrones en los datos de estudio, generación de indices representativos y disminución de la dimensión del data set inicial.PCA forma parte del aprendizaje no supervisado y suele ser uno de los pasos previos a realizar antes de la creación de un modelo.

Vamos a elegir un set de datos y aplicar la técnica PCA para su estudio,una clasificación académica de universidades del mundo publicado por el Times Higher Education (THE).

Esta publicación afiliada al periódico The Times se publica semanalmente en forma de revista. Cuenta además, con una página web que permite consultar sus informes. (THE)

El data set obtenido en formato csv, cuenta con 2.063 observacionesy 14 variables.


Relizaremos una exploración incial de los datos con visualizaciones que nos ayuden a entender la naturaleza de los datos.

Los paquetes que utilizaremos para ello son:
  • dplyr: Procesamiento de data frames.
  • corrplot Visualización de la matriz de correlaciones.
  • highcharter Wrapper de la librería D3 de javascript para visualizaciones dinámicas.
  • DT Generación de tablas dinámicas.


    library(dplyr)
    library(corrplot)
    library(highcharter)
    library(DT)
    library(webshot)
    library(magrittr)


    Existen al menos 50 observaciones con variables no informadas.Descartamos asignar un valor a los mismos que requeriría un modelo de imputación, así que la eliminamos de nuestro data set.


    Por medio del paquete DT creamos una tabla que permita ordenar y visualizar las características del set de datos de forma cómoda.

    options(DT.options = list(pageLength = 5),autoWidth = TRUE)
    datatable(fileUniversity,rownames=FALSE,filter='top',
              class = 'cell-border stripe',colnames=c('Nombre','Pais','Teaching Rating','Inter Outlook Rating','Research Rating','Citations Rating','Industry Income Rating','Numero Estudiantes','Student/Staff','% Estudiantes Internacionales','% Mujeres','Region'))


    Una primera valoración a la vista de los valores recogidos en la tabla es que hay un grupo de universidades excepcionales o élite, fundamentalmente americanas y británicas en ámbitos relacionados con la investigación y calidad en la enseñanza. Las variables ratio de enseñanza, de investigación,de citas en papers de investigación, el retorno económico generado deben guardar correlación. Utilizamos el paquete higcharter para hacer una visualizción multidimensional de estas variables.


    library(magrittr)
    library(highcharter)
    highchart() %>% 
      hc_title(text = "Teaching Rating vs Research Rating",align = "left",
               style = list(color = "#173B0B", fontWeight = "bold")) %>% 
      hc_subtitle(text = "Tamaño del punto por retorno industrial de las investigaciones",
                  align = "left",
                  style = list(color = "#2b908f", fontWeight = "bold")) %>%
      hc_add_series_scatter(fileUniversity$Teaching_Rating,
                            fileUniversity$Research_Rating,
                            fileUniversity$Industry_Income_Rating, 
                            fileUniversity$Citations_Rating,
                            dataLabels = list(
                            enabled = FALSE,
                            format = "{point.label}")) %>%
      hc_yAxis(title = list(text = "Rating en Investigación"),
               labels = list(format = "{value}%"), max = 110) %>%
      hc_xAxis(title = list(text = "Rating en Eseñanza"),
               labels = list(format = "{value}%"), max = 100) %>%
      hc_chart(zoomType = "xy") %>% 
      hc_tooltip(useHTML = TRUE,
                 headerFormat = "<table>",
                 pointFormat = paste("<tr><th colspan=\"1\"><b>{point.label}</b></th></tr>",
                                     "<tr><th>Teaching Rating:</th><td>{point.x} </td></tr>",
                                     "<tr><th>Research Rating:</th><td>{point.y} </td></tr>",
                                     "<tr><th>Nº Estudiantes:</th><td>{point.z} </td></tr>",
                                     "<tr><th>Nº Citas:</th><td>{point.valuecolor}</td></tr>"),
                 footerFormat = "</table>") %>%
      hc_credits(enabled = TRUE, text = "Source: Times Higher Education (THE)",
                 href = "https://www.timeshighereducation.com",
                 style = list(fontSize = "12px"))


    La grafica evidencia una correlacion entre las variables Research y Teaching Rating, el color del punto evidencia la correlación igualmente con la variable citations rating pero no tan destacadamente con el tamaño del punto, el retorno industrial. Agrupamos por país la variable research rating para cuantificar la primera impresión de dominio de las universidades americanas.Representamos la distribución de las 100 mejores universidades de esa variable por país. EEUU sobresale de forma destacada.España no coloca ninguna de sus universidades entre las 100 primeras.

    fileUniversity<-fileUniversity %>% 
                    arrange(desc(fileUniversity$Research_Rating))
    ranking<-fileUniversity[c(1:100),] %>%
             group_by(Country) %>%
             summarise(N.Universidades=n()) %>%
             arrange(N.Universidades)
    
    highchart() %>% 
      hc_title(text = "Número de Universidades por país entre las 100 primeras") %>%
      hc_subtitle(text = "Ranking por research") %>% 
      hc_add_series_labels_values(ranking$Country, ranking$N.Universidades, name = "Pie",
                                  colorByPoint = TRUE, type = "column") %>%
      hc_yAxis(title = list(text = "Número de Universidades"),
               labels = list(format = "{value}"), max = 50) %>% 
      hc_xAxis(categories = ranking$Country) %>% 
      hc_legend(enabled = FALSE) %>% 
      hc_tooltip(pointFormat = "{point.y}") 


    Ampliamos la representación anterior a las 250 primeras universidades en el rating de investigación. España coloca a dos universidades. La Univeridad Autónoma y la Pompeu Fabra de Barcelona.

    data(worldgeojson)
    fileRanking<-fileUniversity %>%
             arrange(desc(Research_Rating))
    ranking<-fileRanking[c(1:250),] %>%
             group_by(Country) %>%
             summarise(N.Universidades=n()) %>%
             arrange(desc(N.Universidades))
    names(ranking)<-c('name','N.Universidades')
    highchart() %>%  
      hc_add_series_map(worldgeojson, ranking,name='Nº Universidades',
                        value = "N.Universidades", joinBy = "name") %>%
      hc_title(text = "Número de Universidades por país entre las 250 primeras",
               align = "right",
               style = list(color = "#173B0B", fontWeight = "bold")) %>%
        hc_subtitle(text = "ranking por ratio de investigacion",
                  align = "right",
                  style = list(color = "#2b908f", fontWeight = "bold")) %>%
       hc_colorAxis(minColor = "#81F79F", maxColor = "#0B610B") %>%
      #hc_colorAxis(stops = color_stops()) %>% 
      hc_legend(valueDecimals = 0, valueSuffix = "%") %>%
      hc_mapNavigation(enabled = TRUE) 


    Analizamos otras dimensiones del set de datos, su tamaño.Número de estudiantes y la proporcion de estos con el staff de la universidad.

    hchart(fileUniversity$Num_Students)
    hchart(fileUniversity$Student.Staff_Ratio)


    Matriz de correlaciones

    fileUniversity$Num_Students <-log(fileUniversity$Num_Students)
    fileUniversity$Student.Staff_Ratio <-log(fileUniversity$Student.Staff_Ratio)
    fileUniversity<-fileUniversity[!is.na(fileUniversity$X._Inter_Students),]
    
    fileUniversitySc <- scale(fileUniversity[,3:11],scale = T) 
    correl <- cor(fileUniversitySc)
    corrplot(correl,order="hclust",tl.col='black',tl.cex=.80)

    PCA

    ######################### PCA
    ###########################
    
    #modo manual
    eS <- eigen(correl)
    
    eigen.val <- eS$values # Eigenvalues
    eigen.val
    ## [1] 3.19093527 1.70168588 1.36231365 1.04090837 0.61125878 0.51388539
    ## [7] 0.36290253 0.15203356 0.06407657
    eigen.vec <- eS$vectors # Eigenvectors
    eigen.vec
    ##              [,1]        [,2]        [,3]       [,4]        [,5]
    ##  [1,] -0.46020720 -0.20876175 -0.23434401  0.2809881  0.18632478
    ##  [2,] -0.40203181  0.33309181  0.29693016 -0.2838026  0.07075998
    ##  [3,] -0.48642468 -0.11282607 -0.29916373  0.1147555  0.13159888
    ##  [4,] -0.43500092  0.06142952 -0.01001537  0.2227012 -0.31203782
    ##  [5,] -0.19556492 -0.32894275 -0.33667624 -0.4971119 -0.57374776
    ##  [6,]  0.04703465  0.38721967 -0.60767575  0.1431341  0.33318631
    ##  [7,]  0.09232108  0.44140665 -0.39566726 -0.4937097 -0.01806369
    ##  [8,] -0.38794056  0.24115486  0.35312765 -0.3195377  0.26876295
    ##  [9,] -0.03558141  0.56376162  0.05868297  0.4053096 -0.57691490
    ##              [,6]        [,7]         [,8]         [,9]
    ##  [1,]  0.05471861 -0.33506994  0.003329906  0.678794826
    ##  [2,]  0.07938086  0.11035881 -0.726968756  0.082084580
    ##  [3,] -0.07123452 -0.33686967 -0.076137476 -0.711563651
    ##  [4,] -0.55233838  0.56742124  0.177395056  0.037726870
    ##  [5,]  0.38295893  0.12200672  0.023280002  0.042524100
    ##  [6,]  0.36067085  0.46561208 -0.011159059 -0.008704135
    ##  [7,] -0.50597963 -0.33120354  0.091738894  0.147924448
    ##  [8,]  0.26631743  0.01417201  0.648461988 -0.026090227
    ##  [9,]  0.27725647 -0.31388556  0.068219937 -0.017525881
    prop.var <- eigen.val / sum(eigen.val) # Proportion of variability
    prop.var
    ## [1] 0.354548363 0.189076209 0.151368183 0.115656486 0.067917643 0.057098376
    ## [7] 0.040322503 0.016892618 0.007119619
    prop.var.accum <- cumsum(eigen.val) / sum(eigen.val) # Proportion of accumulated variability
    prop.var.accum
    ## [1] 0.3545484 0.5436246 0.6949928 0.8106492 0.8785669 0.9356653 0.9759878
    ## [8] 0.9928804 1.0000000
    # obtenemos los scores
    Z <- fileUniversitySc %*% eigen.vec[, 1:4]
    
    # representacion de la matriz de correlaciones
    Z.2 <- fileUniversitySc %*% eigen.vec[, 1:9]
    corrplot(cor(Z.2),order="hclust",tl.col='black',tl.cex=.80)

    par(mfrow=c(1,2))
    barplot(prop.var,col='blue',las=2,main='% de var. por comp.ppal')
    plot(prop.var,type='l',col='blue',main='screeplot')

    dev.off()
    ## null device 
    ##           1
    h<-Z[,1:2]
    h<-data.frame(h)
    h$name<-fileUniversity[,1]
    h$region<-fileUniversity[,12]
    
    #scatterplot <- ggplot(h, 
            #              aes(x = X1, 
             #                 y = X2,
              #                color=region,
               #               alpha=0.5)) 
    #scatterplot + geom_point()