MI ESTUDIO SOBRE R

Valores NA (missing values) y cómo rellenarlos

Partimos de uno de los datasets de pruebas que hay en RStudio.

#library(help="datasets") #lista completa (cogeremos el dataset llamado attenu, por ejemplo)
data(attenu)
attach(attenu)

The following objects are masked from attenu (pos = 8):

    accel, dist, event, mag, station

The following objects are masked from attenu (pos = 13):

    accel, dist, event, mag, station

The following objects are masked from attenu (pos = 27):

    accel, dist, event, mag, station

help(attenu)
boxplot(attenu, cex.axis=0.5) 

Vemos que la variable dist es asimetrica y tiene outliers. Vamos que las variables mag y accel son leptocúrticas (apretadas en torno a la media)

summary(attenu)

     event            mag           station         dist            accel        
 Min.   : 1.00   Min.   :5.000   117    :  5   Min.   :  0.50   Min.   :0.00300  
 1st Qu.: 9.00   1st Qu.:5.300   1028   :  4   1st Qu.: 11.32   1st Qu.:0.04425  
 Median :18.00   Median :6.100   113    :  4   Median : 23.40   Median :0.11300  
 Mean   :14.74   Mean   :6.084   112    :  3   Mean   : 45.60   Mean   :0.15422  
 3rd Qu.:20.00   3rd Qu.:6.600   135    :  3   3rd Qu.: 47.55   3rd Qu.:0.21925  
 Max.   :23.00   Max.   :7.700   (Other):147   Max.   :370.00   Max.   :0.81000  
                                 NA's   : 16                                     

Matriz de parejas de features (CORRELACIONES)

Hacemos una análisis exploratorio de las correlaciones. Matriz de correlaciones (entre variables) y funcion de densidad de probabilidad (univariable):

library(car)
scatterplotMatrix(x=attenu)

O también construirla usando una función personalizada por nosotros:

pairs(attenu, gap=0, lower.panel = panel.smooth, upper.panel = function(x,y){
  panel.smooth(x,y)
  par(usr=c(0,1,0,1))
  correlacion <- cor(x,y,use="complete.obs")
  text(0.6,0.7,col="blue", cex=1.2, round(correlacion, digits=2))
}  )

Es importante ver si hay alguna correlación fácilmente visible (lineal…).

IDs de los elementos que tienen valores NA en cualquiera de sus features (columnas):

which(is.na(attenu$event))

integer(0)

which(is.na(attenu$mag))

integer(0)

which(is.na(attenu$station))

 [1]  79  81  94  96  99 107 108 114 116 118 123 126 128 155 156 160

which(is.na(attenu$dist))

integer(0)

which(is.na(attenu$accel))

integer(0)

Vemos que sólo la variable station tiene valores NA. También se pueden coger sólo las filas en las que todas las features están rellenas (sin ningún valor NA):

attenu[complete.cases(attenu),]

Si queremos ELIMINAR aquellas filas (margin=1 significa ROWS) que tengan más de un valor NA: 1. Aplicar la función SUM sobre filas, contando el número de elementos NA. 2. Si hay uno o más, saco sus índices (which) 3. Creo un vector con esos índices y se lo resto al dataset attenu.

filas_con_na <- apply(attenu, MARGIN = 1, FUN = function(x){ sum(is.na(x)) >=1})
indices_de_filas_con_na <- which(filas_con_na)
datos_sinhuecos=attenu[ -c( indices_de_filas_con_na ) , ]

Comprobamos que ya no hay HUECOS (con algún valor NA):

attenu[!complete.cases(datos_sinhuecos),]

ANALISIS DE COMPONENTES PRINCIPALES (PCA)

Si tenemos muchas features, conviene reducirlas para que el análisis de correlaciones sea más sencillo. Es MUY importante indicar que trabajamos con la matriz de correlaciones (en vez de la de covarianzas): cor=TRUE Además, PCA sólo trabaja con features que sean NUMERICAS, así que debemos comprobar que todas las columnas (features) son numericas

drops <- c("station")
datos_para_pca <- attenu[ , !(names(attenu) %in% drops)]
datos_sin_na <- na.omit(datos_para_pca) # limpiar los NA
sapply(datos_sin_na, class) #clase/modo de cada feature

    event       mag      dist     accel 
"numeric" "numeric" "numeric" "numeric" 

sapply(datos_sin_na, typeof) #tipo de cada feature

   event      mag     dist    accel 
"double" "double" "double" "double" 

pca_attenu <- princomp(datos_sin_na, cor=TRUE)  
pca_attenu

Call:
princomp(x = datos_sin_na, cor = TRUE)

Standard deviations:
   Comp.1    Comp.2    Comp.3    Comp.4 
1.4752709 1.0373500 0.6701831 0.5462007 

 4  variables and  182 observations.

plot(pca_attenu)

summary(pca_attenu)

Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4
Standard deviation     1.4752709 1.0373500 0.6701831 0.54620074
Proportion of Variance 0.5441061 0.2690237 0.1122864 0.07458381
Cumulative Proportion  0.5441061 0.8131298 0.9254162 1.00000000

biplot(pca_attenu) #PC2 vs PC1

#Los pesos de cada individuo (fila) proyectado en las componentes (PC1, PC2...)
pca_pesos <- pca_attenu$scores
pca_pesos_pc1pc2 <- pca_pesos[, 1:2] #Pesos sólo en las componentes que más influyen en la varianza total
plot(pca_pesos_pc1pc2[,1], pca_pesos_pc1pc2[,2]) #misma gráfica que veiamos en el biplot

Vemos que la componente 1 tiene un peso del 54.41% sobre las varianzas; la 2 tiene 26.9%; etc. Cogemos la PC1 y PC2, que suman un 81% de la influencia en la varianza total.

Vemos también los pesos que tienen las variables de entrada (features) dentro de las componentes compuestas (PC1, PC2…)., mirando el plot de flechas:

• PC1 (eje horizontal) está afectada por las variables: event, dist, mag y accel. (es decir, por todas las features de entrada, pero en proporción diferente).
• PC2 (eje vertical): está afectada por mag y accel, principalmente.

CLUSTERING (habiendo hecho PCA antes)

Vamos a mostrar un subtipo de “clustering jerárquico” llamado Agglomerative Nesting (AGNES).

Hay tres tipos de métodos de unión en clustering:

