Machine learning

Referencia principal

An Introduction to Machine Learning with R, Laurent Gatto (2020)

El aprendizaje automático tiene cuatro fases:

Recogida de datos
Preprocesamiento de los datos (limpieza, factorización)
Generación de un modelo
Evaluación del modelo

Es un proceso iterativo, en el cuál se repiten los pasos previos a medida que se va teniendo más información sobre los datos originales y el rendimiento del modelo.

En el machine learning supervisado (SML), los inputs están categorizados e indican los resultados que se quieren obtener.

El machine learning no supervisado (UML) tiene un input numérico; sirve para explorar la existencia de patrones y estructuras dentro de los datos. No hay objetivos fijados previos al análisis.

Es una técnica fundalmente descriptiva que se usa en data mining [1].

flowchart LR
  I{input} --- |categórico| SML("<b>SML</b><br>Supervisado")
  I --- |numérico| UML("<b>UML</b><br>No supervisado")
  UML --- |categórico| KM("<em>Clustering</em>") --- k-clusters
  KM --- jerárquico
  UML --- |numérico| RD("Reducción<br>dimensional") --- PCA
  RD --- t-SNE  
  SML --- |categórico| CL("Clasificación") --- KNN
  CL --- SVM
  SML --- |numérico| RE("Regresión") --- lm[Regresión lineal]
  RE --- rl[Regresión logística]

Tareas básicas de machine learning

Al abordar un proyecto hay que pensar en el tipo de machine learning que va a ser necesaria:

Clasificación
Predicción numérica
Detección de patrones
Clustering

Supervised Machine Learning (SML)

En SML, el input es categórico. En estos modelos hay objetivos en la predicción que se modelan a partir de categorías previas o muestras de entrenamiento.

flowchart LR
  in((<em>in</em>)) --- |categórico| out(("<em>out</em>")) --- |numérico | cat(Regresión)
  out --- |categórico| clas(Clasificación)
  clas --- bin[binaria]
  clas --- mul[multi-variable]

Modelos de regresión

regresión lineal
regresión logística

Modelos de clasificación

máquinas de vectores de soporte (support vector machines, SVN)

SML - Regresión

Regresión lineal simple

Cuando buscamos correlaciones entre variables es útil examinar todas las variables de form simultánea. La función pairs() sirve para visualizar variables de par en par en busca de relaciones entre ellas:

bw <- MASS::birthwt
pairs(~ bw$age + bw$bwt + bw$race)

El coeficiente de correlación de Pearson se puede obtener con cor():

# Una correlación simple:
#           cor(bw$age, bw$bwt)
#
# Una matriz de correlacion:
cor(bw[, c("age", "bwt", "race")])

             age         bwt       race
age   1.00000000  0.09031781 -0.1728180
bwt   0.09031781  1.00000000 -0.1947135
race -0.17281795 -0.19471349  1.0000000

El modelo lineal de regresión hace una estimación por el método de estimación por mínimos cuadrados:

modelo <- lm(age ~ bwt, data = bw)
summary(modelo)


