R y Mineria de Datos

techfest.uc3m.es

job: null

knit: slidify::knit2slides

mode: selfcontained

highlighter: highlight.js

subtitle: techfest.uc3m.es

— .class #id ## Realizacion de las transparencias

• Transparencias realizadas con slidify, un paquete de R que permite mezclar código, resultados y texto, con conversion a html y pdf.

• Otros paquetes similares: knitr (con salida a html, pdf, LaTeX, …).

— .class #id ## Indice

• Mineri­a de datos
• Visualizacion
• Otras cuestiones

— .class #id ## Mineria de datos

• Construccion de modelos a partir de datos para hacer predicciones en el futuro
• Ejemplo de conjuntos de datos: prediccion de crecimiento de algas (a1) en los rios a partir de la estacion, size, compuestos quimicos, etc.

print(xtable(head(algae)), type="html")

	season	size	speed	mxPH	mnO2	Cl	NO3	NH4	oPO4	PO4	Chla	a1
1	winter	small	medium	8.00	9.80	60.80	6.24	578.00	105.00	170.00	50.00	0.00
2	spring	small	medium	8.35	8.00	57.75	1.29	370.00	428.75	558.75	1.30	1.40
3	autumn	small	medium	8.10	11.40	40.02	5.33	346.67	125.67	187.06	15.60	3.30
4	spring	small	medium	8.07	4.80	77.36	2.30	98.18	61.18	138.70	1.40	3.10
5	autumn	small	medium	8.06	9.00	55.35	10.42	233.70	58.22	97.58	10.50	9.20
6	winter	small	high	8.25	13.10	65.75	9.25	430.00	18.25	56.67	28.40	15.10

— .class #id ## Data Wrangling, Data Munging

• Transformar datos “sin procesar” en datos “limpios” que se puedan utilizar en minei­a de datos
• Es una parte de la mineri­a de datos que lleva un tiempo substancial
• Paso de formato en texto libre a matriz de datos (mediante expresiones regulares, …)
• Eliminar NAs, cambios de formato de matriz ancho a largo (o al contrario), agregacion (compactacion de datos mediante la media, mediana u otras, …)
• Paquetes utiles en R: plyr, dplyr, tydir, reshape2, …

— .class #id ## Ejemplo: Data Wrangling, Data Munging

• Vamos a compactar los datos, dejando solo tres columnas: Season, PO4 y la salida a1, pero de season solo apareceran sus cuatro valores (primavera, verano, otono e invierno), y por cada valor, en PO4 y a1 aparecera la media para cada estacion:

library(plyr)
algaeCompactado = ddply(algae, "season", summarize, PO4Mean = mean(PO4), a1Mean = mean(a1)) 
print(xtable(head(algaeCompactado)), type="html")

	season	PO4Mean	a1Mean
1	autumn	102.95	16.70
2	spring	172.15	12.68
3	summer	138.45	15.32
4	winter	158.73	16.66

— .class #id ## Visualizacion de dependencia entre PO4 y la variable de salida (a1)

Visualizacion de la dependencia de la concentracion de algas (a1) frente a la de fosfato (PO4)

plot(a1 ~ PO4, algae)

— .class #id ## A probar: es relevante la variable NH4? y mnO2

# Primero nos traemos los datos de entrenamiento de las algas
load(url("http://www.dcc.fc.up.pt/~ltorgo/DataMiningWithR/DataSets/algae.RData"))
# Despues le quitamos los datos con valores faltantes
algae = na.omit(algae[,1:12])
# Y ahora, haced un plot que muestre si la variable NH4 
# es relevante para predecir la salida a1 (idem para mnO2)

— .class #id ## Visualizacion de dependencias

entre todos los atributos y la variable de salida (a1)

library("corrplot")
correlaciones = cor(algae[,sapply(algae, is.numeric)])
corrplot(correlaciones, method = "ellipse")

— .smallcode12 ## Creacion de un modelo lineal

modelo = lm(a1 ~ . , algae)
summary(modelo)