• Single linkage: según la MINIMA distancia entre dos elementos de clusters diferentes.
• Complete linkage: según la MAXIMA distancia entre dos elementos de clusters diferentes.
• Average linkage: según la distancia MEDIA entre dos elementos de clusters diferentes.

library(cluster)
agnes_single <- agnes(x = datos_sin_na, method = "single")
plot(agnes_single) # Cuando veamos espacios entre las muestras, es la separacion entre clusters. Es DIFICIL de ver.

Con Single, vemos que hay 3 clusters (si height=50), fijándose en los huecos que hay en el gráfico “banner”, que separan a los clusters.

library(cluster)
agnes_complete <- agnes(x = datos_sin_na, method = "complete")
plot(agnes_complete) # Cuando veamos espacios entre las muestras, es la separacion entre clusters. Es DIFICIL de ver.

Con Complete, se ven claramente 3 clusters (si height=200).

library(cluster)
agnes_average <- agnes(x = datos_sin_na, method = "average")
plot(agnes_average) # Cuando veamos espacios entre las muestras, es la separacion entre clusters. Es DIFICIL de ver.

Con Average, se ven 3 clusters (si height=100). Donde más claro lo vemos es con el método “complete”, así que creamos los clusters vistos con COLORES. Esos colores los pintamos sobre los pesos que ya habiamos calculado para PCA (de las componentes PC1 y PC2).

complete_3_clusters=cutree(agnes_complete,3)
plot(pca_pesos,col=complete_3_clusters)

REGRESIÓN LINEAL SIMPLE

Link: http://rpubs.com/joser/RegresionSimple Datos entrada:

grasas_tabla <- read.table('http://verso.mat.uam.es/~joser.berrendero/datos/EdadPesoGrasas.txt', header = TRUE)
names(grasas_tabla)

[1] "peso"   "edad"   "grasas"

pairs(grasas_tabla)

cor(grasas_tabla)

            peso      edad    grasas
peso   1.0000000 0.2400133 0.2652935
edad   0.2400133 1.0000000 0.8373534
grasas 0.2652935 0.8373534 1.0000000

Vemos correlación (0.83) entre las variables edad y grasas, como era de esperar a simple vista. El gráfico muestra también la relación lineal (a ojo).

Recta de mínimos cuadrados: Usamos un MODELO LINEAL (lm) para calcular una regresión lineal muy sencilla, de grasas respecto de la edad: \(y =102.575+5.321*x\)

O mejor aun: \(y = 102.575 (+/-29.6376) + 5.321(+/-0.7243) * x\)

regresion <- lm(grasas ~ edad, data = grasas_tabla) #grasas dependiente de la edad
summary(regresion)

Call:
lm(formula = grasas ~ edad, data = grasas_tabla)

Residuals:
    Min      1Q  Median      3Q     Max 
-63.478 -26.816  -3.854  28.315  90.881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 102.5751    29.6376   3.461  0.00212 ** 
edad          5.3207     0.7243   7.346 1.79e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 43.46 on 23 degrees of freedom
Multiple R-squared:  0.7012,    Adjusted R-squared:  0.6882 
F-statistic: 53.96 on 1 and 23 DF,  p-value: 1.794e-07

plot(grasas_tabla$edad, grasas_tabla$grasas, xlab='Edad', ylab='Grasas')
abline(regresion)

plot(regresion)

Vemos que el R² es 0.7012 –> Bondad de la recta de ajuste.

PREDICCION PARA NUEVAS ENTRADAS:

Usar la recta:

nuevas.edades <- data.frame(edad = seq(30, 50))
prediccion_nueva<-predict(regresion, nuevas.edades)
prediccion_nueva

       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15 
262.1954 267.5161 272.8368 278.1575 283.4781 288.7988 294.1195 299.4402 304.7608 310.0815 315.4022 320.7229 326.0435 331.3642 336.6849 
      16       17       18       19       20       21 
342.0056 347.3263 352.6469 357.9676 363.2883 368.6090 

Explicación del modelo de regresión simple:

• La columna t-value es el estadístico t –> Cociente entre cada estimador y su error típico. Permite realizar los contrastes de hipótesis nula (que sean 0) –> H₀:β₀=0 y H₀:β₁=0
• Los p-valores aparecen en la columna Pr(>|t|) –> Probabilidad de superar el umbral –> Como son muy pequeños (menores que 0.01 ó 0.05, utilizados normalmente) ==> Se rechazan ambas hipótesis nulas.
• El estimador de la desviación típica (σ) de los errores aparece como Residual standard error y su valor en el ejemplo es 43.5
• Los intervalos de confianza para los parámetros se obtienen con el comando confint. El parámetro level permite elegir el nivel de confianza (por defecto es 0.95):

confint(regresion)

                2.5 %     97.5 %
(Intercept) 41.265155 163.885130
edad         3.822367   6.818986

confint(regresion, level = 0.90)

                  5 %       95 %
(Intercept) 51.780153 153.370132
edad         4.079335   6.562018

Los “intervalos de confianza para la respuesta media” y los “intervalos de predicción para la respuesta” se pueden obtener usando el comando predict. Por ejemplo, el siguiente código calcula y representa los dos tipos de intervalos para el rango de edades que va de 20 a 60 años (los de predicción en rojo):

nuevas.edades <- data.frame(edad = seq(20, 60))
# Grafico de dispersion y recta
plot(grasas_tabla$edad, grasas_tabla$grasas, xlab='Edad', ylab='Grasas')
abline(regresion)

# Intervalos de confianza de la respuesta media:
# ic es una matriz con tres columnas: la primera es la prediccion, las otras dos son los extremos del intervalo
ic <- predict(regresion, nuevas.edades, interval = 'confidence')
lines(nuevas.edades$edad, ic[, 2], lty = 2, col = 'green')
lines(nuevas.edades$edad, ic[, 3], lty = 2, col = 'green')

# Intervalos de prediccion
ic <- predict(regresion, nuevas.edades, interval = 'prediction')
lines(nuevas.edades$edad, ic[, 2], lty = 2, col = 'red')
lines(nuevas.edades$edad, ic[, 3], lty = 2, col = 'red')

Análisis de la varianza:

anova(regresion)

Analysis of Variance Table

Response: grasas
          Df Sum Sq Mean Sq F value    Pr(>F)    
