-
Notifications
You must be signed in to change notification settings - Fork 0
/
Practica_2.Rmd
1074 lines (828 loc) · 54 KB
/
Practica_2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: 'Tipología y ciclo de vida de los datos - Práctica 2 - Limpieza y análisis de datos'
author: "Autores: Alberto Quesada / Fernando Sevilla"
date: "Mayo 2020"
output:
html_document:
highlight: default
number_sections: yes
theme: cosmo
toc: yes
toc_depth: 3
includes:
in_header: 75.584-PEC-header.html
pdf_document:
highlight: zenburn
toc: yes
word_document: default
---
<style>
body {
text-align: justify}
</style>
```{r setup, include=FALSE}
knitr::opts_chunk$set(eval=T, echo=T)
```
******
# Detalles de la actividad
******
## Decripción
A lo largo de esta práctica se ha realizado un caso práctico sobre el estudio de los datasets Red/White Wine Quality con el fin de aprender a identificar los datos relevantes de un proyecto analítico, así como utilizar las técnicas apropiadas para la integración, limpieza, validación y análisis de los datos.
## Competencias
A continuación, se muestran las competencias desarrolladas en esta práctica del Máster de Data Science.
* Capacidad de analizar un problema en el nivel de abstracción adecuado a cada situación y aplicar las habilidades y conocimientos adquiridos para abordarlo y resolverlo.
* Capacidad para aplicar las técnicas específicas de tratamiento de datos (integración, transformación, limpieza y validación) para su posterior análisis.
## Objetivos
Como se especifica en el enunciado de esta práctica, los objetivos son los siguientes:
* Aprender a aplicar los conocimientos adquiridos y su capacidad de resolución de problemas en entornos nuevos o poco conocidos dentro de contextos más amplios o multidisciplinares.
* Saber identificar los datos relevantes y los tratamientos necesarios (integración, limpieza y validación) para llevar a cabo un proyecto analítico.
* Aprender a analizar los datos adecuadamente para abordar la información contenida en los datos.
* Identificar la mejor representación de los resultados para aportar conclusiones sobre el problema planteado en el proceso analítico.
* Actuar con los principios éticos y legales relacionados con la manipulación de datos en Tipología y ciclo de vida de los datos Práctica 2 pág 2 función del ámbito de aplicación.
* Desarrollar las habilidades de aprendizaje que les permitan continuar estudiando de un modo que tendrá que ser en gran medida autodirigido o autónomo.
* Desarrollar la capacidad de búsqueda, gestión y uso de información y recursos en el ámbito de la ciencia de datos.
## Nota: Propiedad intelectual
> A menudo es inevitable, al producir una obra multimedia, hacer uso de recursos creados por terceras personas. Es por lo tanto comprensible hacerlo en el marco de una práctica de los estudios de Informática, Multimedia y Telecomunicación de la UOC, siempre y cuando esto se documente claramente y no suponga plagio en la práctica.
> Por lo tanto, si se hicera uso de recursos ajenos, se presentará un documento detallando y especificando el nombre de cada recurso, su autor, el lugar dónde se obtuvo y su estatus legal: si la obra está protegida por el copyright o se acoge a alguna otra licencia de uso (Creative Commons, licencia GNU, GPL ...).
Como autores de la práctica nos aseguraremos que la licencia no impide específicamente su uso en el marco de la práctica. En caso de no encontrar la información correspondiente tendrá que asumir que la obra está protegida por copyright.
******
# Resolución
******
## Introducción
### Descripción del dataset
Como ya hemos comentado anteriormente, para el desarrollo de la práctica se ha analizado los datasets Red/White Wine Quality que podemos encontrar en UCI Machine Learning Repository. En estos ficheros cada registro representa un ensayo de calidad fisicoquímico sobre un determinado vino en base a una serie de atributos. Ambos ficheros presentan la misma configuración y el mismo número de atributos y registros como se muestra a continuación.
| Formato | CSV |
|---------|-----|
| Separador | ; |
| Header | True |
| Número de registros | 4898 |
| Número de atributos | 12 |
A continuación, se muestra una descripción de los atributos de los ficheros en estudio.
| Nombre | Descripción | Tipo |
|----------------------|----------------------------------------|----------|
| fixed acidity | Acidez fija | Numérico |
| volatile acidity | Acidez volátil | Numérico |
| citric acid | Ácido cítrico | Numérico |
| residual sugar | Azúcar residual | Numérico |
| chlorides | Cloruros | Numérico |
| free sulfur dioxide | Dióxido de azufre libre | Numérico |
| total sulfur dioxide | Dióxido de azufre total | Numérico |
| density | Densidad | Numérico |
| pH | PH | Numérico |
| sulphates | Sulfatos | Numérico |
| alcohol | Porcentaje de alcohol en volumen | Numérico |
| quality | Valor de calidad resultante del ensayo | Numérico |
El primer paso antes de realizar el análisis será unir ambos datasets para crear un único conjunto de datos.
### Importancia y objetivos del análisis
El objetivo es realizar un análisis completo del dataset de modo que trabajemos todas las fases de la analítica avanzada.
```{r, out.width = "50%", fig.align="center", fig.retina = 2, echo=FALSE, fig.cap="Tipos de análisis", chunck1}
library(knitr)
include_graphics("Imagen_Analisis.png")
```
* Análisis descriptivo. Nos permite responder a la pregunta ¿Qué ha pasado? En otras palabras, nos ofrece un resumen de los datos históricos.
* Análisis de diagnóstico. Nos permite responder a la pregunta ¿Por qué ha pasado? En esta fase se identifican las relaciones entre las variables en estudio.
* Análisis predictivo. Nos permite responder a la pregunta ¿Qué pasará? En esta fase se generan los modelos y se validan.
* Análisis prescriptivo. Nos permite responder a la pregunta ¿Cómo podemos hacer que pase? Se aplica el modelo a una situación real para buscar puntos óptimos y responder a preguntas estratégicas para el diseño de nuevos productos.
## Análisis del dataset
### Instalación de librerías necesarias para exportar
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck2}
library(readr) # Providing a fast and friendly way to read rectangular data ('csv').
library(knitr) # A General-Purpose Package for Dynamic Report Generation in R.
library(tinytex) # Helper Functions to Install and Maintain 'TeX Live', and Compile 'LaTeX' Documents.
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax. Build complex HTML or 'LaTeX' tables
library(tidyverse) # Collection of R packages designed for data science.
library(scales) # Graphical scales map data to aesthetics, and provide methods for automatically determining breaks and labels for axes
library(dplyr) # tool for working with data frame like objects, both in memory and out of memory.
library(tables) # Computes and displays complex tables of summary statistics
library(ggplot2) # A system for 'declaratively' creating graphics, based on "The Grammar of Graphics"
library(ggExtra) # Package for adding marginal histograms to ggplot2
library(magrittr) # Provides a mechanism for chaining commands with a new forward-pipe operator, %>%.
library(corrplot) # The corrplot package is a graphical display of a correlation matrix, confidence interval....
library(grid) # grid adds an nx by ny rectangular grid to an existing plot.
library(gridExtra) # Provides a number of user-level functions to work with "grid" graphics, notably to arrange multiple grid-based plots on a page, and draw tables.
library(visdat) # Create preliminary exploratory data visualisations of an entire dataset to identify problems or unexpected features using 'ggplot2'.
library(moments) # Compute skewness of a univariate distribution.
library(caret) # Pre-processing transformation (centering, scaling etc.) can be estimated from the training data and applied to any data set with the same variables.
library(nortest) # Performs the Lilliefors (Kolmogorov-Smirnov) test for the composite hypothesis of normality,
library(car) # VIF Regression: A Fast Regression Algorithm For Large Data
library(e1071) #Functions for latent class analysis, short time Fourier transform, fuzzy clustering, support vector machines, shortest path computation, bagged clustering, naive Bayes classifier.
library(party)# A computational toolbox for recursive partitioning.
library(rpart) # Recursive partitioning for classification, regression and survival trees
library(rpart.plot) # Plot an rpart model
library(ROCR) #ROC graphs, sensitivity/specificity curves, lift charts, and precision/recall plots
library(randomForest) #Classification and regression based on a forest of trees using random inputs, based on Breiman (2001)
library(cluster) # Methods for Cluster analysis.
```
### Análisis Descriptivo
En primer lugar vamos a proceder a cargar ambos archivos. En este caso se han extraído directaemnte de la url, de modo que no tengan que pasar por nuestro dispositivo, optimizando el proceso. No obstante, se adjunta la línea de código comentada para cargar los archivos desde la carpeta raíz, por si hubiera alguna modificación de la url.
```{r message=FALSE, warning=FALSE, chunck3}
# Se carga el primer conjunto de datos
white_url <- "http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-white.csv"
whitewine <- read.csv(white_url, header = TRUE, sep = ";")
# whitewine <- read.csv(winequality-white.csv, header = TRUE, sep = ";")
# Se carga el segundo conjunto de datos
red_url <- "http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"
redwine <- read.csv(white_url, header = TRUE, sep = ";")
#redwine <- read.csv(winequality-red.csv, header = TRUE, sep = ";")
#Añadimos una nueva columna con el color del vino antes de unir los dos conjuntos de datos.
whitewine$color<-"white"
redwine$color<-"red"
```
Unión de los conjuntos de datos
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck4}
# Mediante la función rbind se unen (join) verticalmente los conjuntos de datos
wine <- rbind(redwine, whitewine)
#Una vez unidos creamos una nueva columna para futuros modelos que evaluarán la calidad
quality.factor <- ifelse(wine$quality >= 6, "Buena", "Mala")
wine <- data.frame(wine, quality.factor)
table(wine$quality.factor)
```
Empezaremos haciendo un breve análisis de los datos ya que nos interesa tener una idea general de los datos que disponemos. Por ello, primero calcularemos las dimensiones de nuestra base de datos y analizaremos qué tipos de atributos tenemos.
Para empezar, calculamos las dimensiones mediante la función dim(). Obtenemos que disponemos de 9796 registros (filas) y 13 variables (columnas).
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck5}
# Dimensiones de la base de datos
dim(wine)
# Estructura de los datos
str(wine)
```
Una vez cargados los datos, vamos a verificar la transmisión de la información.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck6}
# Hacemos una primera visualización
head(wine)
# Comprobamos el tipo de variable en los datos
sapply(wine,class)
# Mostramos un resumen de los datos originales
summary(wine)
```
Factorizamos la variable "color" para futuros análisis.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck7}
wine$color<-as.factor(wine$color)
wine$quality.factor<-as.factor(wine$quality.factor)
```
Se eliminan los registros duplicados si los hubiese
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck8}
wine <- wine[!duplicated(wine), ]
dim(wine)
```
Comprobamos sie existen valores NA
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck9}
vis_miss(wine)
```
Todos los valores están presentes. No existen valores NA.
Como podemos ver, todas las variables han sido identificadas con el tipo correspondiente y no se aprecian anomalías a simple vista en los datos.
### Análisis Descriptivo visual
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck10}
attach(wine)
#Iniciamos el cuadrante para el gráfico
par(mfrow=c(1,2))
# Representamos histogramas y boxplot para las variables numéricas
hist(fixed.acidity , main="fixed.acidity histogram")
boxplot(fixed.acidity, main="fixed.acidity boxplot")
hist(volatile.acidity , main="volatile.acidity histogram")
boxplot(volatile.acidity , main="volatile.acidity boxplot")
hist(citric.acid , main="citric.acid histogram")
boxplot(citric.acid , main="citric.acid boxplot")
hist(residual.sugar , main="residual.sugar histogram")
boxplot(residual.sugar , main="residual.sugar boxplot")
hist(chlorides , main="chlorides histogram")
boxplot(chlorides , main="chlorides boxplot")
hist(free.sulfur.dioxide , main="free.sulfur.dioxide histogram")
boxplot(free.sulfur.dioxide , main="free.sulfur.dioxide boxplot")
hist(total.sulfur.dioxide , main="total.sulfur.dioxide histogram")
boxplot(total.sulfur.dioxide , main="total.sulfur.dioxide boxplot")
hist(density , main="density histogram")
boxplot(density , main="density boxplot")
hist(sulphates , main="sulphates histogram")
boxplot(sulphates , main="sulphates boxplot")
hist(alcohol , main="alcohol histogram")
boxplot(alcohol , main="alcohol boxplot")
hist(quality , main="quality histogram")
boxplot(quality , main="quality boxplot")
# Representamos gráfico de barras de las variables nominales
barplot(table(color), main="Color histogram")
barplot(table(quality.factor), main="quality.factor histogram")
detach(wine)
```
Todas las variables presentan outliers.
* Quality tiene la mayoría de los valores concentrados en las categorías 5, 6 y 7. Sólo una pequeña proporción está en las categorías [3, 4] y [8, 9] y ninguna en las categorías [1, 2, 10].
* Fixed acidity, volatile.acidity y citric.acid tienen valores atípicos. Si se eliminan esos valores atípicos, la distribución de las variables puede considerarse simétrica.
* Residual.sugar residual tiene una distribución positivamente sesgada; incluso después de eliminar los valores atípicos, la distribución seguirá siendo sesgada.
* Algunas de las variables, por ejemplo, free.sulphur.dioxide y density,, tienen unos pocos valores atípicos, pero éstos son muy diferentes del resto.
* Alcohol tiene una distribución irregular pero no tiene valores atípicos pronunciados.
### Comprobación de la asimetría
Comprobamos a continuación la asimetría de las variables verificar si los datos se distribuyen normalmente o no.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck11}
attach(wine)
skewness(fixed.acidity)
skewness(volatile.acidity)
skewness(citric.acid)
skewness(residual.sugar)
skewness(chlorides)
skewness(free.sulfur.dioxide)
skewness(total.sulfur.dioxide)
skewness(density)
skewness(pH)
skewness(sulphates)
skewness(alcohol)
skewness(quality)
detach(wine)
```
Regla de la asimetría:
* Si la asimetría es 0, los datos son perfectamente simétricos.
* Si la asimetría es menor que -1 o mayor que 1, la distribución es altamente sesgada.
* Si la asimetría está entre -0,5 y 0,5, la distribución es aproximadamente simétrica.
Fuente:*[GoodData.Documentation](https://help.gooddata.com/doc/en/reporting-and-dashboards/maql-analytical-query-language/maql-expression-reference/aggregation-functions/statistical-functions/predictive-statistical-use-cases/normality-testing-skewness-and-kurtosis)
### Transformación de las variables
Se aplica transformación de Boxcox y volvemos a comprobar la asimetría.
Fuente: *[Box-Cox Transformations](http://onlinestatbook.com/2/transformations/box-cox.html)
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck12}
wine_procesado <- preProcess(wine[,1:11], c("BoxCox", "center", "scale"))
new_wine <- data.frame(trans = predict(wine_procesado, wine))
colnames(new_wine)
```
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck13}
attach(new_wine)
skewness(trans.fixed.acidity)
skewness(trans.volatile.acidity)
skewness(trans.citric.acid)
skewness(trans.residual.sugar)
skewness(trans.chlorides)
skewness(trans.free.sulfur.dioxide)
skewness(trans.total.sulfur.dioxide)
skewness(trans.density)
skewness(trans.quality)
skewness(trans.pH)
skewness(trans.sulphates)
skewness(trans.alcohol)
detach(new_wine)
```
Si el coeficiente de asimetría no se aproxima a 0, entonces su conjunto de datos no se distribuye normalmente.
Ahora la mayoría las distribuciones son aproximadamente simétricas (la asimetría está entre -0,5 y 0,5). Las únicas variables que muestran aún una asimestría son:
* trans.citric.acid
* trans.density
Se a va a comprobr mediante gráfico de cuantiles teóricos (Gráficos Q-Q) la normalidad de las dos variables que muestran aín coeficientes de simetría alejados de 0.
Para el resto de variables no es necesario comprobarlo ya que sus coeficientes se distribuyen alrededor de 0 y por tanto se distribuyen normalmente.
Se realiza una comprobación de la normalidad para la variable __trans.citric.acid__.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck14}
# Para comprobar la normalidad de la variable trans.citric.acid se va a usar:
# 1. Gráfico de cuantiles teóricos (Gráficos Q-Q): Consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta.
# 2. Un histograma
# 3. Diagrama de cajas
par(mfrow = c(1, 3))
hist(new_wine$trans.citric.acid, las=1, main="Distribución normal", font.main=1, ylab="Frecuencia", xlab = "trans.citric.acid")
qqnorm(new_wine$trans.citric.acid, las=1, pch=18, main="Simetria", font.main=1, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
qqline(new_wine$trans.citric.acid)
boxplot(new_wine$trans.citric.acid, las=1, main="Valores extremos", font.main=1, xlab="Datos de distribution normal", horizontal=F)
```
No se trata de una distribución normal ya que presenta una asimetría en ambas colas, aunque la cola a la derecha de la media es más larga que la izquierda.
Para asegurarnos completamente se realiza la pruena de Lilliefors (Kolmogorov-Smirnov) para muestras grandes (n>50).
El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para contrastar la normalidad.
Es la alternativa al test de Shapiro-Wilk cuando el número de observaciones es mayor de 50. La función lillie.test() del paquete nortest permite aplicarlo.
__Planteamos la Hipótesis:__
$\ H_0:$ La muestra proviene de una distribución normal.
$\ H_1:$ La muestra no proviene de una distribución normal.
__Nivel de significancia:__
El nivel de significancia que se trabajará es de 0.05. $\alpha=0.05$
__Criterio de decisón:__
Si $p<\alpha$ se rechaza $\ H_0:$
Si $p>\alpha$ no se rechaza $\ H_0:$
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck15}
lillie.test(x=new_wine$trans.citric.acid)
```
Dado que $p-value<\alpha$ se rechaza $\ H_0:$. La muestra no proviene de una distribución normal. Pero, puesto que el tamaño de la muestra es mayor a 30, según el teorema del límite central, asumimos la normalidad.
Se realiza una comprobación de la normalidad para la variable __trans.density__.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck16}
# Para comprobar la normalidad de la variable trans.citric.acid se va a usar:
# 1. Gráfico de cuantiles teóricos (Gráficos Q-Q): Consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta.
# 2. Un histograma
# 3. Diagrama de cajas
par(mfrow = c(1, 3))
hist(new_wine$trans.density, las=1, main="Distribución normal", font.main=1, ylab="Frecuencia", xlab = "trans.density")
qqnorm(new_wine$trans.density, las=1, pch=18, main="Simetria", font.main=1, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
qqline(new_wine$trans.density)
boxplot(new_wine$trans.density, las=1, main="Valores extremos", font.main=1, xlab="Datos de distribution normal", horizontal=F)
```
Presenta una asimetría en en la cola izquierda.
Para asegurarnos completamente se realiza la pruena de Lilliefors (Kolmogorov-Smirnov) para muestras grandes (n>50).
El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para contrastar la normalidad.
Es la alternativa al test de Shapiro-Wilk cuando el número de observaciones es mayor de 50. La función lillie.test() del paquete nortest permite aplicarlo.
__Planteamos la Hipótesis:__
$\ H_0:$ La muestra proviene de una distribución normal.
$\ H_1:$ La muestra no proviene de una distribución normal.
__Nivel de significancia:__
El nivel de significancia que se trabajará es de 0.05. $\alpha=0.05$
__Criterio de decisón:__
Si $p<\alpha$ se rechaza $\ H_0:$
Si $p>\alpha$ no se rechaza $\ H_0:$
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck17}
lillie.test(x=new_wine$trans.density)
```
Dado que $p-value<\alpha$ se rechaza $\ H_0:$. La muestra no proviene de una distribución normal. Pero, puesto que el tamaño de la muestra es mayor a 30, según el teorema del límite central, asumimos la normalidad.
### Eliminación de los Outliers
Se analiza cada variable independientemente para determinar qué valores atípicos, si los hay, deben ser eliminados.
__trans.fixed.acidity__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck18}
boxplot(trans.fixed.acidity ~ trans.quality , data=new_wine,
xlab="Quality", ylab="Fixed Acidity")
summary(new_wine$trans.fixed.acidity)
quantile(new_wine$trans.fixed.acidity, c(.999))
new_wine$trans.quality[new_wine$trans.fixed.acidity > 3.4]
```
Se observan bastantes valores atípicos, especialmente para la calidad de 6. El mayor es tres veces más grande que el valor del tercer cuartil.
El 99,9% de las observaciones tienen una acidez fija tras la normalización de los datos de 3,34 o menos.
Hay 6 outliers. De nuevo, la mayoría de ellas están clasificadas con una calidad de 6. Esto no es sorprendente, porque muchas observaciones en el conjunto de datos están clasificadas como calidad 6.
__trans.volatile.acidity__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck19}
boxplot(trans.volatile.acidity ~ trans.quality , data=new_wine,
xlab="Quality", ylab="volatile acidity")
summary(new_wine$trans.volatile.acidity)
quantile(new_wine$trans.volatile.acidity, c(.999))
new_wine$trans.quality[new_wine$trans.volatile.acidity > 3.7]
```
Esta vez la mayoría de los outliers tienen la calidad 4. Esto es interesante porque no hay muchas observaciones con la calidad 4. La acidez puede ser un factor clave para predecir los vinos malos, por ello puede que estos no sean en realidad valores atípicos.
__trans.citric.acid__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck20}
boxplot(trans.citric.acid ~ trans.quality , data=new_wine,
xlab="Quality", ylab="citric acid")
summary(new_wine$trans.citric.acid)
quantile(new_wine$trans.citric.acid, c(.999))
new_wine$trans.quality[new_wine$trans.citric.acid > 5.5]
```
La mayoría de outliers están clasificados como 6 en calidad. El boxplot es interesante porque tanto el vino de calidad 3 como el de calidad 9 no tienen outliers y todos los demás rangos sí. Esto podría ser sólo una casualidad. También podría ser que la variación en el ácido cítrico está de alguna manera relacionada con los vinos de nivel medio.
__trans.residual.sugar__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck21}
boxplot(trans.residual.sugar ~ trans.quality , data=new_wine,
xlab="Quality", ylab="residual sugar")
summary(new_wine$trans.residual.sugar)
quantile(new_wine$trans.residual.sugar, c(.999))
new_wine$trans.quality[new_wine$trans.residual.sugar > 2.0]
```
EL boxplot no muestra valores atípicos. La calidad 6 muestra 6 valores altos.
__trans.chlorides__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck22}
boxplot(trans.chlorides ~ trans.quality , data=new_wine,
xlab="Quality", ylab="chlorides")
summary(new_wine$trans.chlorides)
quantile(new_wine$trans.chlorides, c(.999))
new_wine$trans.quality[new_wine$trans.chlorides > 3.8]
```
Los vinos con calidad 9 están muy bien distribuidos y no hay outliers. Además, la media es más baja para este grupo que para cualquier otro. Es posible que un nivel muy bajo de cloruros sea indicativo de un muy buen vino. También podría ser debido al número extremadamente pequeño de vinos calidad 9.
__trans.free.sulfur.dioxide__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck23}
boxplot(trans.free.sulfur.dioxide ~ trans.quality , data=new_wine,
xlab="Quality", ylab="free sulfur dioxide")
summary(new_wine$trans.free.sulfur.dioxide)
quantile(new_wine$trans.free.sulfur.dioxide, c(.999))
new_wine$trans.quality[new_wine$trans.free.sulfur.dioxide > 3.9]
```
Todos las calidades tienen al menos un atípico y los medios son bastante similares con la excepción de 4. Parece que a medida que la calidad mejora, el rango de variación tiende a disminuir.
__trans.total.sulfur.dioxide__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck24}
boxplot(trans.total.sulfur.dioxide ~ trans.quality , data=new_wine,
xlab="Quality", ylab="total sulfur dioxide")
summary(new_wine$trans.total.sulfur.dioxide)
quantile(new_wine$trans.total.sulfur.dioxide, c(.999))
new_wine$trans.quality[new_wine$trans.total.sulfur.dioxide > 3.5]
```
La variación tiende a disminuir con la mejora de la calidad del vino.
__trans.density__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck25}
boxplot(trans.density ~ trans.quality , data=new_wine,
xlab="Quality", ylab="trans density")
summary(new_wine$trans.density)
quantile(new_wine$trans.density, c(.999))
new_wine$trans.quality[new_wine$trans.density > 3.0]
```
La calidad 6 muestra un claro outlier. Dada la cantidad de vinos con calidad 6, parece poco probable que una densidad de este nivel sea precisa o representativa del grupo.
__trans.pH__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck26}
boxplot(trans.pH ~ trans.quality , data=new_wine,
xlab="Quality", ylab="trans pH")
summary(new_wine$trans.pH)
quantile(new_wine$trans.pH, c(.999))
new_wine$trans.quality[new_wine$trans.pH > 3.2]
```
Se observa una mayor dispersión en la calidad 3, pero no hay valores atípicos, mientras que la calidad 9 muestra una menor dispersión.
__trans.sulphates__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck27}
boxplot(trans.sulphates ~ trans.quality , data=new_wine,
xlab="Quality", ylab="trans sulphates")
summary(new_wine$trans.sulphates)
quantile(new_wine$trans.sulphates, c(.999))
new_wine$trans.quality[new_wine$trans.sulphates > 3.0]
```
Las medias de las calidades estaán todas en el mismo nivel. Es difícil diferenciar, a partir de las medias, si los sulfatos pueden llegar a influenciar en la calidad.
__trans.alcohol__
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck28}
boxplot(trans.alcohol ~ trans.quality , data=new_wine,
xlab="Quality", ylab="trans alcohol")
summary(new_wine$trans.alcohol)
quantile(new_wine$trans.alcohol, c(.999))
new_wine$trans.quality[new_wine$trans.alcohol > 2.21]
new_wine$trans.quality[new_wine$trans.alcohol < -2.8]
```
Parece haber una relación entre el alcohol y la calidad.
Procedemos a eliminar los valores outliers encontardos en los gráficos boxplot anteriores.
```{r echo=TRUE, message=FALSE, warning=FALSE,chunck29}
# Mostramos los valores identificados como outliers
attach(new_wine)
f_a_out<-boxplot.stats(trans.fixed.acidity)$out
v_a_out<-boxplot.stats(trans.volatile.acidity)$out
c_a_out<-boxplot.stats(trans.citric.acid)$out
r_s_out<-boxplot.stats(trans.residual.sugar)$out
ch_out<-boxplot.stats(trans.chlorides)$out
f_s_d_out<-boxplot.stats(trans.free.sulfur.dioxide)$out
t_s_d_out<-boxplot.stats(trans.total.sulfur.dioxide)$out
d_t_out<-boxplot.stats(trans.density)$out
ph_out<-boxplot.stats(trans.pH )$out
sh_out<-boxplot.stats(trans.sulphates)$out
ac_out<-boxplot.stats(trans.alcohol)$out
qa_out<-boxplot.stats(trans.quality)$out
detach(new_wine)
```
Las variables trans.alcohol, trans.residual.sugar y trans.color no presentan outliers.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck30}
# Eliminamos las instancias con outlier
new_wine<-new_wine[-which(new_wine$trans.fixed.acidity %in% f_a_out),]
new_wine<-new_wine[-which(new_wine$trans.volatile.acidity %in% v_a_out),]
new_wine<-new_wine[-which(new_wine$trans.citric.acid %in% c_a_out),]
new_wine<-new_wine[-which(new_wine$trans.chlorides %in% ch_out),]
new_wine<-new_wine[-which(new_wine$trans.free.sulfur.dioxide %in% f_s_d_out),]
new_wine<-new_wine[-which(new_wine$trans.total.sulfur.dioxide %in% t_s_d_out),]
new_wine<-new_wine[-which(new_wine$trans.density %in% d_t_out),]
new_wine<-new_wine[-which(new_wine$trans.pH %in% ph_out),]
new_wine<-new_wine[-which(new_wine$trans.sulphates %in% sh_out),]
new_wine<-new_wine[-which(new_wine$trans.quality %in% qa_out),]
```
### Análisis de diagnóstico
Tras haber eliminado los valores atípicos, revisamos la correlación entre las distintas variables. El objetivo es centrar el análisis descriptivo en aquellas variables que presentan correlación.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck31}
#Se omite la variable color para pocer representar la matriz de correlación
corrplot(cor(new_wine[1:12]))
```
Las correlaciones positivas se muestran en azul y las negativas en rojo. La intensidad del color y el tamaño del círculo son proporcionales a los coeficientes de correlación. En la parte derecha del correlograma, el color de la leyenda muestra los coeficientes de correlación y los colores correspondientes.
Observando los resultados se verifica que:
1. La calidad del vino no se encuentra correlacionada con las variables "citric.acid", "free.sulfur.dioxide", los "sulphates" y "color".
2. Es coherente la alta correlación positiva entre free.sulfur.dioxide y total.sulfur.dioxide. A mayor cantidad de sulfuro más sulfuro libre hay.
3. Es coherente la alta correlación negativa entre fixed.acidity y PH. Cuanto mayor es el PH en la escala de ácido-base más básico es el compuesto.Por lo tanto a valores altos de PH menor es la acidez.
4. Es coherente la alta correlación negativa entre el azucar y el alcohol.El alcohol es el resultado de la oxidación del azucar, por tanto si hay menos azucar el grado de alcohol aumenta.
A continuación se estudian visaulamente las posibles correlaciones entre:
### Variables residual.sugar y density
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck32}
# Calculamos la correlación y la mostramos por pantalla
RS_D_cor = cor(new_wine$trans.residual.sugar,new_wine$trans.density)
print(RS_D_cor)
# Generamos el scatter plot y añadimos la linea de regresión para una mejor visualización
plot(new_wine$trans.residual.sugar,new_wine$trans.density, main="residual.sugar vs density",xlab="residual.sugar", ylab="density", pch=19)
abline(lm(new_wine$trans.density~new_wine$trans.residual.sugar), col="red") # regression line (y~x)
```
Como podemos apreciar en la gráfica existe una fuerte relación lineal entre ambas variables y todos los puntos se ajustan bastante bien a la linea. El resultado de la correlación, 0.75, refuerza lo comentado anteriormente.
### Variables density y alcochol
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck33}
# Calculamos la correlación y la mostramos por pantalla
D_A_cor = cor(new_wine$trans.density,new_wine$trans.alcohol)
print(D_A_cor)
# Generamos el scatter plot y añadimos la linea de regresión para una mejor visualización
plot(new_wine$trans.density,new_wine$trans.alcohol, main="alcochol vs density",xlab="alcochol", ylab="density", pch=19)
abline(lm(new_wine$trans.alcohol~new_wine$trans.density), col="red") # regression line (y~x)
```
Al igual que antes, existe una fuerte relación lineal entre ambas variables y los puntos se ajustan a la linea. El resultado de la correlación, -0.80, refuerza lo comentado anteriormente.
### Todas las variables y quality
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck34}
# Calculamos la correlación y la mostramos por pantalla
par(mfrow = c(1,2))
for (i in c(1:11)) {
plot(new_wine[, i], jitter(new_wine[, "trans.quality"]), xlab = names(new_wine)[i],
ylab = "trans.quality", col = "firebrick", cex = 0.8, cex.lab = 1.3)
abline(lm(new_wine[, "trans.quality"] ~ new_wine[ ,i]), lty = 2, lwd = 2)
}
par(mfrow = c(1, 1))
```
Las variables trans.alcohol, trans.sulphates y trans.pH parecen tener una correlación positiva con trans.quality, mientras que ltrans.volatile.acidity, trans.chlorides y trans.total.sulfur.dioxidel parecen tener una relación negativa con la calidad.
### Contraste de hipótesis para la diferencia de medias
Dado que nuestro dataset contiene los valores para los vinos blancos y tintos, vamos a realizar un contraste de hipótesis para la diferencia de medias (trans.citric.acid) de las dos muestras. Planteamos la sigueinte hipótesis:
¿Podemos aceptar con un nivel de confianza del 95% que los vinos blancos (color) tienen una acidez (trans.citric.acid) que supera a la acidez de los vinos tintos (color)?.
Planteamos la Hipótesis nula y alternativa
$$
\left\{
\begin{array}{ll}
H_{0}: & \mu_{si}=\mu_{no}\\
H_{1}: & \mu_{si}>\mu_{no}
\end{array}
\right.\\
\ Hipótesis\ unilateral\ a \ la \ dereccha
$$
EL método que se aplica es un contraste de hipótesis para dos muestras independientes. El test es paramétrico, asumiendo el teorema del límite central. Es un test unilateral, con varianzas desconocidas, pero iguales, ya que el test de varianzas no rechaza la hipótesis nula de varianzas iguales.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck35}
#Test de igualdad de varianzas
#H0: F (ratio de varianzas) = 1
#H1: F diferente de 1
var.test(new_wine$trans.citric.acid[new_wine$trans.color=="white"], new_wine$trans.citric.acid[new_wine$trans.color=="red"],
alternative = "greater")
```
A continuación se lleva a cabo el contraste de hipótesis:
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck36}
#Se define una función myttest.
#x1, x2: muestras
#CL: confidence level (nivel de confianza)
#equalvar TRUE/FALSE
#alternative: "bilateral", "less" (x1 less than x2), "greater" (x1 greater than x2).
myttest <- function( x1, x2, CL=0.95,equalvar=TRUE, alternative="bilateral" ){ # z test
mean1<-mean(x1)
n1<-length(x1)
sd1<-sd(x1)
mean2<-mean(x2)
n2<-length(x2)
sd2<-sd(x2)
if (equalvar==TRUE){
s <-sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2 )/(n1+n2-2) )
Sb <- s*sqrt(1/n1 + 1/n2)
df<-n1+n2-2
}
else{ #equalvar==FALSE
Sb <- sqrt( sd1^2/n1 + sd2^2/n2 )
denom <- ( (sd1^2/n1)^2/(n1-1) + (sd2^2/n2)^2/(n2-2))
df <- ((sd1^2/n1 + sd2^2/n2)^2) / denom
}
alfa <- (1-CL)
t<- (mean1-mean2) / Sb
if (alternative=="bilateral"){
tcritical <- qt( alfa/2, df, lower.tail=FALSE ) #two sided
pvalue<-pt( abs(t), df, lower.tail=FALSE )*2 #two sided
}
else if (alternative=="less"){
tcritical <- qt( alfa, df, lower.tail=TRUE )
pvalue<-pt( t, df, lower.tail=TRUE )
}
else{ #(alternative=="greater")
tcritical <- qt( alfa, df, lower.tail=FALSE )
pvalue<-pt( t, df, lower.tail=FALSE )
}
#Guardamos en resultado en un data frame
info<-data.frame(t,df,tcritical,pvalue)
info %>% kable() %>% kable_styling()
return (info)
}
```
Se visualizan tanto el estadístico como el valor p:
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck37}
info<-myttest( new_wine$trans.citric.acid[new_wine$trans.color=="white"],
new_wine$trans.citric.acid[new_wine$trans.color=="red"],
alternative="less",
equalvar="TRUE")
info
```
Dado que hemos aplicado una función propia comprobamos el resultado con la fución original de R.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck38}
#Comparación con t.test de R
t.test( new_wine$trans.citric.acid[new_wine$trans.color=="white"], new_wine$trans.citric.acid[new_wine$trans.color=="red"], alternative = "less",
var.equal = TRUE)
```
El valor p es 0.5, mayor que el nivel de significación es 0.05 (95%). **No podemos rechazar la hipótesis nula** de igualdad de acidez media entre los vinos blancos y tintos. En otras palabras, no podemos asumir que la media de acidez de vinos blancos sea mayor a la de vinos tintos.
### Contraste de hipótesis para la relación entre variables
* Hipótesis nula: H<sub>0</sub>: La variable trans.color y la variable trans.quality son independientes.
* Hipótesis alternativa: H<sub>1</sub>: La variable trans.color y la variable trans.quality no son independientes.
Realizad los cálculos manualmente del test. Para ello, se recomienda construir una función my.chisq que realice estos cálculos y que podáis reusar más adelante.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck39}
# Realizamos el test con la función chisq.test para comprobar los resultados.
chisq.test(new_wine$trans.color,new_wine$trans.quality)
```
El valor crítico es mayor que el chi-cuadrado, por tanto, no podemos rechazar la hipótesis nula de que las variables son independientes.
Del mismo modo, el p valor es mayor al nivel de confianza 0.05 por tanto también debemos aceptar la hipótesis nula y podemos afirmar que no hay una dependencia entre el color y si la calidad del vino con un 95% de confianza.
### Análisis predictivo
Vamos a aplicar un modelo de regresión lineal multivariable. El objetivo consistirá en predecir la calidad del vino en función de las otras variables.
Antes de elegir el mejor modelo lineal vamos a dividir el conjunto de datos en dos:
1. Set de entrenamiento: Un subconjunto para entrenar un modelo (80% de los datos)
2. Set de prueba:Un subconjunto para probar el modelo entrenado (20% de lod datos)
```{r, out.width = "50%", fig.align="center", fig.retina = 2, echo=FALSE, fig.cap="División del conjunto de datos", chunck40}
library(knitr)
include_graphics("División Datos.png")
```
Esta división se ha llevado a cabo para cumplir con los siguientes condiciones:
* Que sea lo suficientemente grande como para generar resultados significativos desde el punto de vista estadístico.
* Que sea representativo de todo el conjunto de datos.
Se divide el conjunto de datos (80% entrenamiento - 20% test)
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck41}
set.seed(100)
#new_wine$trans.color<-as.numeric(new_wine$trans.color)
trainingRowIndex <- sample(1:nrow(new_wine), 0.8*nrow(new_wine))
new_wine_train <- new_wine[trainingRowIndex, ]
new_wine_test <- new_wine[-trainingRowIndex, ]
```
__Regresión lineal Múltiple.__
Realizamos el primer modelo
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck42}
# Incluimos todas las variables
modelo_1 <- lm(trans.quality ~ . , new_wine_train)
summary(modelo_1)
```
Todas las variables presentan significación muy alta, esto puede ser consecuencia de **multicolinealidad** entre las variables. Para detectarla y minimizarla vamos a seguir los siguientes pasos:
1. Análisis de la multicolinealidad
2. Modelo con variables que no muestran colinealidad
3. Análisis de la multicolinealidad
4. Modelo con variables más significativas
4. Predicción de valores
5. Matriz de confusión
**Multicolinealidad**: Se dice que existe multicolinealidad entre las variables explicativas cuando existe algún tipo de dependencia lineal entre ellas, o lo que es lo mismo, si existe una fuerte correlación entre las mismas. La correlación no solamente se refiere a las distintas variables dos a dos, sino a cualquier de ellas con cualquier grupo de las restantes. Por esta razón no es suficiente (aunque sí necesaria) que en la matriz de correlaciones bivariadas haya correlaciones altas.
El principal inconveniente de la multicolinealidad consiste en que se incrementan la varianza de los coeficientes de regresión estimados hasta el punto que resulta prácticamente imposible establecer su significación estadística, ya que como se sabe, el valor de t para un determinado coeficiente de regresión es el valor de dicho coeficiente dividido por su desviación tipica. Si este es grande, el valor de t será bajo y no llegara a la significación.
Hay varios procedimientos para detectar la multicolinealidad entre los predictores. El primero de ellos, basado en la correlación múltiple de un determinado regresor con los restantes se denomina Tolerancia de dicho regresor. Su valor es $\ 1-R^2_i$. Siendo $\ R^2_i$ la correlación múltiple al cuadrado de dicho regresor con los restantes.
Para que haya multicolinealidad dicha correlación ha de ser alta, o lo que es lo mismo la tolerancia baja.
Otro índice relacionado con éste y que nos da una idea del grado de aumento de la varianza se denomina Factor de Inflación de la Varianza, y es precisamente el recíproco de la tolerancia. Su valor es: $\ VIF=1/(1-R^2_y)$. Para que no haya multicolinealidad el denominador tiene que valer cerca de la unidad, por tanto un poco más de 1 el valor de VIF. **Cuanto mayor sea de este valor mayor multicolinealidad habrá y por tanto más inestable es el modelo lineal**; en otras palabras, es conveniente eliminar del modelo las variables con un factor VIF grande.
Si todos los VIF son 1, no hay multicolinealidad, pero si algunos VIF son mayores que 1, los predictores están correlacionados. Cuando un VIF es > 5, el coeficiente de regresión para ese término no se estima adecuadamente.
Fuentes:
*[REGRESIÓN LINEAL MÚLTIPLE](https://personal.us.es/vararey/adatos2/multiple.pdf)
*[Multicolinealidad](https://www.uv.es/uriel/material/multicolinealidad3.pdf)
Analizamos la multicolinealidad mediante el análisis del valor VIF
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck43}
vif(modelo_1)
```
Se observa multicolinealidad en la variable trans.density y se procede a eliminarla del modelo.
Ajustamos el modelo de nuevo
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck44}
# excluímos la variable trans.density
modelo_2 <- lm(trans.quality ~ . -trans.density, new_wine_train)
summary(modelo_2)
```
Analizamos la multicolinealidad mediante el análisis del valor VIF
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck45}
vif(modelo_2)
```
Los valores VIF son menores de 5 para todas las variables, en otras palabras la colinealidad es muy baja con el modelo 2.
Ajustamos un modelo nuevo con las variables que presentan significación:
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck46}
modelo_3 <- lm(trans.quality ~ trans.volatile.acidity +trans.free.sulfur.dioxide +trans.alcohol , new_wine_train)
summary(modelo_3)
```
Los valores VIF son menores de 1.1 para todas las variables, en otras palabras la no hay colinealidad en el modelo 3.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck47}
vif(modelo_3)
```
__Árbol de decición.__
Al igual que para el modelo realizado anteriormente, para la futura evaluación del árbol de decisión, es necesario dividir el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba.
Esta vez tras la división se llevará a cabo una normalización de los conjuntos de datos para evitar que alguna variable con valores más altos pueda influenciar a las otras.
No hay ninguna proporción fijada con respecto al número relativo de componentes de cada subconjunto, pero la más utilizada acostumbra a ser 2/3 para el conjunto de entrenamiento y 1/3, para el conjunto de prueba.
La variable por la que clasificaremos es el campo de si el vino es de alta o baja calidad, que está en la última columna.
Se divide el conjunto de datos en conjuntos de entrenamiento y prueba. La relación entrenamiento / prueba es normalmente de 70/30 ((2*6486/3=4.324). Queremos asegurarnos de que cada división representa el conjunto de datos en términos de las proporciones de la variable objetivo. La variable objetivo aquí es “Coste”
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck48}
set.seed(1)
new_wine_arbol<-new_wine[,1:12]
train <- new_wine_arbol[1:4324,]
test <- new_wine_arbol[4325:6486,]
# Filas y Columnas
n = nrow(new_wine_arbol)
p = ncol(new_wine_arbol)
```
Normalizamos el set de entrenamiento para que obtener un rango entre 0-1
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck49}
normalizado_train = function(x) (x - min(x))/(max(x) - min(x))
train.norm = data.frame(apply(train[,-p], 2, normalizado_train),
quality = train[,p])
summary(train.norm)
```
Normalizamos el set de prueba para que obtener un rango entre 0-1
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck50}
normalizado_test = function(x) (x - min(x))/(max(x) - min(x))
test.norm = data.frame(apply(test[,-p], 2, normalizado_test),
quality = test[,p])
summary(test.norm)
```
Construímos el árbol de decisión
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck51}
arbol = rpart(quality~., data=train.norm)
rpart.plot(arbol)
```
__Modelo modelo no supervisado y basado en el concepto de distancia__
Vamos a aplicar un modelo no supervisado. Para consegirlo no usaremos la columna quality.factor, que es la variable que se quiere predecir. Por lo tanto, intentaremos encontrar agrupaciones usando el resto de atributos que caracterizan a cada vino.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck52}
new_wine_cluster <- new_wine[,c(1:12)]
```
Como inicialmente no conocemos el número óptimo de clústers, probamos con varios valores
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck53}
d <- daisy(new_wine_cluster)
resultados <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
fit <- kmeans(new_wine_cluster, i)
y_cluster <- fit$cluster
sk <- silhouette(y_cluster, d)
resultados[i] <- mean(sk[,3])
}
```
Mostramos en un gráfica los valores de las siluetas media de cada prueba para comprobar que número de clústers es el mejor
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck54}
plot(2:10,resultados[2:10],type="o",col="blue",pch=0,xlab="Numbero de clusters",ylab="Silueta")
```
Aunque el valor esperado es k=2, dado que el conjunto original tiene 2 factores de calidada (Bueno y Malo), el mejor valor que se obtiene es k=4.
Otro forma de evaluar cual es el mejor número de clústers es considerar el mejor modelo, aquel que ofrece la menor suma de los cuadrados de las distancias de los puntos de cada grupo con respecto a su centro (withinss), con la mayor separación entre centros de grupos (betweenss). Como se puede comprobar es una idea conceptualmente similar a la silueta. Una manera común de hacer la selección del número de clústers consiste en aplicar el método elbow (codo), que no es más que la selección del número de clústers en base a la inspección de la gráfica que se obtiene al iterar con el mismo conjunto de datos para distintos valores del número de clústers. Se seleccionará el valor que se encuentra en el “codo” de la curva
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck55}
resultados <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
fit <- kmeans(new_wine_cluster, i)
resultados[i] <- fit$tot.withinss
}
plot(2:10,resultados[2:10],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="tot.tot.withinss")
```
En este caso el número óptimo de clústers no está tan claro dado que la curva no llega a estabilizarse.
También se puede usar la función kmeansruns del paquete fpc que ejecutá el algoritmo kmeans con un conjunto de valores y selecciona el valor del número de clústers que mejor funcione de acuerdo a dos criterios la silueta media (“asw”) y Calinski-Harabasz (“ch”).
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck56}
library(fpc)
#x$duración_corregida<-as.numeric(x$duración_corregida)
fit_ch <- kmeansruns(new_wine_cluster, krange = 1:10, criterion = "ch")
fit_asw <- kmeansruns(new_wine_cluster, krange = 1:10, criterion = "asw")
```
Podemos comprobar el valor con el que se ha obtenido el mejor resultado y también mostrar el resultado obtenido para todos los valores de k usando ambos criterios
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck57}
fit_ch$bestk
fit_asw$bestk
```
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck58}
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio Calinski-Harabasz")
```
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck59}
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio silueta media")
```
Los resultados son muy parecidos a los que hemos obtenido anteriormente. Con el criterio de la silueta media se obtienen dos clústers y con el Calinski-Harabasz se obtiene un peperfil parecido que muestra 2 clúster.
Como se ha comprobado, conocer el número óptimo de clústers no es un problema fácil. Tampoco lo es la evaluación de los modelos de agregación.
Como en el caso que estudiamos sabemos que los datos pueden ser agrupados en 2 factores de calidad, vamos a ver como se ha comportado el kmeans en el caso de pedirle 2 clústers. Para eso comparamos visualmente los campos dos a dos, con el valor real que sabemos está almacenado en el campo “quality.factor” del dataset original.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck60}
new_wine_cluster3 <- kmeans(new_wine_cluster, 2)
par(mfrow = c(1, 2))
# trans.residual.sugar y trans.density
plot(new_wine_cluster[c(4,8)], col=new_wine_cluster3$cluster)
# trans.residual.sugar y trans.alcohol
plot(new_wine_cluster[c(4,11)], col=new_wine_cluster3$cluster)
```
El atributo trans.residual.sugar parece hacer un mejor trabajo para dividir los vinos entre los dos factores de calidad.
### Análisis prescriptivo
__Regresión Lineal__
Predecimos valores utilizando el set de prueba para el modelo lineal.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck61}
# Se redondea el resultado para eliminar los decimales en la predicción
modelo_3_Pred = round(predict(modelo_3, newdata=new_wine_test))
head(modelo_3_Pred)
```
Creamos la matriz de confusión para comprobar la precisión del modelo. Se enfrentan los valores predichos por el modelo con los valores reales del set de prueba.
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck62}
confusionMatrix(table(modelo_3_Pred, new_wine_test$trans.quality))
```
Nuestro modelo tiene una precisión del 52.23%
Diagnosticamos el modelo a continuación
```{r echo=TRUE, message=FALSE, warning=FALSE, chunck63}
#Revisamos los residuos para asegurar que están distribuidos normalmente
hist(modelo_3$residuals)
#Comprobamos la media de los residuos para asegurarnos que ésta se encuentra cerca de 0
mean(modelo_3$residuals)
#Visualizamos los diagramas
layout(matrix(c(1,2,3,4),2,2))
plot(modelo_3)
```