## 
## Call:
## lm(formula = a1 ~ ., data = algae)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -37.67 -10.69  -2.69   6.14  67.71 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  28.109289  29.710249    0.95   0.3454   
## seasonspring -1.085417   4.306550   -0.25   0.8013   
## seasonsummer -0.917019   4.080509   -0.22   0.8225   
## seasonwinter  1.820816   3.983706    0.46   0.6482   
## sizemedium    2.604712   3.828343    0.68   0.4972   
## sizesmall     8.957460   4.236788    2.11   0.0360 * 
## speedlow      1.699436   4.927303    0.34   0.7306   
## speedmedium  -2.240937   3.407871   -0.66   0.5117   
## mxPH         -1.064961   3.490481   -0.31   0.7607   
## mnO2          0.754383   0.708709    1.06   0.2887   
## Cl           -0.032487   0.033091   -0.98   0.3276   
## NO3          -1.595442   0.550627   -2.90   0.0043 **
## NH4           0.001783   0.000993    1.80   0.0744 . 
## oPO4         -0.015050   0.039573   -0.38   0.7042   
## PO4          -0.040151   0.030821   -1.30   0.1945   
## Chla         -0.107572   0.082105   -1.31   0.1919   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.3 on 168 degrees of freedom
## Multiple R-squared:  0.338,  Adjusted R-squared:  0.279 
## F-statistic: 5.72 on 15 and 168 DF,  p-value: 2.09e-09

— .class #id ## La formula de R

Supongamos que queremos crear un modelo en terminos de PO4 unicamente. La formula de R permite especificar que atributos se usan:

• “.” : todos
• PO4: unicamente PO4
• PO4+NH4: unicamente estos dos atributos

modelo = lm(a1 ~ PO4, algae)

— .class #id ## La formula de R

plot(a1 ~ PO4, algae)
abline(modelo, col="red")

— .class #id ## La formula de R

Si queremos crear el modelo en terminos de solo dos variables:

modelo2 =  lm(a1 ~ PO4+NH4, algae)

— .class #id ## Evaluacion del modelo en un conjunto de test

Error cuadratico medio: \(RMSE = \sqrt{\frac{\sum_{i=1}^{i=N} (y_i-modelo(x_i))^2}{N}}\)

rmse = function(error) sqrt(mean(error^2))
modelo = lm(a1 ~ ., algae)
predicciones_entrenamiento = predict(modelo, algae)
predicciones_test = predict(modelo, algaeTest)
error_entrenamiento = rmse(predicciones_entrenamiento-algae$a1)
error_test = rmse(predicciones_test-algaeTest$a1)
(errores_lineal = c(entrenamiento=error_entrenamiento, test=error_test))

## entrenamiento          test 
##         16.49         15.10

— .class #id ## Modelos que se pueden aprender

• Arboles de decision (rpart)
• Redes de neuronas (nnet, neuralnet, RSNNS, amore)
• Maquinas de vectores de soporte (e1071, kernlab)
• Random Forests (Random Forests)
• Gradient Boosting (gbm)
• …
• El paquete caret incluye a todos estos y muchos mas

— .class #id ## Aprendizaje de modelos con gbm

library(gbm)
modelo = gbm(a1 ~ ., data=algae, n.trees=200, interaction.depth=3)

## Distribution not specified, assuming gaussian ...

predicciones_entrenamiento = predict(modelo, algae, n.trees=200)
predicciones_test = predict(modelo, algaeTest, n.trees=200)
error_entrenamiento = rmse(predicciones_entrenamiento-algae$a1)
error_test = rmse(predicciones_test-algaeTest$a1)
(errores_gbm=c(entrenamiento=error_entrenamiento, test=error_test))

## entrenamiento          test 
##         18.50         16.11

— .class #id ## Aprendizaje de modelos con caret

caret es un paquete R que incluye la mayor parte de los algoritmos de aprendizaje disponibles en R, y permite automatizar muchas tareas de mineria de datos, tales como el ajuste de parametros, pero también la selección de atributos, etc.

http://caret.r-forge.r-project.org/

— .class #id ## Aprendizaje de modelos con MLR

MLR es otro paquete R que automatiza muchas tareas de mineria de datos (seleccion de atributos, validacion cruzada, validacion de modelos, etc.).

https://github.com/berndbischl/mlr

— .class #id ## Ajuste de parametros con caret y gbm

library(caret)

## Loading required package: ggplot2
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:survival':
## 
##     cluster

fitControl <- trainControl(method = "cv", ## 10-fold CV
                           number = 10)

modelo <- train(a1 ~ ., data = algae,
                 method = "gbm",
                 trControl = fitControl,
                 verbose = FALSE)

— .smallcode12 ## Ajuste de parametros con caret y gbm

modelo