edad       1 101933  101933  53.964 1.794e-07 ***
Residuals 23  43444    1889                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Diagnóstico del modelo:

residuos <- rstandard(regresion)
valores.ajustados <- fitted(regresion)
plot(valores.ajustados, residuos)

#Prueba de normalidad (sencilla) para saber si la muestra tiene una distribucion NORMAL
qqnorm(residuos)
qqline(residuos)

Como tiene pocas muestras, podríamos haber pintado una t-Student encima, para poder comparar.

Otras pruebas de normalidad –> Kolmogorov-Smirnov y Shapiro-Wilk.

Ejemplo interesante:

x <- 1:20
w <- 1 + sqrt(x)/2 #desviaciones estandar (pesos)
y <- x + w*rnorm(x)
dummy <-- data.frame(x=x, y=y)
fm <- lm(y~x, data=dummy) #MODELO simple linear regression
summary(fm)

Call:
lm(formula = y ~ x, data = dummy)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4670 -1.4875 -0.3557  1.7932  7.1525 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.7369     1.4954   0.493    0.628    
x             1.0565     0.1248   8.464 1.09e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.219 on 18 degrees of freedom
Multiple R-squared:  0.7992,    Adjusted R-squared:  0.788 
F-statistic: 71.63 on 1 and 18 DF,  p-value: 1.089e-07

fm1 <- lm(y~x, data=dummy, weight=1/w^2) #MODELO con weighted regression
summary(fm1)

Call:
lm(formula = y ~ x, data = dummy, weights = 1/w^2)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.8149 -0.7074 -0.1920  0.7042  2.3323 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.0708     1.0537   1.016    0.323    
x             1.0880     0.1074  10.126 7.36e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.164 on 18 degrees of freedom
Multiple R-squared:  0.8507,    Adjusted R-squared:  0.8424 
F-statistic: 102.5 on 1 and 18 DF,  p-value: 7.359e-09

Resto del análisis:

attach(dummy) #hacemos visibles las columnas del dataframe, para trabajar directamente con sus nombres.

The following objects are masked _by_ .GlobalEnv:

    x, y

lrf <-lowess(x,y, f = 2/3, iter = 3, delta = 0.01 * diff(range(x))) #Funcion de regresion local no parametrica: LOWESS smoother (locally-weighted polynomial regression)
summary(lrf)

  Length Class  Mode   
x 20     -none- numeric
y 20     -none- numeric

plot(x,y)
lines(x, lrf$y, col="blue") #pintamos la linea de la REGRESION LINEAL

abline(0,1, lty=3, col="magenta") #linea de regresion verdadera (intercept 0, slope=1) (linea de puntos)
abline(coef(fm), col="green") #linea de regresion sin pesos

abline(coef(fm1), col="red") #linea de regresion con pesos

detach()
plot( fitted(fm), resid(fm), xlab="Fitted values", ylab="Residuals", main="Residuals vs Fitted" ) #Analisis de heterodicidad

qqnorm( resid(fm), main="Residuals Rankit Plot" ) #Analisis de skewness, kurtosis y outliers (no muy util en este caso)

#rm( fm, fm1, lrf, x, dummy) #limpieza de variables

Manual de CRAN-INTRO.ap.8.2

Lo básico:

attach(faithful)

The following objects are masked from faithful (pos = 8):

    eruptions, waiting

The following objects are masked from faithful (pos = 13):

    eruptions, waiting

The following objects are masked from faithful (pos = 24):

    eruptions, waiting

summary(eruptions) #Como he hecho attach, evito hacer esto: summary(faithful$eruptions)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.600   2.163   4.000   3.488   4.454   5.100 

fivenum(eruptions) #Pinta un "dibujo" hecho con numeros, que parece un histograma

[1] 1.6000 2.1585 4.0000 4.4585 5.1000

stem(eruptions)

  The decimal point is 1 digit(s) to the left of the |

  16 | 070355555588
  18 | 000022233333335577777777888822335777888
  20 | 00002223378800035778
  22 | 0002335578023578
  24 | 00228
  26 | 23
  28 | 080
  30 | 7
  32 | 2337
  34 | 250077
  36 | 0000823577
  38 | 2333335582225577
  40 | 0000003357788888002233555577778
  42 | 03335555778800233333555577778
  44 | 02222335557780000000023333357778888
  46 | 0000233357700000023578
  48 | 00000022335800333
  50 | 0370

Dibujos:

hist(eruptions,seq(from=1.6, to=5.2, by=0.2), probability = TRUE)
lines(density(eruptions, bw=0.1))

rug(eruptions)

La función de densidad (CDF) es:

plot(ecdf(eruptions), do.points=FALSE, verticals = TRUE)

Si cogemos sólo aquellos con X>3:

erupciones_largas<-eruptions[eruptions>3]
hist(erupciones_largas)

plot(ecdf(erupciones_largas), do.points=FALSE, verticals = TRUE)
x<-seq(3, 5.4, 0.01)
#Usando la distribucion NORMAL (norm)
lines(x, pnorm(x, mean = mean(erupciones_largas), sd=sqrt(var(erupciones_largas))), lty=3)  #linea de puntos(lty=3)

PRUEBAS DE NORMALIDAD:

1. Q-Q (Quantile-Quantile)

par(pty="s") #pinta la caja cuadrada (en vez de que sea rectangular por defecto)
qqnorm(erupciones_largas)
qqline(erupciones_largas)

#La COMPARAMOS con una t-Student de 5 grados de libertad:
puntos_t_student <- rt(250, df = 5)
qqnorm(puntos_t_student)
qqline(puntos_t_student)

#Pintamos el qqplot
qqplot(qt( ppoints(250), df=5 ), puntos_t_student, xlab = "QQplot para la densidad t-Student")
qqline(puntos_t_student)

2. Test Shapiro-Wilk:

shapiro.test(erupciones_largas)

    Shapiro-Wilk normality test

data:  erupciones_largas
W = 0.97934, p-value = 0.01052

3. Test Kolmogorov-Smirnov:

ks.test(erupciones_largas, "pnorm", mean=mean(erupciones_largas), sd=sqrt(var(erupciones_largas)))

ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  erupciones_largas
D = 0.066133, p-value = 0.4284
alternative hypothesis: two-sided

REGRESION LOGISTICA