Call:
lm(formula = age ~ bwt, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.8920 -3.9614 -0.7409  3.1685 20.4196 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.131e+01  1.605e+00   13.27   <2e-16 ***
bwt         6.563e-04  5.292e-04    1.24    0.216    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.291 on 187 degrees of freedom
Multiple R-squared:  0.008157,  Adjusted R-squared:  0.002853 
F-statistic: 1.538 on 1 and 187 DF,  p-value: 0.2165

lm() devuelve un objeto devuelve toda la información sobre el análisis, con un coeficiente de correlación que habla sobre la bondad del ajuste ($R^{2}$) - RMSE (root mean squared error). Es una medida de cuanto en qué medida impacta la variable independiente sobre la variabilidad de los datos observada.

Cómo reportar R²

Cuando $R^{2}$ se obtiene entre dos variables, $\sqrt{R^{2}} = r$ ($r$ siendo el coeficiente de correlación de Pearson).

Esto devuelve una información para una recta, donde $y$ es la variable dependiente, $x$ es la variable independiente.

\[ y = \alpha x + \beta \]

$\alpha$ es la pendiente y $\beta$ es el punto en el que la línea se cruza con el eje Y (intercepción).

plot(bw$bwt, bw$age, xlab = "bwt", ylab = "age")
abline(modelo)

Regresión lineal múltiple

Se pueden explorar todas las variables de un dataframe al mismo tiempo

cor(x = bw, method = "pearson")

              low         age         lwt         race       smoke          ptl
low    1.00000000 -0.11893928 -0.16962694  0.137792751  0.16140431  0.196087267
age   -0.11893928  1.00000000  0.18007315 -0.172817953 -0.04434618  0.071606386
lwt   -0.16962694  0.18007315  1.00000000 -0.165048544 -0.04417908 -0.140029003
race   0.13779275 -0.17281795 -0.16504854  1.000000000 -0.33903074  0.007951293
smoke  0.16140431 -0.04434618 -0.04417908 -0.339030745  1.00000000  0.187557063
ptl    0.19608727  0.07160639 -0.14002900  0.007951293  0.18755706  1.000000000
ht     0.15237025 -0.01583700  0.23636040  0.019929917  0.01340704 -0.015399579
ui     0.16904283 -0.07515558 -0.15276317  0.053602088  0.06215900  0.227585340
ftv   -0.06296026  0.21539394  0.14052746 -0.098336254 -0.02801314 -0.044429660
bwt   -0.78480516  0.09031781  0.18573328 -0.194713487 -0.19044806 -0.154653390
               ht          ui         ftv         bwt
low    0.15237025  0.16904283 -0.06296026 -0.78480516
age   -0.01583700 -0.07515558  0.21539394  0.09031781
lwt    0.23636040 -0.15276317  0.14052746  0.18573328
race   0.01992992  0.05360209 -0.09833625 -0.19471349
smoke  0.01340704  0.06215900 -0.02801314 -0.19044806
ptl   -0.01539958  0.22758534 -0.04442966 -0.15465339
ht     1.00000000 -0.10858506 -0.07237255 -0.14598189
ui    -0.10858506  1.00000000 -0.05952341 -0.28392741
ftv   -0.07237255 -0.05952341  1.00000000  0.05831777
bwt   -0.14598189 -0.28392741  0.05831777  1.00000000

multimodel <- lm(bwt ~ low + 
                 age + 
                 lwt + 
                 race + 
                 smoke + 
                 ptl + 
                 ht + 
                 ui + 
                 ftv, data = bw)
summary(multimodel)


Call:
lm(formula = bwt ~ low + age + lwt + race + smoke + ptl + ht + 
    ui + ftv, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-991.22 -300.96   -5.39  277.74 1637.80 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3612.508    229.457  15.744  < 2e-16 ***
low         -1131.217     73.957 -15.296  < 2e-16 ***
age            -6.245      6.347  -0.984 0.326416    
lwt             1.051      1.133   0.927 0.355085    
race         -100.905     38.544  -2.618 0.009605 ** 
smoke        -174.116     72.000  -2.418 0.016597 *  
ptl            81.340     68.552   1.187 0.236980    
ht           -181.955    137.661  -1.322 0.187934    
ui           -336.776     93.314  -3.609 0.000399 ***
ftv            -7.578     30.992  -0.245 0.807118    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 433.7 on 179 degrees of freedom
Multiple R-squared:  0.6632,    Adjusted R-squared:  0.6462 
F-statistic: 39.16 on 9 and 179 DF,  p-value: < 2.2e-16

Hecho esto, se pueden seleccionar las mejores correlaciones utilizando el criterio de información de Aikaike o AIC utilizando step

step(object = multimodel, 
     directio = "both",
     trace = 1)

Start:  AIC=2305.09
bwt ~ low + age + lwt + race + smoke + ptl + ht + ui + ftv

        Df Sum of Sq      RSS    AIC
- ftv    1     11246 33682899 2303.2
- lwt    1    161709 33833361 2304.0
- age    1    182159 33853812 2304.1
- ptl    1    264837 33936490 2304.6
- ht     1    328638 34000291 2304.9
<none>               33671653 2305.1
- smoke  1   1100074 34771727 2309.2
- race   1   1289179 34960832 2310.2
- ui     1   2450185 36121838 2316.4
- low    1  44009294 77680946 2461.1

Step:  AIC=2303.15
bwt ~ low + age + lwt + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- lwt    1    154400 33837299 2302.0
- age    1    205543 33888442 2302.3
- ptl    1    268799 33951698 2302.7
- ht     1    319845 34002744 2302.9
<none>               33682899 2303.2
+ ftv    1     11246 33671653 2305.1
- smoke  1   1095243 34778142 2307.2
- race   1   1280252 34963151 2308.2
- ui     1   2442561 36125460 2314.4
- low    1  44045331 77728230 2459.2

Step:  AIC=2302.02
bwt ~ low + age + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- age    1    158960 33996259 2300.9
- ht     1    229750 34067049 2301.3
- ptl    1    232303 34069602 2301.3
<none>               33837299 2302.0
+ lwt    1    154400 33682899 2303.2
+ ftv    1      3937 33833361 2304.0
- smoke  1   1132887 34970186 2306.2
- race   1   1427486 35264785 2307.8
- ui     1   2527757 36365056 2313.6
- low    1  45611646 79448944 2461.3

Step:  AIC=2300.9
bwt ~ low + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- ptl    1    189792 34186051 2299.9
- ht     1    228505 34224763 2300.2
<none>               33996259 2300.9
+ age    1    158960 33837299 2302.0
+ lwt    1    107817 33888442 2302.3
+ ftv    1     19622 33976637 2302.8
- smoke  1   1054854 35051113 2304.7
- race   1   1301334 35297593 2306.0
- ui     1   2450537 36446796 2312.1
- low    1  45476560 79472818 2459.4

Step:  AIC=2299.95
bwt ~ low + race + smoke + ht + ui

        Df Sum of Sq      RSS    AIC
- ht     1    236483 34422534 2299.3
<none>               34186051 2299.9
+ ptl    1    189792 33996259 2300.9
+ age    1    116450 34069602 2301.3
+ lwt    1     85589 34100462 2301.5
+ ftv    1     21990 34164061 2301.8
- smoke  1    941319 35127371 2303.1
- race   1   1268492 35454543 2304.8
- ui     1   2279585 36465636 2310.2
- low    1  45501243 79687294 2457.9

Step:  AIC=2299.26
bwt ~ low + race + smoke + ui

        Df Sum of Sq      RSS    AIC
<none>               34422534 2299.3
+ ht     1    236483 34186051 2299.9
+ ptl    1    197771 34224763 2300.2
+ age    1    114599 34307935 2300.6
+ lwt    1     24448 34398086 2301.1
+ ftv    1     12780 34409754 2301.2
- smoke  1    935878 35358412 2302.3
- race   1   1269413 35691947 2304.1
- ui     1   2122280 36544815 2308.6
- low    1  48004723 82427257 2462.3


Call:
lm(formula = bwt ~ low + race + smoke + ui, data = bw)

Coefficients:
(Intercept)          low         race        smoke           ui  
    3586.50     -1139.20       -97.34      -157.42      -303.19

El modelo óptimo es el que tiene el AIC más bajo.

Akaike information criterion

Este algoritmo intenta minimizar la pérdida de información de un modelo (bondad del ajuste o goodness of fit), intentando al mismo tiempo que sea lo más simple posible.

Utiliza la función de verosimilitud para la bondad de ajuste, penalizando el número de variables.

Wikipedia

Podemos tomar la información del análisis previo y generar un nuevo modelo

multimodel2 <- lm(bwt ~ low +
                  race + 
                  smoke + 
                  ui, data = bw)
summary(multimodel2)


Call:
lm(formula = bwt ~ low + race + smoke + ui, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-1025.8  -351.0    30.8   285.8  1500.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3586.50      86.68  41.379  < 2e-16 ***
low         -1139.20      71.12 -16.019  < 2e-16 ***
race          -97.34      37.37  -2.605 0.009942 ** 
smoke        -157.42      70.38  -2.237 0.026510 *  
ui           -303.19      90.02  -3.368 0.000922 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 432.5 on 184 degrees of freedom
Multiple R-squared:  0.6557,    Adjusted R-squared:  0.6482 
F-statistic: 87.59 on 4 and 184 DF,  p-value: < 2.2e-16

Evaluación del modelo

La forma de evaluar la bondad del modelo generado es utilizando la raíz del error cuadrático medio (RMSE; root mean squared edrror) en regresión; y precisión predictiva del modelo para la clasificación.

Para calcular estas medidas no se deben utilizar muestras que contengan elementos incluidos en la muestra de entrenamiento porque sobreestiman la bondad del modelo.

Rendimiento en regresión

El paquete caret (Classification And REgression Training) contiene muchas herramientas útiles en SML, como caret::predict.

in-sample y out-of-sample error

Cuando el modelo se valida contra sí mismo decimos que la raíz del error cuadrático medio es in-sample error - pertenece a la misma muestra con la que se ha entrando el modelo y es mucho menos precisa:

library(caret)

# Generamos un modelo con una regresión lineal
diam <- lm(price ~ ., diamonds)

# Generamos la predicción
diam_pred <- predict(diam, diamonds)

# Obtenemos RMSE 
error <- diam_pred - diamonds$price
rmse_in <- sqrt(mean(error^2))
rmse_in

[1] 1129.843

Para poner a prueba el modelo con una muestra ajena a la generación del mismo(out-of-sample error), entrenamos al modelo con una fracción de la muestra y la validamos con la fracción restante:

# 80% de la muestra para entrenar y 20% para validar:
# Crea números secuenciales del primer al último elemento
all <- 1:length(diamonds$price)

# Selecciona el 80% para entrenar el modelo
trn <- sample(length(diamonds$price), 
              length(diamonds$price)*0.8)

# Entrena el modelo con la muestra de entrenamiento...
diam2 <- lm(price ~ ., diamonds[trn,])

# ... lo valida con la reserva
diam2_pred <- predict(diam2, diamonds[-trn,])

# ... y obtiene el RMSE
error2 <- diam2_pred - diamonds$price[-trn]
rmse_out <- sqrt(mean(error2^2))
rmse_out

[1] 1134.417

Muestras pequeñas

Estos métodos no son fiables en muestras pequeñas. La presencia de un solo outlier en estas circunstancias puede influir de forma desmesurada en el modelo.

La validación cruzada es un método más robusto para demostrar la bondad de un modelo.

Validación cruzada

La validación cruzada realiza el mismo proceso anterior un número $n$ de veces; una vez se demuestra la bondad de los $n$ modelos parciales, se genera un modelo definitivo con toda la muestra.

Para este proceso se utiliza caret::train.

Note

caret también se encarga de que las proporciones de las clases se mantengan constantes en las muestras que se utilizan para generar los sub-modelos (folds).

`caret::train`

model <- train(price ~ ., diamonds,
  # algoritmo del modelo
  method = "lm",       
       # 'cv' - cross-validation
  trControl = trainControl(
       method = "cv",
       # número de modelos parciales que se crean
       number = 10,
       # esta opción puede estar en TRUE en modo interactivo
       verboseIter = FALSE))
model

Linear Regression 

53940 samples
    9 predictor

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 48545, 48547, 48547, 48546, 48545, 48547, ... 
Resampling results:

  RMSE      Rsquared   MAE     
  1131.048  0.9196972  740.5421

Tuning parameter 'intercept' was held constant at a value of TRUE

diam_cv <- predict(model, diamonds)
error <- diam_cv - diamonds$price
rmse_cv <- sqrt(mean(error^2))
rmse_cv

[1] 1129.843

SML - Clasificación

Diferencias entre SVN y KNN

Investigar un poco… hay artículos en internet que lo explican

El modelo entrenado con una muestra de entrenamiento se pone a prueba con una muestra distinta de la misma población (muestra de prueba).

En ocasiones se segmenta la muestra original (una parte se utiliza para entrenar el modelo y la otra para comprobar su bondad).

flowchart LR
     original["Muestra o<br>población"] ==> train
     original --> test
     train["Datos de<br>entrenamiento"] ==> m(("Modelo"))
     test("Datos de prueba") -.-> m -.-> testc("Datos de prueba<br>categorizados")
     test ---> eval(("Evaluación"))
     testc -.-> eval

Máquinas de vectores de soporte (SVM)

Respecto a KNN

útil en espacios multi-dimensionales
útil cuando hay más dimensiones que casos
evitar overfitting
más sensible al ruido
computacionalmente complejo

SVM es útil en espacios multi-dimensionales o donde el número de dimensiones sea mayor que el de casos.

Generar training y testing data

caTools::sample.split() recibe un vector factorizado y un porcentaje. Devuelve dos grupos donde las proporciones del factor están preservadas (sirve para crear grupos de entrenamiento y comprobación).

split <- caTools::sample.split(iris$Species, SplitRatio = 0.7)
data_train <- iris[split, ]
data_test <- iris[!split, ]

Modelado

El paquete e1071 contiene una función para modelar SVM

svm_iris <- e1071::svm(Species ~ ., data = data_train)

Se realiza las predicciones

pred_iris <- predict(svm_iris, newdata = data_test)
pred_iris

         2          4         19         26         28         29         30 
    setosa     setosa     setosa     setosa     setosa     setosa     setosa 
        32         34         36         39         40         45         47 
    setosa     setosa     setosa     setosa     setosa     setosa     setosa 
        48         52         55         56         57         63         64 
    setosa versicolor versicolor versicolor versicolor versicolor versicolor 
        65         69         70         72         85         86         92 
versicolor versicolor versicolor versicolor versicolor versicolor versicolor 
        93         94        102        104        106        108        116 
versicolor versicolor  virginica  virginica  virginica  virginica  virginica 
       121        122        126        131        133        136        137 
 virginica  virginica  virginica  virginica  virginica  virginica  virginica 
       142        145        146 
 virginica  virginica  virginica 
Levels: setosa versicolor virginica

Evaluación del modelo

caret

El paquete caret (Classification And REgression Training) contiene muchas herramientas útiles en SML

Después comprobamos las predicciones con las etiquetas contenidas en los datos de prueba:

conf_matr_iris <- caret::confusionMatrix(pred_iris, data_test$Species)
conf_matr_iris

Confusion Matrix and Statistics

            Reference
Prediction   setosa versicolor virginica
  setosa         15          0         0
  versicolor      0         15         0
  virginica       0          0        15

Overall Statistics
                                     
               Accuracy : 1          
                 95% CI : (0.9213, 1)
    No Information Rate : 0.3333     
    P-Value [Acc > NIR] : < 2.2e-16  
                                     
                  Kappa : 1          
                                     
 Mcnemar's Test P-Value : NA         

Statistics by Class:

                     Class: setosa Class: versicolor Class: virginica
Sensitivity                 1.0000            1.0000           1.0000
Specificity                 1.0000            1.0000           1.0000
Pos Pred Value              1.0000            1.0000           1.0000
Neg Pred Value              1.0000            1.0000           1.0000
Prevalence                  0.3333            0.3333           0.3333
Detection Rate              0.3333            0.3333           0.3333
Detection Prevalence        0.3333            0.3333           0.3333
Balanced Accuracy           1.0000            1.0000           1.0000

k-nearest neighbours (KNN)

Respecto a SVM:

el modelo se construye durante la fase de predicción
menos sensible al ruido
sensible a outliers
peor rendimiento en espacios multidimensionales

La clasificación se realiza midiendo la distancia entre observaciones no categorizadas e infieriendo sus clases por sus vecinos más cercanos

La función knn entrena el modelo y hace las predicciones en un mismo tiempo

knn(train,    # el modelo de entrenamiento
    test,     # los datos de prueba
    cl,       # las etiquetas de 'train'
    k = 1)    # el número de 'vecinos' que se consideran

# Selecciona 50 filas aleatorias de 'iris' para los datos de entrenamiento...
trn <- sample(150,50)

# ... y para poner a prueba el modelo
new <- sample(150,50)

# Entrena el modelo
iris_knn <- class::knn(iris[trn, -5], 
                       iris[new,-5],
                       iris$Species[trn],
                       k = 20,
                       prob = TRUE)

# Comprueba:
table(iris_knn,             # - las conclusiones del modelo
      iris$Species[new])    # - con la realidad

            
iris_knn     setosa versicolor virginica
  setosa         17          0         0
  versicolor      0         15         4
  virginica       0          1        13

# Esto devuelve la tasa de aciertos: 
mean(iris_knn == iris$Species[new])

[1] 0.9

Random forest

Los modelos random forest son modelos no-lineales robustos que previene el over-fitting, pero requiere un ajuste manual de los parámetros.

Están basados en árboles de decisiones, que segmentan los datos y los clasifican en cada rama. Un número excesivo de ramas lleva a over-fitting (considerar que el ruido de la muestra es relevante).

Para evitar esto se puede “podar” el árbol limitando el número de iteraciones posibles, o limitando el número de ramas que se consideran una vez generado el modelo.

Unsupervised Machine Learning (UML)

El objetivo del aprendizaje no supervisado es exploratorio. Encuentra estructuras y patrones ocultos en los datos.

El clustering agrupa los datos en base a similitudes entre ellos. La reducción de dimensionalidad trata de simplificar los datos para generar hipótesis o visualizarlos de forma más sencilla.

Clustering no jerárquico

Agrupa datos de la muestra en clusters en base a similitudes entre sus elementos individuales. El modelo intenta generar la menor suma de cuadrados de cada observación a su centro (total within cluster sum of squares).

k-means clustering

stats::kmeans(data,
              centers = 3,
              iter.max = 10,
              nstart = 10)

center = n genera n centros, tres valores que actúan de media para cada grupo. Cada elemento es asignado a uno de los grupos de forma aleatoria
El modelo asigna cada elemento al grupo que tenga su media más cerca; genera una nueva media y repite el proceso hasta que los centros dejan de moverse (o se alcanza un número máximo de iteraciones, iter.max = n)
nstart = n modifica el elemento aleatorio de kmeans, que marca la posición inicial de cada centro; realiza n posiciones iniciales para mejorar el modelo

Generar el modelo

x <- subset(iris, select = c(Petal.Length, Sepal.Length))
cl <- stats::kmeans(x, centers = 3, nstart = 100)
cl

K-means clustering with 3 clusters of sizes 51, 58, 41

Cluster means:
  Petal.Length Sepal.Length
1     1.492157     5.007843
2     4.393103     5.874138
3     5.678049     6.839024

Clustering vector:
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
 21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
 41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
  1   1   1   1   1   1   1   1   1   1   3   2   3   2   2   2   2   2   2   2 
 61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   3   3   2   2 
 81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2 
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
  3   2   3   3   3   3   2   3   3   3   3   3   3   2   2   3   3   3   3   2 
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
  3   2   3   2   3   3   2   2   3   3   3   3   3   3   3   3   3   3   2   3 
141 142 143 144 145 146 147 148 149 150 
  3   3   2   3   3   3   2   3   3   2 

Within cluster sum of squares by cluster:
[1]  9.893725 23.508448 20.407805
 (between_SS / total_SS =  90.5 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"

Método del codo

Para determinar el punto que mejor define las agrupaciones con este método se puede dibujar un gráfico que represente el número de grupos en el eje x y suma total de los cuadrados en el eje y. El punto en el que el gráfico se acoda indica el punto de mayor discriminación.

elbow <- c()
for (c in 1:15) {
    elbow[c] <- stats::kmeans(x,
                              centers = c)$tot.withinss
}

plot(elbow, type="b")

Resultados

ggplot(data = x, aes(x = Sepal.Length, y = Petal.Length,
       colour = factor(cl$cluster))) +
geom_point() + 
geom_point(data = cl$centers,
           aes(x = Sepal.Length, y = Petal.Length),
               color= "black", size = 7, shape = 4) +
scale_colour_discrete(name = "Species") + 
theme_classic()
ggplot(data = iris, aes(x = Sepal.Length, y = Petal.Length,
       colour = Species)) + 
geom_point() +
theme_classic()

Clustering jerárquico

En el clustering jerárquico tiene dos aproximaciones; la aglomerativa y la divisva.

La agrupación aglomerativa es útil cuando no se conoce con certeza el la estructura de la muestra y se quieren determinar cuantos grupos existen; considera cada elemento como un grupo, al cuál va añadiendo elementos similares de forma progresiva (hclust, agnes).

La agrupación divisiva es útil cuando se quiere controlar el número de grupos que se van a generar - segmenta la muestra original de forma progresiva hasta obtenerlos (diana).

Agrupamiento jerárquico aglomerativo

Se define un cluster con cada uno de los $n$ puntos de la muestra. En cada iteración se fusionan los dos cluster más próximos entre sí hasta que se consume la muestra. Las distancias utilizan el vector multidimensional de cada elemento.

Este proceso genera un dendrograma, que relaciona los grupos y las distancias entre ellos. El último definie la distancia que se considera relevante para diferenciar los grupos generados.

dist() calcula las distancias entre elementos (la distancia Euclidiana; devuelve una matriz de distancias). hclust() aplica el clustering jerárquico. Se puede representar después con plot():

iris_dist <- dist(iris[,1:4])
iris_clst <- hclust(iris_dist)
plot(iris_clst, main = "Agrupamiento jerárquico aglomerativo")

Definir los clusters

cutree corta el árbol definiendo h (la altura) o k (el número de clusters deseado) - devuelve un vector con un factor.

plot(iris_clst, main="", sub="", xlab="", ylab ="")
abline(h = 3.9, col="red")

identical(cutree(iris_clst, h = 3.9),
          cutree(iris_clst, k = 3))

[1] TRUE

Agrupamiento jerárquico divisivo

hclust también permite hacer clustering divisivo:

iris_div <- hclust(iris_dist, 
       method = "ward.D2")
plot(iris_div)

cutree(iris_div, k = 4)

  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 2 3 2 3 2 3 3 2 3 2 3 2 3 3 2 3 2 2 2 2
 [75] 2 2 2 4 2 3 3 3 3 2 3 2 2 2 3 3 3 2 3 3 3 3 3 2 3 3 4 2 4 4 4 4 3 4 4 4 4
[112] 4 4 2 2 4 4 4 4 2 4 2 4 2 4 4 2 2 4 4 4 4 4 2 2 4 4 4 2 4 4 4 2 4 4 4 2 4
[149] 4 2

iris_div_clust <- data.frame(iris[1:4], 
                             Clusters = cutree(iris_div, k = 4))
iris_div_clust[sample(1:150, 5), ]

    Sepal.Length Sepal.Width Petal.Length Petal.Width Clusters
102          5.8         2.7          5.1         1.9        2
16           5.7         4.4          1.5         0.4        1
63           6.0         2.2          4.0         1.0        3
9            4.4         2.9          1.4         0.2        1
8            5.0         3.4          1.5         0.2        1

Reducción dimensional

La reducción dimensional transforma los datos y los transfiere a un nuevo espacio dimensional más reducido.

Las nuevas dimensiones resumen las características de los datos originales. Esta simplificación permite

visualizar los datos de forma sencilla
encontrar patrones y estructuras en los datos
realizar pre-procesamiento de los datos y las simplifica. genera componentes principales (PC) independientes entre sí y los numera(PC1, PC2, PC3…)

Principal Component Analysis (PCA)

PCA es una técnica de reducción dimensional lineal que genera nuevas dimensiones de las dimensiones originales. Las nuevas dimensiones

son combinaciones lineales de las dimensiones originales
la primera dimensión (PC1) recoge la mayor variabilidad de los datos
- después PC2, PC3, PC4…
los componentes son independientes entre sí

`prcomp` (principal componente analysis)

prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE,
       tol = NULL, rank. = NULL, ...)

scale.= aparece como FALSE por defecto; pero en general tiene sentido escalar las varianzas antes (y help(prcomp) lo recomienda).

Preparar los datos

Muchos de los modelos requieren preparar los datos previamente

estandarizando todas las variables
eliminando las que no sean numéricas

Para la estandarización, escalamos las variables para obtener una distribución normal $X \sim N(0,1)$:

\[ z = \dfrac{x - \mu}{\sigma} \]

Aunque algunas funciones (como prcomp) permiten escalar los datos antes del análisis pasando la opción scale. = TRUE (o algo parecido), es buena idea revisarlos y prepararlos de antemano

dd <- survival::diabetic
head(dd,3)

  id laser age   eye trt risk  time status
1  5 argon  28  left   0    9 46.23      0
2  5 argon  28 right   1    9 46.23      0
3 14 xenon  12  left   1    8 42.50      0

dd_std <- scale(dd[, sapply(dd,is.numeric)])
head(dd_std,3)

            id        age        trt       risk      time     status
[1,] -1.752093  0.4873237 -0.9987302 -0.4731892 0.4987246 -0.8042944
[2,] -1.752093  0.4873237  0.9987302 -0.4731892 0.4987246 -0.8042944
[3,] -1.733930 -0.5928762  0.9987302 -1.1511404 0.3240656 -0.8042944

La presencia de un factor como variable del conjunto de datos puede alterar el análisis (se interpreta como numérico aunque no lo sea). Hay soluciones a este problema.

Los factores ordinales pueden ser aceptables.

En iris, difícil visualizar las relaciones entre grupos utilizando las cuatro dimensiones numéricas; utilizando PCA reducimos las dimensiones y facilitamos la interpretación:

pairs(iris[, -5], col = iris[,5], pch=19)

Las cuatro dimensiones numéricas de `iris`

Primero estandarizamos la muestra

head(iris,3)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa

iris_std <- scale(iris[, sapply(iris, is.numeric)])
iris_std <- as.data.frame(iris_std)
iris_std$Species <- iris$Species
head(iris_std,3)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1   -0.8976739   1.0156020    -1.335752   -1.311052  setosa
2   -1.1392005  -0.1315388    -1.335752   -1.311052  setosa
3   -1.3807271   0.3273175    -1.392399   -1.311052  setosa

Aplicamos PCA

irispca <- prcomp(iris_std[,-5], scale.=TRUE)
summary(irispca)

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

Proportion of Variance devuelve la cantidad de la varianza explicada por cada PC. En iris, PC1 y PC2 explican el 96% de la varianza (73% y 23% respectivamente).

plot(irispca)

pca$x contiene los vectores de cada PC para poder examinarlos posteriormente:

pcx <- irispca$x
plot(pcx[,1], pcx[,2], col = iris[,5], pch =19,
     xlab = "PC1", ylab = "PC2")
pcx <- irispca$x
plot(pcx[,3], pcx[,4], col = iris[,5], pch =19,
     xlab = "PC3", ylab = "PC4")

PC1 y PC2 explican el 96% de la varianza

biplot() devuelve los puntos orientados en PC1 y PC2, con unos vectores superpuestos que representan los datos originales:

biplot(irispca)

pca$rotation contiene información sobre cómo contribuyen las distintas variables a cada PC (en signo y cuantía), los cargos

irispca$rotation

                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

En este ejemplo, Sepal.Width contribuye de forma muy importante en PC2, de forma negativa.

t-SNE (t-distributed Stochastic Neighbour Embedding)

t-SNE es una técnica de reducción dimensional no lineal que intenta condensar una variable multidimensional a dos o tres dimensiones. Agrupa elementos próximos e ignora elementos distantes¹. Es útil para variables con un gran número de dimensiones.

`Rtsne::Rtsne`

Primero hay que eliminar duplicados de los datos con unique().

Los resultados de Rtsne son estocásticos - cada vez que se ejecuta devuelve resultados distintos.

uiris <- unique(iris[, 1:5])
iristsne <- Rtsne::Rtsne(uiris[,1:4])
plot(iristsne$Y, col = uiris$Species)

Escalar los datos

Al igual que en PCA, es posible escalar las dimensiones antes del análisis con pca_scale = TRUE. Los datos están centrados en 0 porque los datos están normalizados por defecto (normalize = TRUE). Ver escalar los datos.

Al igual que otros algoritmos (principalmente algoritmos de clasificación), t-SNE tiene dos parámetros que modifican sustancialmente los resultados:

Perplejidad: ??? Iteraciones; el número de iteraciones hasta que el proceso de clustering finaliza

Incompleto

An Introduction to Machine Learning with R, Laurent Gatto (2020)

Deep learning

El aprendizaje profundo utiliza aprendizaje automático basado en redes neuronales artificiales, inspiradas en el cerebro humano [song_inferring_2024?]. El término profundo indica la presencia de múltiples capas, al contrario que las redes neuronales clásicas.

Las redes neuronales convolucionales (CNN) están diseñadas para analizar datos estructurados como una cuadrícula (imágenes, secuencias temporales). Son muy utilziadas en computer vision y manipulación de imágenes.

Paquetes para CNN

keras
tensorflow

Las redes neuronales recurrentes (RNN) están optimizadas para el procesamiento de datos secuenciales y temporales - procesamiento del lenguaje natural, texto, audio…

Lantz B. Machine learning with R: Learn techniques for building and improving machine learning models, from data preparation to model tuning, evaluation, and working with big data, fourth edition. 4th ed. Place of publication not identified: Packt Publishing; 2023.

Chollet F. Deep learning with R. Second edition. Shelter Island: Manning; 2022.

Gatto L. An Introduction to Machine Learning with R. 2020.

Song Y, Millidge B, Salvatori T, Lukasiewicz T, Xu Z, Bogacz R. Inferring neural activity before plasticity as a foundation for learning beyond backpropagation. Nat Neurosci. 2024 Feb;27(2):348–58.

Jones E, Harden S, Crawley MJ. The R book. Third edition. Hoboken, NJ: Wiley; 2022.

Field A, Miles J, Field Z. Discovering statistics using R. Repr. Los Angeles, CA, USA: Sage; 2014.

Yap BW, Sim CH. Comparisons of various types of normality tests. Journal of Statistical Computation and Simulation. 2011 Dec;81(12):2141–55.

Perezgonzalez JD. Fisher, Neyman-Pearson or NHST? A tutorial for teaching data testing. Front Psychol. 2015 Mar;6:223.

Carmona Pontaque F. Álgebra Matricial en Estadística. Análisis Multivariante. Fundació Universitat Oberta de Catalunya; 2024.

10.

Carmona Pontaque F. Modelos lineales. Ediciones de la Universidad de Barcelona; 2004.

11.

Faraway JJ. Linear models with R. Second edition. Boca Raton London: CRC Press, Taylor & Francis; 2014. (Texts in statistical science).

https://en.wikipedia.org/wiki/T-distributed_stochastic_neighbor_embedding↩︎

Tareas básicas de machine learning

Supervised Machine Learning (SML)

SML - Regresión

Regresión lineal simple

Regresión lineal múltiple

Evaluación del modelo

Rendimiento en regresión

in-sample y out-of-sample error

Validación cruzada

caret::train

SML - Clasificación

Máquinas de vectores de soporte (SVM)

Modelado

Evaluación del modelo

k-nearest neighbours (KNN)

Random forest

Unsupervised Machine Learning (UML)

Clustering no jerárquico

k-means clustering

Generar el modelo

Método del codo

Resultados

Clustering jerárquico

Agrupamiento jerárquico aglomerativo

Definir los clusters

Agrupamiento jerárquico divisivo

Reducción dimensional

Principal Component Analysis (PCA)

prcomp (principal componente analysis)

t-SNE (t-distributed Stochastic Neighbour Embedding)

Rtsne::Rtsne

Deep learning

`caret::train`

`prcomp` (principal componente analysis)

`Rtsne::Rtsne`