## Stochastic Gradient Boosting 
## 
## 184 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## 
## Summary of sample sizes: 165, 165, 166, 165, 165, 165, ... 
## 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  RMSE  Rsquared  RMSE SD  Rsquared SD
##   1                   50      15    0.5       5        0.2        
##   1                  100      15    0.5       5        0.2        
##   1                  150      15    0.5       5        0.2        
##   2                   50      15    0.5       5        0.2        
##   2                  100      15    0.5       5        0.2        
##   2                  150      15    0.5       5        0.2        
##   3                   50      15    0.5       5        0.2        
##   3                  100      15    0.5       5        0.2        
##   3                  150      16    0.5       5        0.2        
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## RMSE was used to select the optimal model using  the smallest value.
## The final values used for the model were n.trees = 50, interaction.depth = 2 and shrinkage = 0.1.

## Stochastic Gradient Boosting 
## 
## 184 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## 
## Summary of sample sizes: 165, 165, 166, 165, 165, 165, ... 
## 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  RMSE  Rsquared  RMSE SD  Rsquared SD
##   1                   50      15    0.5       5        0.2        
##   1                  100      15    0.5       5        0.2        
##   1                  150      15    0.5       5        0.2        
##   2                   50      15    0.5       5        0.2        
##   2                  100      15    0.5       5        0.2        
##   2                  150      15    0.5       5        0.2        
##   3                   50      15    0.5       5        0.2        
##   3                  100      15    0.5       5        0.2        
##   3                  150      16    0.5       5        0.2        
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## RMSE was used to select the optimal model using  the smallest value.
## The final values used for the model were n.trees = 50, interaction.depth = 2 and shrinkage = 0.1.

— .class #id ## Ajuste de parametros con caret y gbm

print(xtable(modelo$results[order(modelo$results$RMSE),]),type="html")

	shrinkage	interaction.depth	n.trees	RMSE	Rsquared	RMSESD	RsquaredSD
4	0.10	2	50.00	14.97	0.49	5.15	0.22
1	0.10	1	50.00	14.99	0.49	5.05	0.21
2	0.10	1	100.00	15.06	0.49	5.13	0.20
7	0.10	3	50.00	15.15	0.47	4.81	0.20
5	0.10	2	100.00	15.17	0.49	5.11	0.21
6	0.10	2	150.00	15.32	0.49	5.04	0.21
3	0.10	1	150.00	15.36	0.47	5.17	0.21
8	0.10	3	100.00	15.45	0.45	5.32	0.24
9	0.10	3	150.00	15.63	0.46	5.35	0.24

— .class #id ## Ajuste de parametros con caret y gbm

trellis.par.set(caretTheme())
plot(modelo)

— .class #id ## Ajuste de parametros con caret y gbm

trellis.par.set(caretTheme())
plot(modelo, plotType = "level",
     scales = list(x = list(rot = 90)))

— .class #id ## Ajuste de parametros con caret y gbm

predicciones_entrenamiento = predict(modelo, algae)
predicciones_test = predict(modelo, algaeTest)
error_entrenamiento = rmse(predicciones_entrenamiento-algae$a1)
error_test = rmse(predicciones_test-algaeTest$a1)
(errores_gbm_ajustado=c(
  entrenamiento=error_entrenamiento, 
  test=error_test))

## entrenamiento          test 
##         12.07         13.81

— .class #id ## Ajuste de parametros con caret y gbm

print(xtable(rbind(errores_lineal, errores_gbm, errores_gbm_ajustado)), type="html")

	entrenamiento	test
errores_lineal	16.49	15.10
errores_gbm	18.50	16.11
errores_gbm_ajustado	12.07	13.81

— .class #id ## A probar (1): ajustar parametros de maquina de vectores de soporte

1. Primero: los datos de entrenamiento ya los teneis disponibles. Descargaros ahora los datos de test:

load(url("http://www.dcc.fc.up.pt/~ltorgo/DataMiningWithR/DataSets/testAlgae.RData"))
load(url("http://www.dcc.fc.up.pt/~ltorgo/DataMiningWithR/DataSets/algaeSols.RData"))
algaeTest = cbind(a1=algae.sols$a1, test.algae)
algaeTest = na.omit(algaeTest[1:12])

— .class #id ## A probar (2): ajustar parametros de maquina de vectores de soporte

2. Segundo: Instalar el paquete de maquinas de vectores de soporte, entrenar un modelo con algunos valores de los parametros y calcular los errores de test:

install.packages("kernlab")
library("kernlab")
rmse = function(error) sqrt(mean(error^2))
modeloSVM = ksvm(a1~.,data=algae, C=5, kpar=list(sigma=0.9))
predicciones_test = predict(modeloSVM, algaeTest)
print(rmse(predicciones_test-algaeTest$a1))

— .class #id ## A probar (3): ajustar parametros de maquina de vectores de soporte

3. Tercero: Ajustar los parametros con caret, de la siguiente manera:

install.packages("caret")
library("caret")

fitControl <- trainControl(method = "cv", ## 10-fold CV
                           number = 10)

modeloSVMCaret <- train(a1 ~ ., data = algae,
                 method = "svmRadial",
                 trControl = fitControl,
                 verbose = FALSE)

modeloSVMCaret
predicciones_test = predict(modeloSVMCaret, algaeTest)
print(rmse(predicciones_test-algaeTest$a1))

— .class #id ## Visualizacion (plots) en R

Varios sistemas:

• Base: comenzar con un plot vacio e ir anadiendo resultados y anotaciones
• lattice: Utiles cuando queremos ver como una variable cambia en terminos de las demas
• ggplot: Lenguaje de los graficos. Construye plots por capas (layers)

— .class #id ## Visualizacion (plots) en R

plot(a1 ~ PO4, algae)
predicciones = predict(modelo, algae)
points(predicciones ~ PO4, algae, col="red")
legend("topright", pch=1, col=c("black", "red"), legend=c("Real", "Modelo"))

— .class #id ## Visualizacion en base R: Histograma

hist(algae$PO4)

— .class #id ## Visualizacion en base R: boxplot

boxplot(a1 ~ size, algae)

— .class #id ## Visualizacion interactiva en base R: manipulate

library(manipulate)
manipulate(hist(algae[,variable], breaks=valor), 
           valor=slider(2,20),
           variable=picker("PO4","Cl","a1"))

— .class #id ## Visualizacion en Lattice: variables condicionales

miplot = xyplot(a1 ~ PO4|size, algae, auto.key=TRUE)
print(miplot)

— .class #id ## Visualizacion en Lattice: variables condicionales

miplot = xyplot(a1 ~ PO4|size*speed, algae, auto.key=TRUE)
print(miplot)

— .class #id ## Visualizacion en Lattice

algaep = cbind(algae, pred=predicciones)
miplot = xyplot(a1+pred ~ PO4|size*speed, algaep, auto.key=TRUE)
print(miplot)

— .class #id ## Visualizacion en Lattice: grupos

algaeg = algae
algaeg$grande = ifelse(algae$a1>median(algae$a1), "GRANDE", "PEQUENO")
miplot = xyplot(Cl ~ PO4, algaeg, groups=algaeg$grande, auto.key=TRUE)
print(miplot)

— .class #id ## Visualizacion en ggplot

library(ggplot2)
miplot=qplot(PO4, Cl, data=algaeg, col=algaeg$grande)
print(miplot)

— .class #id ## Visualizacion en ggplot

Construye un plot elemento a elemento, basado en el “lenguaje de los graficos”. Los principales elementos son:

• Aesthetics: que papel juega cada variable en el plot? (posicion de los puntos, color de los puntos, forma de los puntos, altura de una barra, …)
• Geoms: que objetos geometricos se muestran? (barras, puntos, lineas, …)
• Statistics: que objetos estadisticos se pueden usar para resumir los resultados (lineas de tendencia, …)

— .class #id ## Visualizacion en ggplot

miplot = ggplot(algaeg, aes(x=PO4, y=a1, col=grande))
print(miplot + geom_point())

— .class #id ## Visualizacion en ggplot

Anadimos un estadi­stico (tendencia):

print(miplot + geom_point() + stat_smooth())

## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.

— .class #id ## Graficos interactivos: Shiny, rCharts

Se pueden hacer interfaces sencillos basados en javascript con Shiny. Estas aplicaciones pueden ser subidas posteriormente a Rpubs:

• Ejemplo Shiny 1
• Ejemplo Shiny 2
• Ejemplo rCharts

— .class #id ## Caracteristicas avanzadas de R

• Paralelismo y multicore: foreach, doMC (multicore), snowfall (clusters), (http://cran.r-project.org/web/views/HighPerformanceComputing.html) [hig performance en general: http://cran.r-project.org/web/views/HighPerformanceComputing.html]
• Graficos interactivos: rCharts, shiny, ggvis, plotly
• Big data: RHadoop, SparkR, h2o (Oxdata)