challenger <- read.table('http://verso.mat.uam.es/~joser.berrendero/datos/challenger.txt', header = TRUE)
table(challenger$defecto)

 0  1 
16  7 

colores <- NULL
colores[challenger$defecto==0] <- 'green'
colores[challenger$defecto==1] <- 'red'
plot(challenger$temp, challenger$defecto, pch = 21, bg = colores, xlab = 'Temperatura', ylab = 'Probabilidad de defectos')
legend('bottomleft', c('No defecto', 'Si defecto'), pch = 21, col = c('green', 'red'))

DISTRIBUCIONES DE PROBABILIDAD

Hay muchas: binomial (binom), Poisson (pois), normal (norm), exponencial (exp), t-Student (t), chi-cuadrado (chisq), F (f)…

MAQUINAS DE VECTOR SOPORTE (SVM)

Busca el hiperplano que separe las observaciones según su clases (clasificación). Buena explicación: https://www.youtube.com/watch?v=ImhjKyc88x0

Link de mi ejemplo: http://rpubs.com/joser/svm

library(MASS)
library(e1071)
load(url('http://www.uam.es/joser.berrendero/datos/practica-svm-io.RData'))
# Mostramos 10 filas
head(breast.cancer2)

Gráficamente, pintamos la correlación entre las 2 features de entrada:

# Prepara los datos
x <- cbind(breast.cancer2$x.smoothness, breast.cancer2$x.concavepoints)
head(x)

        [,1]     [,2]
[1,] 0.09779 0.047810
[2,] 0.10750 0.031100
[3,] 0.10240 0.020760
[4,] 0.08983 0.029230
[5,] 0.08600 0.005917
[6,] 0.10310 0.027490

y <- breast.cancer2$y
head(y)

[1] 0 0 0 0 0 0
Levels: 0 1

n0 <- sum(y==0)
n1 <- sum(y==1)
# Para que los graficos queden mas bonitos (rojo = maligno, verde = benigno)
colores <- c(rep('green',n0),rep('red',n1))
pchn <- 21
# Diagrama de dispersion
plot(x, pch = pchn, bg = colores, xlab='smoothness', ylab='concavepoints')

SVM lineal

Siendo la variable y un factor, queremos predecir usando el resto de features.

Para calcular la regla de clasificación SVM lineal con C=10, se usa la función svm del paquete e1071. El primer argumento y~. indica que la variable y (que debe ser necesariamente un factor) se desea predecir en términos del resto de variables del fichero. La sintaxis es similar a la que utilizaríamos para ajustar un modelo lineal o un modelo de regresión logística. El segundo argumento indica el fichero en el que están las variables que vamos a usar. El argumento kernel corresponde al núcleo que representa el producto escalar que queremos utilizar. La opción linear corresponde a k(x,y)=x′y. El argumento cost determina la penalización que ponemos a los errores de clasificación. Con el fin de estimar la probabilidad de clasificar erróneamente una observación se puede utilizar validación cruzada, dividiendo la muestra en, por ejemplo, dos partes. Ello se consigue fijando cross=2. Finalmente, scale=FALSE se usa para usar los datos no estandarizados (por defecto, sí se estandarizan).

C <- 10
svm.lineal <- svm(y~., data=breast.cancer2, kernel='linear', cost=C, cross=2, scale=FALSE) #Su type es C-classification porque y es un factor
summary(svm.lineal)

Call:
svm(formula = y ~ ., data = breast.cancer2, kernel = "linear", cost = C, cross = 2, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 
      gamma:  0.5 

Number of Support Vectors:  258

 ( 129 129 )


Number of Classes:  2 

Levels: 
 0 1

2-fold cross-validation on training data:

Total Accuracy: 88.92794 
Single Accuracies:
 89.43662 88.42105 

x.svm <- x[svm.lineal$index,]
w <- crossprod(x.svm, svm.lineal$coefs)
w0 <- svm.lineal$rho
plot(x, pch = pchn, bg = colores, xlab='smoothness', ylab='concavepoints')
abline(w0/w[2], -w[1]/w[2], lwd=2, col='blue')

#Estadisticas: precisión y cobertura
table(predict(svm.lineal), breast.cancer2$y, dnn=c("Prediction", "Actual"))

          Actual
Prediction   0   1
         0 342  43
         1  15 169

library(caret)
confusionMatrix(breast.cancer2$y, predict(svm.lineal))

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 342  15
         1  43 169
                                          
               Accuracy : 0.8981          
                 95% CI : (0.8702, 0.9217)
    No Information Rate : 0.6766          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.776           
 Mcnemar's Test P-Value : 0.0003922       
                                          
            Sensitivity : 0.8883          
            Specificity : 0.9185          
         Pos Pred Value : 0.9580          
         Neg Pred Value : 0.7972          
             Prevalence : 0.6766          
         Detection Rate : 0.6011          
   Detection Prevalence : 0.6274          
      Balanced Accuracy : 0.9034          
                                          
       'Positive' Class : 0               
                                          

Indica que hay 258 vectores soporte. Usando validación cruzada con la muestra dividida en dos partes se estima una probabilidad de acierto en la clasificación de aproximadamente el 88%. Podemos cambiar el parámetro de penalización para ver si estos valores aumentan o disminuyen.

SVM cuadrático

svm.cuadratico <- svm(formula = y~., data = breast.cancer2, kernel='polynomial', degree=2, gamma=1, coef0=1, cost=C, cross=10, scale=FALSE)
summary(svm.cuadratico)

Call:
svm(formula = y ~ ., data = breast.cancer2, kernel = "polynomial", degree = 2, gamma = 1, coef0 = 1, cost = C, 
    cross = 10, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  polynomial 
       cost:  10 
     degree:  2 
      gamma:  1 
     coef.0:  1 

Number of Support Vectors:  211

 ( 106 105 )


Number of Classes:  2 

Levels: 
 0 1

10-fold cross-validation on training data:

Total Accuracy: 90.50967 
Single Accuracies:
 87.5 87.7193 91.22807 92.98246 92.98246 92.98246 91.22807 87.7193 87.7193 92.98246 

#Estadisticas: precisión y cobertura
table(predict(svm.cuadratico), breast.cancer2$y, dnn=c("Prediction", "Actual"))

          Actual
Prediction   0   1
         0 340  36
         1  17 176

library(caret)
confusionMatrix(breast.cancer2$y, predict(svm.cuadratico))

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 340  17
         1  36 176
                                          
               Accuracy : 0.9069          
                 95% CI : (0.8799, 0.9294)
    No Information Rate : 0.6608          
    P-Value [Acc > NIR] : < 2e-16         
                                          
                  Kappa : 0.7971          
 Mcnemar's Test P-Value : 0.01342         
                                          
            Sensitivity : 0.9043          
            Specificity : 0.9119          
         Pos Pred Value : 0.9524          
         Neg Pred Value : 0.8302          
             Prevalence : 0.6608          
         Detection Rate : 0.5975          
   Detection Prevalence : 0.6274          
      Balanced Accuracy : 0.9081          
                                          
       'Positive' Class : 0               
                                          

SVM radial

svm.radial <- svm(formula = y~., data = breast.cancer2, kernel='radial', degree=2, gamma=1, coef0=1, cost=C, cross=10, scale=FALSE)
summary(svm.radial)

Call:
svm(formula = y ~ ., data = breast.cancer2, kernel = "radial", degree = 2, gamma = 1, coef0 = 1, cost = C, cross = 10, 
    scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  10 
      gamma:  1 

Number of Support Vectors:  212

 ( 106 106 )


Number of Classes:  2 

Levels: 
 0 1

10-fold cross-validation on training data:

Total Accuracy: 90.15817 
Single Accuracies:
 92.85714 91.22807 87.7193 89.47368 92.98246 94.73684 87.7193 89.47368 84.21053 91.22807 

#Estadisticas: precisión y cobertura
table(predict(svm.radial), breast.cancer2$y, dnn=c("Prediction", "Actual"))

          Actual
Prediction   0   1
         0 340  36
         1  17 176

library(caret)
confusionMatrix(breast.cancer2$y, predict(svm.radial))

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 340  17
         1  36 176
                                          
               Accuracy : 0.9069          
                 95% CI : (0.8799, 0.9294)
    No Information Rate : 0.6608          
    P-Value [Acc > NIR] : < 2e-16         
                                          
                  Kappa : 0.7971          
 Mcnemar's Test P-Value : 0.01342         
                                          
            Sensitivity : 0.9043          
            Specificity : 0.9119          
         Pos Pred Value : 0.9524          
         Neg Pred Value : 0.8302          
             Prevalence : 0.6608          
         Detection Rate : 0.5975          
   Detection Prevalence : 0.6274          
      Balanced Accuracy : 0.9081          
                                          
       'Positive' Class : 0               
                                          

Se ve que es mejor el kernel LINEAL (grado), pero debemos mirar el accuracy y cobertura –> Gana CUADRATICO, parece.

ARBOLES DE DECISIÓN

Link: https://www.youtube.com/watch?v=mgh4KdbYHv0 Es un algoritmo de CLASIFICACIÓN supervisada. La variable dependiente predicha puede ser:

• Cuantitativa –> Árbol de REGRESIÓN
• Cualitativa –> Árbol de CLASIFICACIÓN

Cargar dataset

library(rpart) #Algoritmo para árbol de decisión
library(rpart.plot) #Plot de árbol de decisión
library("C50") #paquete que tiene el dataset de entrada y el algoritmo (rpart)
data(churn) #DATASET de ejemplo sobre CHURN (clientes que abandonan una empresa)

Preparar set de datos, dividiendo en entrenamiento y test

churn_juntos <- rbind(churnTest, churnTrain) #junta filas
head(churn_juntos, n=3L) #FEATURES disponibles (5000x20). Mostramos 3 filas

churn <- churn_juntos[, c(4,7,8,16,19,17,20)] #Usaremos solo algunas columnas (features)
names(churn) # nombres en ingles

[1] "international_plan"            "total_day_minutes"             "total_day_calls"               "total_intl_minutes"           
[5] "number_customer_service_calls" "total_intl_calls"              "churn"                        

names(churn) <- c("Tiene plan internacional", "Minutos al dia", "Llamadas al dia", "Minutos internacionales", "Reclamaciones", "Llamadas internacionales", "Cancelacion") #los renombro a nombres en espanhol
ind <- sample(2, nrow(churn), replace=TRUE, prob=c(0.6, 0.4)) #train (60%) y test (40%)
trainData <- churn[ind==1, ] #train
testData <- churn[ind==2, ] #test

Crear árbol de decisión

ArbolRpart <- rpart(Cancelacion ~ ., method="class", data=trainData) #Formula: variable dependiente (Cancelacion) depende todas las otras

Gráficos

print(ArbolRpart)                         

n= 2965 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 2965 425 no (0.14333895 0.85666105)  
   2) Minutos al dia>=253.9 255 121 yes (0.52549020 0.47450980)  
     4) Minutos al dia>=316.35 19   0 yes (1.00000000 0.00000000) *
     5) Minutos al dia< 316.35 236 115 no (0.48728814 0.51271186)  
      10) Minutos al dia>=265.4 150  65 yes (0.56666667 0.43333333) *
      11) Minutos al dia< 265.4 86  30 no (0.34883721 0.65116279) *
   3) Minutos al dia< 253.9 2710 291 no (0.10738007 0.89261993)  
     6) Reclamaciones>=3.5 221 108 no (0.48868778 0.51131222)  
      12) Minutos al dia< 160.3 88  13 yes (0.85227273 0.14772727) *
      13) Minutos al dia>=160.3 133  33 no (0.24812030 0.75187970) *
     7) Reclamaciones< 3.5 2489 183 no (0.07352350 0.92647650)  
      14) Tiene plan internacional=yes 226  80 no (0.35398230 0.64601770)  
        28) Minutos internacionales>=13.05 42   0 yes (1.00000000 0.00000000) *
        29) Minutos internacionales< 13.05 184  38 no (0.20652174 0.79347826)  
          58) Llamadas internacionales< 2.5 34   0 yes (1.00000000 0.00000000) *
          59) Llamadas internacionales>=2.5 150   4 no (0.02666667 0.97333333) *
      15) Tiene plan internacional=no 2263 103 no (0.04551480 0.95448520) *

rpart.plot(ArbolRpart,extra=4)  # extra=4:probabilidad de observaciones por clase

printcp(ArbolRpart)             # estadísticas de resultados

Classification tree:
rpart(formula = Cancelacion ~ ., data = trainData, method = "class")

Variables actually used in tree construction:
[1] Llamadas internacionales Minutos al dia           Minutos internacionales  Reclamaciones            Tiene plan internacional

Root node error: 425/2965 = 0.14334

n= 2965 

        CP nsplit rel error  xerror     xstd
1 0.058824      0   1.00000 1.00000 0.044896
2 0.049412      3   0.82353 0.85647 0.042046
3 0.030588      6   0.64471 0.67529 0.037883
4 0.010000      8   0.58353 0.62353 0.036551

plotcp(ArbolRpart)              # evolución del error a medida que se incrementan los nodos

# Podado del árbol
pArbolRpart<- prune(ArbolRpart, cp= ArbolRpart$cptable[which.min(ArbolRpart$cptable[,"xerror"]),"CP"])
pArbolRpart<- prune(ArbolRpart, cp= 0.011111)
printcp(pArbolRpart)

Classification tree:
rpart(formula = Cancelacion ~ ., data = trainData, method = "class")

Variables actually used in tree construction:
[1] Llamadas internacionales Minutos al dia           Minutos internacionales  Reclamaciones            Tiene plan internacional

Root node error: 425/2965 = 0.14334

n= 2965 

        CP nsplit rel error  xerror     xstd
1 0.058824      0   1.00000 1.00000 0.044896
2 0.049412      3   0.82353 0.85647 0.042046
3 0.030588      6   0.64471 0.67529 0.037883
4 0.010000      8   0.58353 0.62353 0.036551

Predice las cancelaciones en testData

# Validamos la capacidad de predicción del árbol con el fichero de validación
testPredRpart <- predict(ArbolRpart, newdata = testData, type = "class")
# Visualizamos una matriz de confusión
table(testPredRpart, testData$Cancelacion)

             
testPredRpart  yes   no
          yes  190   61
          no    92 1692

Estadística del modelo

# Calculamos el % de aciertos
sum(testPredRpart == testData$Cancelacion) / length(testData$Cancelacion)*100

[1] 92.48157

SVM (otro ejemplo)

Modelo SVM sobre datos de iris:

data(iris)
attach(iris)

The following objects are masked from iris (pos = 8):

    Petal.Length, Petal.Width, Sepal.Length, Sepal.Width, Species

The following objects are masked from iris (pos = 13):

    Petal.Length, Petal.Width, Sepal.Length, Sepal.Width, Species

The following objects are masked from iris (pos = 19):

    Petal.Length, Petal.Width, Sepal.Length, Sepal.Width, Species

## classification mode
# default with factor response:
model <- svm(Species ~ ., data = iris)
# alternatively the traditional interface:
x <- subset(iris, select = -Species)
y <- Species
model <- svm(x, y) 
print(model)

Call:
svm.default(x = x, y = y)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.25 

Number of Support Vectors:  51

summary(model)

Call:
svm.default(x = x, y = y)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.25 

Number of Support Vectors:  51

 ( 8 22 21 )


Number of Classes:  3 

Levels: 
 setosa versicolor virginica

Testeo con datos de entrenamiento (claro, el modelo está muy ajustado para ellos):

# test with train data
pred <- predict(model, x)
# (same as:)
#pred <- fitted(model)
# Check accuracy:
table(pred, y)

            y
pred         setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         2
  virginica       0          2        48

# compute decision values and probabilities:
pred <- predict(model, x, decision.values = TRUE)
attr(pred, "decision.values")[1:4,]

  setosa/versicolor setosa/virginica versicolor/virginica
1          1.196152         1.091757            0.6708810
2          1.064621         1.056185            0.8483518
3          1.180842         1.074542            0.6439798
4          1.110699         1.053012            0.6782041

# visualize (classes by color, SV by crosses):
distancias_escaladas <- cmdscale(dist(iris[,-5]))
plot(distancias_escaladas,
     col = as.integer(iris[,5]),
     pch = c("o","+")[1:150 %in% model$index + 1])

Modelo de REGRESION en 2 dimensiones:

## try regression mode on two dimensions
# create data
x <- seq(0.1, 5, by = 0.05)
y <- log(x) + rnorm(x, sd = 0.2)
# estimate model and predict input values
m   <- svm(x, y)
summary(m)

Call:
svm.default(x = x, y = y)


Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  1 
    epsilon:  0.1 


Number of Support Vectors:  79

new <- predict(m, x)
# visualize
plot(x, y)
points(x, log(x), col = 2)

points(x, new, col = 4)

## density-estimation

Creamos otras features con distribución normal y modelamos:

# create 2-dim. normal with rho=0:
X <- data.frame(a = rnorm(1000), b = rnorm(1000))
attach(X)

The following objects are masked from X (pos = 7):

    a, b

The following objects are masked from X (pos = 8):

    a, b

The following objects are masked from X (pos = 12):

    a, b

The following objects are masked from X (pos = 13):

    a, b

The following objects are masked from X (pos = 18):

    a, b

The following objects are masked from X (pos = 19):

    a, b

# MODELO:
#traditional way:
#m <- svm(X, gamma = 0.1)
# formula interface:
m <- svm(~., data = X, gamma = 0.1)
# Otra forma:X <- data.frame(a = rnorm(1000), b = rnorm(1000))
attach(X)

The following objects are masked from X (pos = 3):

    a, b

The following objects are masked from X (pos = 8):

    a, b

The following objects are masked from X (pos = 9):

    a, b

The following objects are masked from X (pos = 13):

    a, b

The following objects are masked from X (pos = 14):

    a, b

The following objects are masked from X (pos = 19):

    a, b

The following objects are masked from X (pos = 20):

    a, b

#m <- svm(~ a + b, gamma = 0.1)
summary(m)

Call:
svm(formula = ~., data = X, gamma = 0.1)


Parameters:
   SVM-Type:  one-classification 
 SVM-Kernel:  radial 
      gamma:  0.1 
         nu:  0.5 

Number of Support Vectors:  500




Number of Classes: 1

# test:
newdata <- data.frame(a = c(0, 4), b = c(0, 4))
predict(m, newdata)

    1     2 
 TRUE FALSE 

# visualize:
plot(X, col = 1:1000 %in% m$index + 1, xlim = c(-5,5), ylim=c(-5,5))
points(newdata, pch = "+", col = 2, cex = 5)

# weights: (example not particularly sensible)
i2 <- iris
levels(i2$Species)[3] <- "versicolor"
summary(i2$Species)

    setosa versicolor 
        50        100 

wts <- 100 / table(i2$Species)
wts

    setosa versicolor 
         2          1 

m <- svm(Species ~ ., data = i2, class.weights = wts)

SVM (ejemplo perros y gatos)

Link: https://www.youtube.com/watch?v=ImhjKyc88x0

library(e1071)
data(cats, package = 'MASS') # Features: sexo (sex), peso (bwt), peso_corazon (hwt)
#Datasets: train y test
ind <- sample(x=2, size = nrow(cats), replace = TRUE, prob = c(0.7,0.3)) #vector 1:144 de índices. Sus valores están en el rango entre 1:x; en este caso, los unos aparecerán un 70% y los doses un 30%.
testset <- cats[ind == 1,] #Coge las filas de CATS en aquellos indices que 'ind' toma valor 1
trainset <- cats[ind == 2,] #Coge las filas de CATS en aquellos indices que 'ind' toma valor 2
#Entrenamos el modelo
modelo_svm <- svm(Sex~., data=trainset, kernel = "radial")
prediccion <- predict(modelo_svm, newdata = testset[-1])
#Resultado y porcentaje de acierto
plot(modelo_svm, cats)

MC <- table(testset[,1], prediccion) #matriz de confusión
MC

   prediccion
     F  M
  F  7 23
  M  4 58

acierto <- ( sum(diag(MC)) )/ (sum(MC))
acierto

[1] 0.7065217

Otro ejemplo completo

Link: https://machinelearningmastery.com/how-to-estimate-model-accuracy-in-r-using-the-caret-package/

Data split:

# load the libraries
library(caret)
library(klaR)
# load the iris dataset
data(iris)
# define an 80%/20% train/test split of the dataset
split=0.80
trainIndex <- createDataPartition(iris$Species, p=split, list=FALSE)
data_train <- iris[ trainIndex,]
data_test <- iris[-trainIndex,]
# train a naive bayes model
model <- NaiveBayes(Species~., data=data_train)
# make predictions
x_test <- data_test[,1:4]
y_test <- data_test[,5]
predictions <- predict(model, x_test)
# summarize results
confusionMatrix(predictions$class, y_test)

Confusion Matrix and Statistics

            Reference
Prediction   setosa versicolor virginica
  setosa         10          0         0
  versicolor      0         10         2
  virginica       0          0         8

Overall Statistics
                                          
               Accuracy : 0.9333          
                 95% CI : (0.7793, 0.9918)
    No Information Rate : 0.3333          
    P-Value [Acc > NIR] : 8.747e-12       
                                          
                  Kappa : 0.9             
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: setosa Class: versicolor Class: virginica
Sensitivity                 1.0000            1.0000           0.8000
Specificity                 1.0000            0.9000           1.0000
Pos Pred Value              1.0000            0.8333           1.0000
Neg Pred Value              1.0000            1.0000           0.9091
Prevalence                  0.3333            0.3333           0.3333
Detection Rate              0.3333            0.3333           0.2667
Detection Prevalence        0.3333            0.4000           0.2667
Balanced Accuracy           1.0000            0.9500           0.9000

Bootstrap:

# load the library
library(caret)
# load the iris dataset
data(iris)
# define training control
train_control <- trainControl(method="boot", number=100)
# train the model
model <- train(Species~., data=iris, trControl=train_control, method="nb")
# summarize results
print(model)

Naive Bayes 

150 samples
  4 predictors
  3 classes: 'setosa', 'versicolor', 'virginica' 

No pre-processing
Resampling: Bootstrapped (100 reps) 
Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
Resampling results across tuning parameters:

  usekernel  Accuracy   Kappa    
  FALSE      0.9518209  0.9270189
   TRUE      0.9537661  0.9300203

Tuning parameter 'fL' was held constant at a value of 0
Tuning parameter 'adjust' was held constant at a value of 1
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were fL = 0, usekernel = TRUE and adjust = 1.

k-fold Cross Validation:

# load the library
library(caret)
# load the iris dataset
data(iris)
# define training control
train_control <- trainControl(method="cv", number=10)
# fix the parameters of the algorithm
grid <- expand.grid(.fL=c(0), .usekernel=c(FALSE))
# train the model
model <- train(Species~., data=iris, trControl=train_control, method="nb")#, tuneGrid=grid
# summarize results
print(model)

Naive Bayes 

150 samples
  4 predictors
  3 classes: 'setosa', 'versicolor', 'virginica' 

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 135, 135, 135, 135, 135, 135, ... 
Resampling results across tuning parameters:

  usekernel  Accuracy  Kappa
  FALSE      0.96      0.94 
   TRUE      0.96      0.94 

Tuning parameter 'fL' was held constant at a value of 0
Tuning parameter 'adjust' was held constant at a value of 1
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were fL = 0, usekernel = FALSE and adjust = 1.

Repeated k-fold Cross Validation:

# load the library
library(caret)
# load the iris dataset
data(iris)
# define training control
train_control <- trainControl(method="repeatedcv", number=10, repeats=3)
# train the model
model <- train(Species~., data=iris, trControl=train_control, method="nb")
# summarize results
print(model)

Naive Bayes 

150 samples
  4 predictors
  3 classes: 'setosa', 'versicolor', 'virginica' 

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 3 times) 
Summary of sample sizes: 135, 135, 135, 135, 135, 135, ... 
Resampling results across tuning parameters:

  usekernel  Accuracy   Kappa    
  FALSE      0.9555556  0.9333333
   TRUE      0.9600000  0.9400000

Tuning parameter 'fL' was held constant at a value of 0
Tuning parameter 'adjust' was held constant at a value of 1
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were fL = 0, usekernel = TRUE and adjust = 1.

Leave One Out Cross Validation:

# load the library
library(caret)
# load the iris dataset
data(iris)
# define training control
train_control <- trainControl(method="LOOCV")
# train the model
model <- train(Species~., data=iris, trControl=train_control, method="nb")
# summarize results
print(model)

Naive Bayes 

150 samples
  4 predictors
  3 classes: 'setosa', 'versicolor', 'virginica' 

No pre-processing
Resampling: Leave-One-Out Cross-Validation 
Summary of sample sizes: 149, 149, 149, 149, 149, 149, ... 
Resampling results across tuning parameters:

  usekernel  Accuracy   Kappa
  FALSE      0.9533333  0.93 
   TRUE      0.9600000  0.94 

Tuning parameter 'fL' was held constant at a value of 0
Tuning parameter 'adjust' was held constant at a value of 1
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were fL = 0, usekernel = TRUE and adjust = 1.

Otro ejemplo SVM vs RPART:

library(e1071)
library(rpart)
#DATOS:
data(Ozone, package="mlbench")
## Split en Datasets (train y test)
index     <- 1:nrow(Ozone)
tercio <- trunc(length(index)/3) #usamos un 30% para test
testindex <- sample(index, tercio, replace = FALSE, prob = NULL)
testset   <- na.omit(Ozone[testindex,]) #TEST Quita las filas que tengan ALGUN valor NA
trainset  <- na.omit(Ozone[-testindex,])#TRAIN Quita las filas que tengan ALGUN valor NA
#Columna TARGET
indice_target <- which( colnames(trainset) == "V4" )
## svm (MAQUINAS VECTOR SOPORTE) para REGRESIÓN NO LINEAL
svm.model <- svm(formula = V4 ~ ., data = trainset, cost = 1000, gamma = 0.0001)
svm.pred  <- predict(svm.model, testset[,-indice_target])
svm_error <- crossprod(svm.pred - testset[,indice_target]) / length(testindex)
svm_error

         [,1]
[1,] 10.00099

## rpart (ARBOLES DE REGRESION) para REGRESIÓN NO LINEAL
rpart.model <- rpart(formula = V4 ~ ., data = trainset)
rpart.pred  <- predict(rpart.model, testset[,-indice_target])
rpart_error <- crossprod(rpart.pred - testset[,indice_target]) / length(testindex)
rpart_error

         [,1]
[1,] 23.82726

ENSEMBLING

Link: https://www.datacamp.com/community/tutorials/ensemble-r-machine-learning

library("SuperLearner")
library(MASS) #Data
# Train and test sets
train <- Pima.tr
test <- Pima.te
head(train)

summary(train)

     npreg            glu              bp              skin            bmi             ped              age         type    
 Min.   : 0.00   Min.   : 56.0   Min.   : 38.00   Min.   : 7.00   Min.   :18.20   Min.   :0.0850   Min.   :21.00   No :132  
 1st Qu.: 1.00   1st Qu.:100.0   1st Qu.: 64.00   1st Qu.:20.75   1st Qu.:27.57   1st Qu.:0.2535   1st Qu.:23.00   Yes: 68  
 Median : 2.00   Median :120.5   Median : 70.00   Median :29.00   Median :32.80   Median :0.3725   Median :28.00            
 Mean   : 3.57   Mean   :124.0   Mean   : 71.26   Mean   :29.21   Mean   :32.31   Mean   :0.4608   Mean   :32.11            
 3rd Qu.: 6.00   3rd Qu.:144.0   3rd Qu.: 78.00   3rd Qu.:36.00   3rd Qu.:36.50   3rd Qu.:0.6160   3rd Qu.:39.25            
 Max.   :14.00   Max.   :199.0   Max.   :110.00   Max.   :99.00   Max.   :47.90   Max.   :2.2880   Max.   :63.00            

help(Pima.tr)

El campo target es type (que es un factor), que indica si el paciente tiene diabetes o no.

#Since the type column was a factor, R will encode it to 1 and 2, but this is not what you want: ideally, you would like to work with the type encoded as 0 and 1, which are "No" and "Yes", respectively. In the above code chunk, you subtract 1 from the whole set to get your 0-1 encoding. R will also encode this in the factor order.
y <- as.numeric(train[,8]) - 1
ytest <- as.numeric(test[,8]) - 1
#Predictores y respuestas
x <- data.frame(train[,1:7])
xtest <- data.frame(test[,1:7])
#MODELOS DISPONIBLES:
listWrappers()

All prediction algorithm wrappers in SuperLearner:

 [1] "SL.bartMachine"      "SL.bayesglm"         "SL.biglasso"         "SL.caret"            "SL.caret.rpart"     
 [6] "SL.cforest"          "SL.dbarts"           "SL.earth"            "SL.extraTrees"       "SL.gam"             
[11] "SL.gbm"              "SL.glm"              "SL.glm.interaction"  "SL.glmnet"           "SL.ipredbagg"       
[16] "SL.kernelKnn"        "SL.knn"              "SL.ksvm"             "SL.lda"              "SL.leekasso"        
[21] "SL.lm"               "SL.loess"            "SL.logreg"           "SL.mean"             "SL.nnet"            
[26] "SL.nnls"             "SL.polymars"         "SL.qda"              "SL.randomForest"     "SL.ranger"          
[31] "SL.ridge"            "SL.rpart"            "SL.rpartPrune"       "SL.speedglm"         "SL.speedlm"         
[36] "SL.step"             "SL.stepAIC"          "SL.step.forward"     "SL.step.interaction" "SL.svm"             
[41] "SL.template"         "SL.xgboost"         

All screening algorithm wrappers in SuperLearner:

[1] "All"
[1] "screen.corP"           "screen.corRank"        "screen.glmnet"         "screen.randomForest"   "screen.SIS"           
[6] "screen.template"       "screen.ttest"          "write.screen.template"

#Modelos, usando cross-validation (con 5 intentos):
set.seed(150)
model <- SuperLearner(y, 
                      x,
                      SL.library=list("SL.ranger",
                                      "SL.ksvm",
                                      "SL.ipredbagg",
                                      "SL.bayesglm"))
# Return the model
summary(model)

                  Length Class  Mode       
call                4    -none- call       
libraryNames        4    -none- character  
SL.library          2    -none- list       
SL.predict        200    -none- numeric    
coef                4    -none- numeric    
library.predict   800    -none- numeric    
Z                 800    -none- numeric    
cvRisk              4    -none- numeric    
family             11    family list       
fitLibrary          4    -none- list       
varNames            7    -none- character  
validRows          10    -none- list       
method              3    -none- list       
whichScreen         7    -none- logical    
control             2    -none- list       
cvControl           4    -none- list       
errorsInCVLibrary   4    -none- logical    
errorsInLibrary     4    -none- logical    
metaOptimizer       8    nnls   list       
env               115    -none- environment
times               3    -none- list       

Predicciones con SuperLearner: