1 Carga y depuración de los datos
1.1 Importación y aproximación a los datos
Antes de nada voy a cargar los datos y a echarles un ojo, además incorporo las librerías que voy a necesitar, para posteriormente elegir las variables.
source("FuncionesAida.R")
library(readxl)
library(questionr)
library(psych)
library(car)
library(corrplot)
library(ggplot2)
library(gridExtra)
library(dplyr)
library(caret)
library(lmSupport)
library(glmnet)
library(pROC)
DatosEleccionesEspaña <- read_excel("DatosEleccionesEspaña.xlsx")
datos<- as.data.frame(DatosEleccionesEspaña)| Name | CodigoProvincia | CCAA | Population | TotalCensus | AbstentionPtge | AbstencionAlta | Izda_Pct | Dcha_Pct | Otros_Pct | Izquierda | Derecha | Age_0-4_Ptge | Age_under19_Ptge | Age_19_65_pct | Age_over65_pct | WomanPopulationPtge | ForeignersPtge | SameComAutonPtge | SameComAutonDiffProvPtge | DifComAutonPtge | UnemployLess25_Ptge | Unemploy25_40_Ptge | UnemployMore40_Ptge | AgricultureUnemploymentPtge | IndustryUnemploymentPtge | ConstructionUnemploymentPtge | ServicesUnemploymentPtge | totalEmpresas | Industria | Construccion | ComercTTEHosteleria | Servicios | ActividadPpal | inmuebles | Pob2010 | SUPERFICIE | Densidad | PobChange_pct | PersonasInmueble | Explotaciones |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Abadía | 10 | Extremadura | 336 | 282 | 20,21 | 0 | 60,44 | 35,55 | 1,78 | 1 | 0 | 3,87 | 18,16 | 55,06 | 26,79 | 44,05 | 0,89 | 79,76 | 0,30 | 19,34 | 2,38 | 54,76 | 42,86 | 4,76 | 9,52 | 11,90 | 73,81 | 15 | 0 | 0 | 0 | 0 | Otro | 216 | 326 | 4507,56 | MuyBaja | 3,07 | 1,56 | 28 |
| Abertura | 10 | Extremadura | 429 | 364 | 25,27 | 0 | 54,78 | 44,12 | 0,37 | 1 | 0 | 1,63 | 13,05 | 56,64 | 30,30 | 50,12 | 1,63 | 90,91 | 2,80 | 7,23 | 16,22 | 32,43 | 51,35 | 8,11 | 8,11 | 10,81 | 67,57 | 11 | 0 | 0 | 0 | 0 | Otro | 382 | 459 | 6270,76 | MuyBaja | -6,54 | 1,12 | 67 |
| Acebo | 10 | Extremadura | 569 | 569 | 27,24 | 0 | 44,20 | 53,14 | 0,97 | 0 | 1 | 1,23 | 9,14 | 54,83 | 36,03 | 49,03 | 0,70 | 78,91 | 0,70 | 18,10 | 8,20 | 36,07 | 55,74 | 22,95 | 9,84 | 13,12 | 49,18 | 49 | 0 | 0 | 0 | 0 | Otro | 918 | 674 | 5702,10 | MuyBaja | -15,58 | 0,62 | 74 |
| Acehúche | 10 | Extremadura | 822 | 704 | 30,11 | 1 | 50,81 | 45,33 | 0,00 | 1 | 0 | 4,26 | 14,96 | 60,10 | 24,94 | 51,09 | 0,12 | 93,92 | 0,49 | 5,11 | 7,41 | 61,11 | 31,48 | 16,67 | 5,56 | 16,67 | 59,26 | 50 | 0 | 0 | 0 | 0 | Otro | 599 | 842 | 9106,46 | MuyBaja | -2,38 | 1,37 | 66 |
| Aceituna | 10 | Extremadura | 623 | 540 | 30,18 | 1 | 44,56 | 49,87 | 0,80 | 0 | 1 | 3,53 | 15,57 | 59,39 | 25,04 | 48,15 | 0,64 | 93,26 | 0,16 | 4,17 | 15,38 | 48,08 | 36,54 | 21,15 | 0,00 | 11,54 | 61,54 | 22 | 0 | 0 | 0 | 0 | Otro | 394 | 625 | 4007,61 | MuyBaja | -0,32 | 1,58 | 96 |
| Ahigal | 10 | Extremadura | 1421 | 1263 | 22,57 | 0 | 53,07 | 45,50 | 0,41 | 1 | 0 | 3,17 | 13,30 | 56,51 | 30,19 | 47,78 | 0,56 | 92,96 | 0,63 | 6,33 | 4,22 | 47,89 | 47,89 | 21,13 | 5,63 | 12,68 | 56,34 | 90 | 5 | 18 | 56 | 11 | ComercTTEHosteleria | 1023 | 1451 | 5207,35 | MuyBaja | -2,07 | 1,39 | 281 |
Ya en primer lugar llama la atención el inmenso número de columnas que tenermos (41), así que nos enfrentamos a un problema de alta dimensionalidad. Como por algo hay que empezar, voy a elegir mis variables objetivo. Las que he elegido para predecir son {Izda_Pct y Derecha} así ambas predicciones a parte de distinguirse por su tipología lo harán por su naturaleza y esto puede resultar en variaciones interesantes entre los modelos. Además voy a generar una variable identificativa combinando el código de provincia y el nombre, que usaré para denominar las filas, de esta forma me puedo deshacer de la variable Name.
ID<- c()
for(i in c(0:nrow(datos))){
ID<- c(ID, paste(datos$CodigoProvincia[i], datos$Name[i], sep='-'))
}
rownames(datos)<- ID
valid.cols = c(c(2:5), c(8), c(12:ncol(datos)))
datos<- datos[,valid.cols]
summary(datos)## CodigoProvincia CCAA Population TotalCensus
## Min. : 1.00 Length:8119 Min. : 5 Min. : 5
## 1st Qu.:13.00 Class :character 1st Qu.: 166 1st Qu.: 140
## Median :26.00 Mode :character Median : 549 Median : 447
## Mean :26.67 Mean : 5742 Mean : 4261
## 3rd Qu.:41.00 3rd Qu.: 2428 3rd Qu.: 1846
## Max. :52.00 Max. :3141991 Max. :2363829
##
## Izda_Pct Derecha Age_0-4_Ptge Age_under19_Ptge
## Min. : 0.00 Min. :0.0000 Min. : 0.000 Min. : 0.000
## 1st Qu.:21.89 1st Qu.:0.0000 1st Qu.: 1.389 1st Qu.: 8.334
## Median :35.16 Median :1.0000 Median : 2.978 Median :13.889
## Mean :34.40 Mean :0.6211 Mean : 3.019 Mean :13.568
## 3rd Qu.:46.03 3rd Qu.:1.0000 3rd Qu.: 4.533 3rd Qu.:19.058
## Max. :94.12 Max. :1.0000 Max. :13.245 Max. :33.696
##
## Age_19_65_pct Age_over65_pct WomanPopulationPtge ForeignersPtge
## Min. : 23.46 Min. : 0.00 Min. :11.77 Min. :-8.96
## 1st Qu.: 53.84 1st Qu.:19.82 1st Qu.:45.73 1st Qu.: 1.06
## Median : 58.66 Median :27.56 Median :48.48 Median : 3.59
## Mean : 57.37 Mean :29.07 Mean :47.30 Mean : 5.62
## 3rd Qu.: 61.82 3rd Qu.:36.91 3rd Qu.:50.00 3rd Qu.: 8.18
## Max. :100.00 Max. :76.47 Max. :72.68 Max. :71.47
##
## SameComAutonPtge SameComAutonDiffProvPtge DifComAutonPtge
## Min. : 0.00 Min. : 0.000 Min. : 0.000
## 1st Qu.: 75.81 1st Qu.: 0.676 1st Qu.: 4.933
## Median : 84.49 Median : 2.190 Median : 8.271
## Mean : 81.63 Mean : 4.337 Mean : 10.729
## 3rd Qu.: 90.46 3rd Qu.: 5.277 3rd Qu.: 13.898
## Max. :127.16 Max. :67.308 Max. :100.000
##
## UnemployLess25_Ptge Unemploy25_40_Ptge UnemployMore40_Ptge
## Min. : 0.000 Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.000 1st Qu.: 28.57 1st Qu.: 41.67
## Median : 5.882 Median : 39.94 Median : 50.00
## Mean : 7.322 Mean : 37.00 Mean : 50.18
## 3rd Qu.: 10.470 3rd Qu.: 46.67 3rd Qu.: 60.04
## Max. :100.000 Max. :100.00 Max. :100.00
##
## AgricultureUnemploymentPtge IndustryUnemploymentPtge
## Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 3.493 Median : 7.143
## Mean : 8.401 Mean : 10.008
## 3rd Qu.: 11.732 3rd Qu.: 14.286
## Max. :100.000 Max. :100.000
##
## ConstructionUnemploymentPtge ServicesUnemploymentPtge totalEmpresas
## Min. : 0.000 Min. : 0.00 Min. : 0.0
## 1st Qu.: 0.000 1st Qu.: 50.00 1st Qu.: 7.0
## Median : 8.333 Median : 62.02 Median : 30.0
## Mean : 10.838 Mean : 58.65 Mean : 398.6
## 3rd Qu.: 14.286 3rd Qu.: 72.12 3rd Qu.: 147.0
## Max. :100.000 Max. :100.00 Max. :299397.0
## NA's :5
## Industria Construccion ComercTTEHosteleria
## Min. : 0.00 Min. : 0.00 Min. : 0.0
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.0
## Median : 0.00 Median : 0.00 Median : 0.0
## Mean : 23.42 Mean : 48.88 Mean : 146.7
## 3rd Qu.: 14.00 3rd Qu.: 25.00 3rd Qu.: 65.0
## Max. :10521.00 Max. :30343.00 Max. :80856.0
## NA's :188 NA's :139 NA's :9
## Servicios ActividadPpal inmuebles Pob2010
## Min. : 0.0 Length:8119 Min. : 6 Min. : 5
## 1st Qu.: 0.0 Class :character 1st Qu.: 180 1st Qu.: 178
## Median : 0.0 Mode :character Median : 486 Median : 582
## Mean : 172.2 Mean : 3246 Mean : 5796
## 3rd Qu.: 40.0 3rd Qu.: 1589 3rd Qu.: 2483
## Max. :177677.0 Max. :1615548 Max. :3273049
## NA's :62 NA's :138 NA's :7
## SUPERFICIE Densidad PobChange_pct
## Min. : 2.58 Length:8119 Min. :-52.2700
## 1st Qu.: 1839.19 Class :character 1st Qu.:-10.4000
## Median : 3487.74 Mode :character Median : -4.9600
## Mean : 6214.70 Mean : -4.8974
## 3rd Qu.: 6893.88 3rd Qu.: 0.0925
## Max. :175022.91 Max. :138.4600
## NA's :9 NA's :7
## PersonasInmueble Explotaciones
## Min. :0.110 Min. : 1
## 1st Qu.:0.850 1st Qu.: 22
## Median :1.250 Median : 52
## Mean :1.296 Mean : 2447
## 3rd Qu.:1.730 3rd Qu.: 137
## Max. :3.330 Max. :99999
## NA's :138En el ‘summary’ detectamos varias cosas que tendremos que tener en cuenta:
Las columnas totalEmpresas, Industria, Construccion, ComercTTEHosteleria, Servicios, inmuebles, Pob2010, SUPERFICIE, PobChange_pct y PersonasInmueble presentan datos ‘missing’ en proporción variable que tendré que tratar. Además las columnas Densidad, ForeignersPtge toman valores sospechosos como ‘?’ y negativos por lo que los pasaré a datos ‘missing’.
Las columnas CCAA, ActividadPpal y Densidad son tipo ‘character’ y seguramente deba pasarlas a tipo ‘factor’.
Varias columnas numéricas tienen un valor máximo muy alejado del 3º cuartil, resultan especialmente sospechosas Industria, Construcción, ComercTTEHostelería y Servicios, donde más de la mitad de los datos son ceros y el máximo es del orden de las decenas de millar. Por su parte en Explotaciones el máximo es 99999, y parece un input incorrecto.
Las variables tipo porcentaje acaban en Ptge y en pct lo que resulta confuso. Además las referidas al desempleo alcanzan un valor máximo del 100% lo que me hace pensar que se trate de medidas relativas.
1.2 Depuración de los datos
1.2.1 Definición de categóricas
Lo primero que hacemos es pasar a factores las columnas tipo ‘character’ y comprobar que no hay columnas numéricas con una cantidad de valores suficientemente baja como para que nos interese convertirlas también.
datos[,c(2,6,28,32)] <- lapply(datos[,c(2,6,28,32)], factor)
# Cuento el n?mero de valores diferentes para las numéricas
sapply(Filter(is.numeric, datos),function(x) length(unique(x)))## CodigoProvincia Population
## 52 3597
## TotalCensus Izda_Pct
## 3310 6569
## Age_0-4_Ptge Age_under19_Ptge
## 3761 5891
## Age_19_65_pct Age_over65_pct
## 6215 6778
## WomanPopulationPtge ForeignersPtge
## 4524 2329
## SameComAutonPtge SameComAutonDiffProvPtge
## 6151 4207
## DifComAutonPtge UnemployLess25_Ptge
## 5574 2342
## Unemploy25_40_Ptge UnemployMore40_Ptge
## 2681 2751
## AgricultureUnemploymentPtge IndustryUnemploymentPtge
## 2525 2538
## ConstructionUnemploymentPtge ServicesUnemploymentPtge
## 2505 2904
## totalEmpresas Industria
## 1226 308
## Construccion ComercTTEHosteleria
## 457 803
## Servicios inmuebles
## 758 3088
## Pob2010 SUPERFICIE
## 3625 8110
## PobChange_pct PersonasInmueble
## 3049 283
## Explotaciones
## 758| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CodigoProvincia | 1 | 8119 | 26,67 | 14,90 | 26,00 | 26,71 | 20,76 | 1,00 | 52,00 | 51,00 | 0,01 | -1,32 | 0,17 |
| Population | 2 | 8119 | 5741,85 | 46215,20 | 549,00 | 1401,04 | 698,30 | 5,00 | 3141991,00 | 3141986,00 | 45,98 | 2814,43 | 512,90 |
| TotalCensus | 3 | 8119 | 4260,67 | 34428,89 | 447,00 | 1079,86 | 560,42 | 5,00 | 2363829,00 | 2363824,00 | 46,49 | 2888,34 | 382,10 |
| Izda_Pct | 4 | 8119 | 34,40 | 16,48 | 35,16 | 34,35 | 17,91 | 0,00 | 94,12 | 94,12 | 0,06 | -0,49 | 0,18 |
| Age_0-4_Ptge | 5 | 8119 | 3,02 | 2,05 | 2,98 | 2,94 | 2,33 | 0,00 | 13,24 | 13,24 | 0,34 | -0,21 | 0,02 |
| Age_under19_Ptge | 6 | 8119 | 13,57 | 6,78 | 13,89 | 13,68 | 7,90 | 0,00 | 33,70 | 33,70 | -0,10 | -0,79 | 0,08 |
| Age_19_65_pct | 7 | 8119 | 57,37 | 6,82 | 58,66 | 57,92 | 5,56 | 23,46 | 100,00 | 76,54 | -0,81 | 2,15 | 0,08 |
| Age_over65_pct | 8 | 8119 | 29,07 | 11,75 | 27,56 | 28,32 | 12,33 | 0,00 | 76,47 | 76,47 | 0,60 | 0,07 | 0,13 |
| WomanPopulationPtge | 9 | 8119 | 47,30 | 4,36 | 48,48 | 47,89 | 2,76 | 11,77 | 72,68 | 60,92 | -1,67 | 5,80 | 0,05 |
| ForeignersPtge | 10 | 8119 | 5,62 | 7,35 | 3,59 | 4,54 | 4,60 | -8,96 | 71,47 | 80,43 | 2,50 | 11,34 | 0,08 |
| SameComAutonPtge | 11 | 8119 | 81,63 | 12,29 | 84,49 | 83,22 | 10,20 | 0,00 | 127,16 | 127,16 | -1,52 | 3,47 | 0,14 |
| SameComAutonDiffProvPtge | 12 | 8119 | 4,34 | 6,39 | 2,19 | 2,97 | 2,95 | 0,00 | 67,31 | 67,31 | 3,29 | 14,55 | 0,07 |
| DifComAutonPtge | 13 | 8119 | 10,73 | 8,85 | 8,27 | 9,34 | 5,94 | 0,00 | 100,00 | 100,00 | 2,42 | 9,65 | 0,10 |
| UnemployLess25_Ptge | 14 | 8119 | 7,32 | 9,41 | 5,88 | 5,81 | 8,72 | 0,00 | 100,00 | 100,00 | 4,15 | 31,63 | 0,10 |
| Unemploy25_40_Ptge | 15 | 8119 | 37,00 | 20,32 | 39,94 | 37,04 | 12,39 | 0,00 | 100,00 | 100,00 | 0,21 | 1,41 | 0,23 |
| UnemployMore40_Ptge | 16 | 8119 | 50,18 | 22,80 | 50,00 | 50,93 | 14,12 | 0,00 | 100,00 | 100,00 | -0,23 | 0,85 | 0,25 |
| AgricultureUnemploymentPtge | 17 | 8119 | 8,40 | 12,96 | 3,49 | 5,69 | 5,18 | 0,00 | 100,00 | 100,00 | 3,23 | 15,56 | 0,14 |
| IndustryUnemploymentPtge | 18 | 8119 | 10,01 | 12,53 | 7,14 | 7,85 | 10,59 | 0,00 | 100,00 | 100,00 | 3,09 | 16,04 | 0,14 |
| ConstructionUnemploymentPtge | 19 | 8119 | 10,84 | 13,28 | 8,33 | 8,44 | 10,08 | 0,00 | 100,00 | 100,00 | 3,09 | 14,61 | 0,15 |
| ServicesUnemploymentPtge | 20 | 8119 | 58,65 | 24,26 | 62,02 | 60,69 | 17,72 | 0,00 | 100,00 | 100,00 | -0,81 | 0,80 | 0,27 |
| totalEmpresas | 21 | 8114 | 398,60 | 4219,37 | 30,00 | 85,36 | 44,48 | 0,00 | 299397,00 | 299397,00 | 53,68 | 3472,00 | 46,84 |
| Industria | 22 | 7931 | 23,42 | 158,61 | 0,00 | 6,39 | 0,00 | 0,00 | 10521,00 | 10521,00 | 44,26 | 2642,01 | 1,78 |
| Construccion | 23 | 7980 | 48,88 | 421,86 | 0,00 | 11,57 | 0,00 | 0,00 | 30343,00 | 30343,00 | 52,56 | 3503,51 | 4,72 |
| ComercTTEHosteleria | 24 | 8110 | 146,74 | 1233,02 | 0,00 | 31,86 | 0,00 | 0,00 | 80856,00 | 80856,00 | 45,40 | 2646,94 | 13,69 |
| Servicios | 25 | 8057 | 172,15 | 2446,81 | 0,00 | 22,01 | 0,00 | 0,00 | 177677,00 | 177677,00 | 57,48 | 3830,74 | 27,26 |
| inmuebles | 26 | 7981 | 3246,16 | 24314,71 | 486,00 | 927,86 | 555,98 | 6,00 | 1615548,00 | 1615542,00 | 44,53 | 2643,65 | 272,17 |
| Pob2010 | 27 | 8112 | 5795,81 | 47535,68 | 582,00 | 1433,66 | 730,92 | 5,00 | 3273049,00 | 3273044,00 | 47,15 | 2939,56 | 527,78 |
| SUPERFICIE | 28 | 8110 | 6214,70 | 9218,19 | 3487,74 | 4414,91 | 3029,12 | 2,58 | 175022,91 | 175020,33 | 6,07 | 62,28 | 102,36 |
| PobChange_pct | 29 | 8112 | -4,90 | 10,38 | -4,96 | -5,11 | 7,83 | -52,27 | 138,46 | 190,73 | 1,50 | 15,08 | 0,12 |
| PersonasInmueble | 30 | 7981 | 1,30 | 0,57 | 1,25 | 1,28 | 0,64 | 0,11 | 3,33 | 3,22 | 0,26 | -0,63 | 0,01 |
| Explotaciones | 31 | 8119 | 2447,20 | 15062,74 | 52,00 | 80,85 | 54,86 | 1,00 | 99999,00 | 99998,00 | 6,32 | 37,95 | 167,17 |
De momento voy a dejarlo así, ya que los datos son funcionales, pero me incomodan la columna CCAA, sobre la que trataré de agrupar categorías con cierto criterio más adelante, y la de CódigoProvincia, esta columna implica una agrupación de los datos, pero contiene demasiadas categorías (52), y puede aportar información aunque haya sido empleada como identificativa, así que voy a conservarla como numérica, pese a que no me entusiasma la idea.
1.2.2 Tratamiento de ‘inputs’ erróneos
Lo siguiente que voy a hacer es convertir en ‘NAs’ esos valores que hemos visto que parecían inputs incorrectos, además voy a echarle un ojo a las variables cualitativas para ver si puedo agrupar categorías poco representadas.
# Missings no declarados (99999, ?)
datos$Explotaciones<-recode.na(datos$Explotaciones,99999,as.numeric = T)
datos$Densidad<-recode.na(datos$Densidad,"?")
# Valores fuera de rango
datos$ForeignersPtge<-replace(datos$ForeignersPtge, which(datos$ForeignersPtge < 0), NA)1.2.3 Agrupación de categorías
Miramos las variables categóricas para detectar categorías poco representadas. La forma de categorizar CCAA no es trivial ya que sabemos que subyacen relaciones geográficas, demográficas, históricas, … que debemos tener en cuenta. En este caso, dado que no podía extenderme en exceso, he decidido crear 4 categorías atendiendo esencialmente a variables demográficas. Primero agrupo el resto de variables categóricas. 
Vemos que lo lógico es agrupar Servicios, Construcción e Industria (y aprovechar para renombrar Hosteleria) para la variable ActividadPpal. Por su parte en densidad vamos a agrupar Alta y Baja como Alta y a renombrar MuyBaja como Baja. De esta forma reducimos las categorías en ambos casos, teniendo categorías más representadas pero que conservan el significado de las categorías originales.
# Agrupamos las categorías poco representadas
# Interesa agrupar Construcción/Industria y Servicios
datos$ActividadPpal<- car::recode(datos$ActividadPpal,
"c('Construccion', 'Industria', 'Servicios')='Const/Indus/Servi';
c('ComercTTEHosteleria') ='Hosteleria'")
datos$Densidad<- car::recode(datos$Densidad,
"c('Alta', 'Baja')='Alta';c('MuyBaja') ='Baja'")
# Renombro además las variables asociadas
names(datos)[26]<- 'Hosteleria'
Las categorías de ActividadPpal tenían asociada su propia columna, que a priori pensé que procedían de generar ‘dummies’ pero toman valores bastante extraños, en Construcción el 61% de los datos son 0, en Servicios el 61%, en Industria el 61% y en Hosteleria el 61%, son porcentajes muy altos para columnas que tienen valores máximos del orden de las decenas de millar, así que voy a generarme una variable nueva en la que diga si el valor es mayor que cero o no, que creo que es la información que realmente podemos extraer de estas columnas.
Abordo ahora la agrupación de la variable Comunidad Autónoma, que he realizado basándome en mi propio criterio sobre una preagrupación en las variables Population, Age_under19_Ptge y ForeignersPtge (lo que sería un ‘clustering’ a ojímetro).
datos$ActPpal.bin <- as.factor(apply(ifelse(datos[,c(24:27)]<=0,0,1), 1, prod, na.rm=T))
# Agrupo las comunidades autónomas
CCAA.grouped <- as.data.frame(datos[,c(2,3,8,12)] %>%
group_by(CCAA) %>%
summarise(
Population = mean(Population, na.rm=T),
Age_under19_Ptge = mean(Age_under19_Ptge, na.rm=T),
ForeignersPtge = mean(ForeignersPtge, na.rm=T)
))
rownames(CCAA.grouped) <- CCAA.grouped$CCAA
CCAA.grouped <- CCAA.grouped[,c(2:4)]
factor_3 <- function (values, desired.labels){
lim0 <- min(values)
lim1 <- mean(values) - sqrt(var(values))
lim3 <- mean(values) + sqrt(var(values))
lim4 <- max(values)
category <- cut(values, c(lim0,lim1,lim3,lim4), include.lowest=c(T,F,F), labels=desired.labels)
}
CCAA.grouped$Population <- factor_3(CCAA.grouped$Population, c("poco poblado", "normal poblado",
"muy poblado"))
CCAA.grouped$Age_under19_Ptge <- factor_3(CCAA.grouped$Age_under19_Ptge, c("poco joven", "normal joven",
"muy joven"))
CCAA.grouped$ForeignersPtge <- factor_3(CCAA.grouped$ForeignersPtge, c("pocos extranjero", "normal extranjero",
"muy extranjero"))
datos$CCAA.grouped<- car::recode(datos$CCAA,
"c('Asturias', 'Cantabria', 'Extremadura', 'Galicia')='PE-NJ';
c('Andalucía', 'Canarias', 'CastillaMancha','Cataluña', 'Ceuta',
'ComValenciana', 'Navarra', 'PaísVasco') ='NE-NJ';
c('Aragón', 'CastillaLeón', 'Rioja')='NE-PJ';
c('Madrid', 'Baleares', 'Melilla', 'Murcia')='ME-NJ'")| Population | Age_under19_Ptge | ForeignersPtge | |
|---|---|---|---|
| Asturias | normal poblado | normal joven | pocos extranjero |
| Cantabria | normal poblado | normal joven | pocos extranjero |
| Extremadura | normal poblado | normal joven | pocos extranjero |
| Galicia | normal poblado | normal joven | pocos extranjero |
| Andalucía | normal poblado | normal joven | normal extranjero |
| Aragón | normal poblado | poco joven | normal extranjero |
| Canarias | normal poblado | normal joven | normal extranjero |
| CastillaLeón | poco poblado | poco joven | normal extranjero |
| CastillaMancha | normal poblado | normal joven | normal extranjero |
| Cataluña | normal poblado | normal joven | normal extranjero |
| Ceuta | muy poblado | muy joven | normal extranjero |
| ComValenciana | normal poblado | normal joven | normal extranjero |
| Navarra | normal poblado | normal joven | normal extranjero |
| PaísVasco | normal poblado | normal joven | normal extranjero |
| Rioja | normal poblado | poco joven | normal extranjero |
| Baleares | normal poblado | normal joven | muy extranjero |
| Madrid | normal poblado | normal joven | muy extranjero |
| Melilla | muy poblado | muy joven | muy extranjero |
| Murcia | normal poblado | normal joven | muy extranjero |
Las agrupaciones las he decidido basándome en la última tabla, el porcentaje de extranjeros define las categorías y dentro de las comunidades ‘normal extranjero’ las que tienen ‘normal joven’ y ‘muy joven’ van juntas y las ‘poco joven’ aparte. Visto así nos quedan los siguientes 4 grupos: (Asturias, Cantabria, Extremadura y Galicia), (Andalucía, Canarias, Castilla la Mancha, Cataluña, Ceuta, Valencia, Navarra y País Vasco), (Aragón, Castilla y León y la Rioja) y (Baleares, Madrid, Melilla y Murcia). Vemos que los grupos son de tamaños razonables y de alguna forma es más preciso que una agrupación basada únicamente en nuestros prejuicios, por supuesto están generados a partir de otros datos, luego recogen la información ya presente en estos, en cualquier caso conservaré tanto la columna original como la agrupada y ya veremos cual aporta mayor información.
1.2.4 Tratamiento de ‘Outliers’
A partir de este momento ya separo mis variables objetivo de mis variables input y continuo trabajando únicamente sobre las variables input. En primer lugar voy a calcular el porcentaje de ‘outliers’ presentes en cada variable, si los porcentajes son pequeños, los sustituiré por valores ‘missing’.
varObjCont<- datos$Izda_Pct
varObjBin<- datos$Derecha
input<- as.data.frame(datos[,-c(5,6)])
# Porcentaje por variable
sapply(Filter(is.numeric, input), function(x) atipicosAmissing(x)[[2]])/nrow(input)## CodigoProvincia Population
## 0.0000000000 0.0992733095
## TotalCensus Age_0-4_Ptge
## 0.0961941126 0.0000000000
## Age_under19_Ptge Age_19_65_pct
## 0.0000000000 0.0029560291
## Age_over65_pct WomanPopulationPtge
## 0.0000000000 0.0025865254
## ForeignersPtge SameComAutonPtge
## 0.0052962187 0.0001231679
## SameComAutonDiffProvPtge DifComAutonPtge
## 0.0203226998 0.0049267151
## UnemployLess25_Ptge Unemploy25_40_Ptge
## 0.0032023648 0.0000000000
## UnemployMore40_Ptge AgricultureUnemploymentPtge
## 0.0000000000 0.0199531962
## IndustryUnemploymentPtge ConstructionUnemploymentPtge
## 0.0059120581 0.0065278975
## ServicesUnemploymentPtge totalEmpresas
## 0.0000000000 0.1050621998
## Industria Construccion
## 0.0874491933 0.0891735435
## Hosteleria Servicios
## 0.0985343023 0.1187338342
## inmuebles Pob2010
## 0.0886808720 0.0975489592
## SUPERFICIE PobChange_pct
## 0.0270969331 0.0012316788
## PersonasInmueble Explotaciones
## 0.0000000000 0.0454489469# Modifico los atípicos como missings
input[,as.vector(which(sapply(input, class)=="numeric"))]<- sapply(Filter(is.numeric, input),
function(x) atipicosAmissing(x)[[1]])Podemos ver que las variables con mayor porcentaje de ‘outliers’ apenas rondan el 10%. Seguramente incorporarlos como ‘missing’ favorezca el comportamiento del modelo. Por supuesto este proceso podría repetirse ad-infinitum ya que si busco ‘outliers’ en el nuevo ‘set’ aparecerán nuevos que antes no lo eran. Tal vez la aplicación de estas líneas sobre el ‘set’ de datos de forma recurrente hasta alcanzar una tolerancia podría resultar en unos datos exquisitamente centrados pero los estarían tan maquillados que habríamos perdido demasiada información. Una única aplicación para eliminar el factor lotería parece lo razonable.
1.2.5 Tratamiento de ‘missings’
Una vez tratadas todas las posibles fuentes de datos ‘missing’ tenemos que abordar estos casos. En primer lugar debemos mirar nuestra estructura de datos, para ello voy a representar la correlación entre los datos perdidos y no perdidos (en binario) de nuestras variables. 
Hay varias variables que aparecen fuertementemente correlacionadas, en principio la correlacion entre Population, TotalCensus y Pob2010 es trivial y debería conservarse al tratar las propias variables que seguramente estén muy correlacionadas, incluso podríamos considerarlas dependientes. Por otro lado es bastante sugerente que allí donde a estas variables les faltan datos les falten también a totalEmpresas, Inmuebles y Industria, Construccion, Hosteleria y Servicios. Deben existir algunas localidades que no recogen estos datos o que no están queriendo publicarlos, así que añadiremos una variable con la proporción de missings de cada localidad, ya que si podría resultar indicativa. De igual forma miraremos el porcentaje de ‘missings’ en cada columna, por si alguna conviniese directamente eliminarla.
input$prop_missings<-apply(is.na(input),1,mean)
summary(input$prop_missings)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.03508 0.02857 0.34286apply(is.na(input), 2, mean)## CodigoProvincia CCAA
## 0.0000000000 0.0000000000
## Population TotalCensus
## 0.0992733095 0.0961941126
## Age_0-4_Ptge Age_under19_Ptge
## 0.0000000000 0.0000000000
## Age_19_65_pct Age_over65_pct
## 0.0029560291 0.0000000000
## WomanPopulationPtge ForeignersPtge
## 0.0025865254 0.0857248430
## SameComAutonPtge SameComAutonDiffProvPtge
## 0.0001231679 0.0203226998
## DifComAutonPtge UnemployLess25_Ptge
## 0.0049267151 0.0032023648
## Unemploy25_40_Ptge UnemployMore40_Ptge
## 0.0000000000 0.0000000000
## AgricultureUnemploymentPtge IndustryUnemploymentPtge
## 0.0199531962 0.0059120581
## ConstructionUnemploymentPtge ServicesUnemploymentPtge
## 0.0065278975 0.0000000000
## totalEmpresas Industria
## 0.1056780392 0.1106047543
## Construccion Hosteleria
## 0.1062938786 0.0996428132
## Servicios ActividadPpal
## 0.1263702426 0.0000000000
## inmuebles Pob2010
## 0.1056780392 0.0984111344
## SUPERFICIE Densidad
## 0.0282054440 0.0113314448
## PobChange_pct PersonasInmueble
## 0.0020938539 0.0169971671
## Explotaciones ActPpal.bin
## 0.0687276758 0.0000000000
## CCAA.grouped prop_missings
## 0.0000000000 0.0000000000Entre las columnas numéricas los valores mayores se encuentran en torno al 10% así que en principio voy a conservarlas todas, eso sí, inputaré valores en esos missings. Para las categóricas solo hay ‘missings’ en Densidad, algo más de un 1%, así que también inputaré los valores. De entre las filas hasta el 75% de las observaciones está por debajo del 4% de características ‘missing’, sin embargo preocupa el máximo, que es del 35%. De momento voy a dejarlas todas ya que ninguna supera el límite del 50%, igual que en las columnas sustituiré estos ‘missings’ por valores similares a los observados en la varible. Dependiendo de como estén distribuidos los datos en las variables inputaré aleatorios, valores iguales a la media o a la mediana, así que primero miro las distribuciones de las variables numéricas que voy a imputar.
| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Population | 1 | 7313 | 1255,87 | 1840,51 | 431,00 | 817,60 | 521,88 | 5,00 | 9208,00 | 9203,00 | 2,18 | 4,42 | 21,52 |
| TotalCensus | 2 | 7338 | 988,49 | 1415,14 | 358,00 | 653,92 | 424,02 | 5,00 | 6958,00 | 6953,00 | 2,15 | 4,25 | 16,52 |
| Age_19_65_pct | 3 | 8095 | 57,44 | 6,60 | 58,67 | 57,94 | 5,55 | 30,36 | 85,72 | 55,36 | -0,75 | 1,13 | 0,07 |
| WomanPopulationPtge | 4 | 8098 | 47,36 | 4,17 | 48,49 | 47,91 | 2,75 | 26,47 | 69,23 | 42,76 | -1,39 | 3,53 | 0,05 |
| ForeignersPtge | 5 | 7423 | 6,09 | 6,32 | 4,10 | 5,03 | 4,48 | 0,00 | 40,09 | 40,09 | 1,86 | 4,46 | 0,07 |
| SameComAutonPtge | 6 | 8118 | 81,64 | 12,26 | 84,49 | 83,22 | 10,20 | 15,41 | 127,16 | 111,74 | -1,50 | 3,30 | 0,14 |
| SameComAutonDiffProvPtge | 7 | 7954 | 3,70 | 4,52 | 2,12 | 2,79 | 2,83 | 0,00 | 25,73 | 25,73 | 2,02 | 4,48 | 0,05 |
| DifComAutonPtge | 8 | 8079 | 10,46 | 7,98 | 8,24 | 9,27 | 5,89 | 0,00 | 55,56 | 55,56 | 1,76 | 4,24 | 0,09 |
| UnemployLess25_Ptge | 9 | 8093 | 7,03 | 7,84 | 5,86 | 5,77 | 8,69 | 0,00 | 75,00 | 75,00 | 2,10 | 8,10 | 0,09 |
| AgricultureUnemploymentPtge | 10 | 7957 | 7,19 | 9,37 | 3,23 | 5,31 | 4,78 | 0,00 | 46,43 | 46,43 | 1,62 | 2,21 | 0,11 |
| IndustryUnemploymentPtge | 11 | 8071 | 9,47 | 10,46 | 7,01 | 7,74 | 10,40 | 0,00 | 82,00 | 82,00 | 1,70 | 4,15 | 0,12 |
| ConstructionUnemploymentPtge | 12 | 8066 | 10,25 | 11,18 | 8,33 | 8,32 | 9,88 | 0,00 | 86,49 | 86,49 | 2,04 | 6,14 | 0,12 |
| totalEmpresas | 13 | 7261 | 72,93 | 113,95 | 23,00 | 45,29 | 34,10 | 0,00 | 566,00 | 566,00 | 2,28 | 4,88 | 1,34 |
| Industria | 14 | 7221 | 6,23 | 11,84 | 0,00 | 3,22 | 0,00 | 0,00 | 56,00 | 56,00 | 2,17 | 4,18 | 0,14 |
| Construccion | 15 | 7256 | 11,15 | 21,17 | 0,00 | 5,77 | 0,00 | 0,00 | 100,00 | 100,00 | 2,17 | 4,16 | 0,25 |
| Hosteleria | 16 | 7310 | 28,67 | 54,07 | 0,00 | 14,95 | 0,00 | 0,00 | 260,00 | 260,00 | 2,19 | 4,34 | 0,63 |
| Servicios | 17 | 7093 | 15,90 | 32,48 | 0,00 | 7,29 | 0,00 | 0,00 | 160,00 | 160,00 | 2,44 | 5,58 | 0,39 |
| inmuebles | 18 | 7261 | 882,30 | 1138,88 | 399,00 | 625,43 | 431,44 | 6,00 | 5811,00 | 5805,00 | 2,12 | 4,30 | 13,37 |
| Pob2010 | 19 | 7320 | 1300,47 | 1875,47 | 460,50 | 857,38 | 549,30 | 5,00 | 9361,00 | 9356,00 | 2,15 | 4,26 | 21,92 |
| SUPERFICIE | 20 | 7890 | 5080,93 | 4910,05 | 3375,00 | 4144,00 | 2859,54 | 2,58 | 27681,20 | 27678,62 | 1,92 | 3,87 | 55,28 |
| Densidad* | 21 | 8027 | 1,80 | 0,40 | 2,00 | 1,87 | 0,00 | 1,00 | 2,00 | 1,00 | -1,50 | 0,24 | 0,00 |
| PobChange_pct | 22 | 8102 | -5,01 | 9,80 | -4,97 | -5,13 | 7,81 | -52,27 | 54,05 | 106,32 | 0,45 | 3,31 | 0,11 |
| PersonasInmueble | 23 | 7981 | 1,30 | 0,57 | 1,25 | 1,28 | 0,64 | 0,11 | 3,33 | 3,22 | 0,26 | -0,63 | 0,01 |
| Explotaciones | 24 | 7561 | 83,48 | 93,74 | 46,00 | 64,17 | 45,96 | 1,00 | 465,00 | 464,00 | 1,86 | 3,11 | 1,08 |
Pese a que todas las distribuciones son leptocúrticas, y por lo tanto están de alguna forma concentradas en torno a la media, no voy a inputar siempre por medidas centrales. Tanto si hay mucha diferencia entre la media y la mediana como si el porcentaje de ‘missings’ es alto inputaré con datos aleatorios, y si la diferencia es pequeña (aqui ‘mucha’ y ‘pequeña’ obedecerán a mi criterio) inputaré mediante la mayor o menor de ellas dependiendo del signo del coeficiente de simetría. Por otro lado, variables cualitativas con ‘missings’ solo hay una, y con un pocentaje de ‘missings’ realmente bajo, así que haré la inputación aleatoria.
Todo aparece dentro de un bucle while porque no se me estaban sustituyendo todos los ‘missing’ en una sola ejecución
## Imputaciones
while(length(which(apply(is.na(input),2, mean) >0))){
# Imputo todas las cuantitativas con la media
inputando <- c("Age_19_65_pct", "SameComAutonDiffProvPtge", "DifComAutonPtge", "SUPERFICIE",
"PobChange_pct", 'WomanPopulationPtge', 'UnemployLess25_Ptge',
'IndustryUnemploymentPtge', 'ConstructionUnemploymentPtge', 'SameComAutonPtge',
'PersonasInmueble')
input[,inputando]<- sapply(input[, inputando], function(x) ImputacionCuant(x,"media"))
# Imputo todas las cuantitativas con aleatorios
inputando <- c("Population", "TotalCensus", "ForeignersPtge", "totalEmpresas", "Industria", "Construccion",
"Hosteleria", "Servicios", "inmuebles", "Pob2010", 'AgricultureUnemploymentPtge',
"Explotaciones")
input[,inputando]<- sapply(input[, inputando], function(x) ImputacionCuant(x,"aleatorio"))
# Imputo todas las cualitativas
input[,as.vector(which(sapply(input, class)=="factor"))]<-sapply(Filter(is.factor, input),
function(x) ImputacionCuali(x,"aleatorio"))
# A veces se cambia el tipo de factor a character al imputar, así que hay que indicarle que es factor
input[,as.vector(which(sapply(input, class)=="character"))] <- lapply(input[,as.vector(which(sapply(input, class)=="character"))] , factor)
}
# Reviso que no queden datos missings
which(apply(is.na(input),2,sum) >0)## named integer(0)Hecho esto doy por terminada la depuración, ya he corregido los errores que traían los datos y dispongo de unos datos tratables estadísticamente.
2 Exploración de los datos
Aunque durante la eliminación de los inputs erróneos y los outliers y la sustitución de los valores missings hemos alcanzado cierta perspectiva sobre el contenido de los datos, vamos ahora a explorarlos y transformalos para intentar facilitar la tarea de predicción al modelo.
2.1 Influencia de las variables
Lo primero que hago es añadir un par de variables aleatorias para establecer que me permita valorar el contenido de información de las demás. Por otro lado veo mediante la V de Cramer la asociación presente entre mis variable objetivo y el resto.
# Primera visualización ----
# Creo la variable aleatoria
input$aleatorio<-runif(nrow(input))
input$aleatorio2<-runif(nrow(input))
# Obtengo la importancia de las variables.
par(mfrow=c(1,2), mar=c(4.1, 14, 4.1, 2.1))
graficoVcramer(input,varObjBin, "V de Cramer para Derecha")
graficoVcramer(input,varObjCont, "V de Cramer para Izda_Pct")
par(mfrow=c(1,1)) Ya vemos que como nos barruntábamos la variable CCAA tiene una importancia muy destacada sobre el resto, nuestra agrupación no es del todo absurda ya que se mantiene en segunda lugar pero presenta una pérdida de información considerable respecto a la variable original. También debe preocuparnos la importancia del código de provincia, que a priori debería ser identificativa (nosotros la hemos usado como tal) pero por lo visto en nuestro país la influencia del entorno es muy fuerte (al menos en política). A la vista de que la variable ForeignersPtge aparece por debajo de las variables aleatorias para Derecha, y que esta variable la usé yo para generar mi agrupación de Comunidades Autónomas, voy a repetir la agrupación, mirando con un ojo la tabla y con el otro cerrado, y voy a ver genero una agrupación que sea más útil.
# Nueva agrupación
input$CCAA.grouped2<- car::recode(input$CCAA,
"c('Asturias', 'Cantabria', 'Rioja', 'Aragón', 'Galicia') = 'norteños';
c('Andalucía', 'Canarias', 'ComValenciana', 'Baleares') = 'sureños';
c('CastillaMancha','CastillaLeón', 'Extremadura', 'Murcia', 'Ceuta',
'Melilla') = 'intereños';
c('Madrid', 'Cataluña', 'PaísVasco', 'Navarra') ='pijeños'")
# Obtengo la importancia de las variables.
par(mfrow=c(1,2), mar=c(4.1, 14, 4.1, 2.1))
graficoVcramer(input,varObjBin, "V de Cramer para Derecha")
graficoVcramer(input,varObjCont, "V de Cramer para Izda_Pct")
par(mfrow=c(1,1))
# La nueva agrupación mejora a la anterior así que elimino la anterior
input$CCAA.grouped<- input$CCAA.grouped2
input$CCAA.grouped2<- NULLDado que la nueva agrupación ha mejorado la antigua será la que conserve. Posiblemente para generar ‘dummies’ pueda usar las Comunidades Autónomas como categorías pero su versión agrupada me será útil para introducir alguna interacción, ya que si quisiera hacer la interacción de CCAA y ActividadPpal (que parece que va a dar juego) tendría 57 combinaciones que supera el número de variables actual.
Mirando también el resto de variables resulta interesante ver que ActividadPpal y ActividadPpal.bin (que generé a partir de los ceros de las variables Industria, Hosteleria, etc.) son prácticamente igual de influyentes, incluso para la variable binaria hay mayor asociación con el factor generado por nosostros, mientras que Hosteleria, Industria, Servicios, … son variables que aparecen a la cola del gráfico.
2.1.1 Sobre la binaria
Visto esto procede ahora ver la distribución de las variables cualitativas y las cuantitativas más significativas sobre nuestras variables objetivo, empiezo por la binaria
# Nombro los factores
factores<- c("CCAA", "ActividadPpal", "Densidad", "ActPpal.bin", "CCAA.grouped")
# Veo el efecto en forma de mosaico (nop lo muestro por no excederme)
# for(name in factores){
# mosaico_targetbinaria(input[[name]],varObjBin,name)
# }
# Y ahora en un diagrama de barras
p1<- barras_targetbinaria(input[["ActividadPpal"]],varObjBin,"ActividadPpal")
p2<- barras_targetbinaria(input[["Densidad"]],varObjBin,"Densidad")
p3<- barras_targetbinaria(input[["ActPpal.bin"]],varObjBin,"ActPpal.bin")
p4<- barras_targetbinaria(input[["CCAA.grouped"]],varObjBin,"CCAA.grouped")
gridExtra::grid.arrange(p1, p2, p3, p4, ncol=4, nrow=1)
Debo indicar que no he representado la variable CCAA porque el elevado número de valores impedía sacar alguna conclusión. En este gráfico podemos ver que la agrupación que hicimos para ActividadPpal pudo haber sido más drástica, juntando Hostelería con Const/Indud/Servi y al menos las predicciones sobre Derecha dudo que cambiasen significativamente. En el resto de casos la distribución de ceros y unos entre las categorías si parece lo suficientemente dispar, salvo entre las categorías intereños y norteños de las recién agrupadas Comunidades Autónomas, pero agrupándolas aglutinaríamos en una categoría demasiados datos. Como el efecto de tener la variable ActividadPpal como binaria parece que ya lo tenemos recogido en ActPpal.bin creo que no introduciré cambios en los datos.
Toca ahora ver la influencia de las variables cuantitativas, para ello vamos a ver si hay diferencia entre como están distribuídos los ceros y los unos de nuestra variable objetivo en las ocho (hay 33, debía elegir una muestra) variables cuantitativas que salieron más relevantes en los gráficos de la V de Cramer.
#Veo gráficamente el efecto de ocho variables cuantitativas sobre la binaria
# dput(names(Filter(is.numeric, input)))
numericas.mainbin<- c("Age_under19_Ptge", "Age_over65_pct","PersonasInmueble",
"Age_0-4_Ptge", "CodigoProvincia", "Age_19_65_pct",
"PobChange_pct", "UnemployMore40_Ptge")
p1<- boxplot_targetbinaria(input$Age_under19_Ptge,varObjBin,"Age_under19_Ptge")
p2<- boxplot_targetbinaria(input$Age_over65_pct,varObjBin,"Age_over65_pct")
p3<- boxplot_targetbinaria(input$PersonasInmueble,varObjBin,"PersonasInmueble")
p4<- boxplot_targetbinaria(input[["Age_0-4_Ptge"]],varObjBin,"Age_0-4_Ptge")
p5<- boxplot_targetbinaria(input$CodigoProvincia,varObjBin,"CodigoProvincia")
p6<- boxplot_targetbinaria(input$Age_19_65_pct,varObjBin,"Age_19_65_pct")
p7<- boxplot_targetbinaria(input$PobChange_pct,varObjBin,"PobChange_pct")
p8<- boxplot_targetbinaria(input$UnemployMore40_Ptge,varObjBin,"UnemployMore40_Ptge")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, ncol=4, nrow=2)
# Histograma para las más interesantes
p1<- hist_targetbinaria(input$Age_under19_Ptge,varObjBin,"Age_under19_Ptge")
p2<- hist_targetbinaria(input$Age_over65_pct,varObjBin,"Age_over65_pct")
p3<- hist_targetbinaria(input$PersonasInmueble,varObjBin,"PersonasInmueble")
p4<- hist_targetbinaria(input[["Age_0-4_Ptge"]],varObjBin,"Age_0-4_Ptge")
p5<- hist_targetbinaria(input$CodigoProvincia,varObjBin,"CodigoProvincia")
p6<- hist_targetbinaria(input$Age_19_65_pct,varObjBin,"Age_19_65_pct")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, ncol=3, nrow=2)
Podemos ver que hasta Age_19_65_pct las distribuciones son bastante diferentes y es a partir de esta variable cuando las cajas comienzan a solaparse en el mismo rango de valores. De esta forma corroboramos que al menos de esas 5 primeras variables cuantitativas podremos extraer bastante información. Con los gráficos de densidad podemos ver de forma más precisa como las distribuciones de esas 5 o 6 variables cuantitativas que encabezaban la V de Cramer tienen amplias regiones de no solapamiento. No obstante tanto el CodigoProvincia como Age_19_65_pct muestran menor poder discriminante que el resto. En el caso CódigoProvincia se deberá seguro a que es una variable en la que los números tienen un valor cualitativo y no cuantitativo, no obstante adopta length(unique(input$CodigoProvincia)) valores diferentes lo que imposibilita un tratamiento diferente.
2.1.2 Sobre la contínua
Miro ahora el impacto de las 8 variables cuantitativas más influyentes según la V de Cramer sobre la variable objetivo contínua (Izda_Pct).
numericas.maincon<- c("UnemployLess25_Ptge", "SameComAutonDiffProvPtge","UnemployMore40_Ptge",
"AgricultureUnemploymentPtge", "Age_under19_Ptge", "Age_over65_pct",
"CodigoProvincia", "Age_0-4_Ptge")
#Principales variables numéricas frente a la objetivo continua
graficoCorrelacion(varObjCont,input[numericas.maincon])
par(mfrow=c(1,2), mar=c(4.1, 4.1, 4.1, 4.1))
corrplot(cor(cbind(varObjCont,Filter(is.numeric, input)), use="pairwise", method="pearson"), method = "ellipse",type = "upper")
corrplot(cor(cbind(varObjCont,input[numericas.maincon]), use="pairwise", method="pearson"), method = "ellipse",type = "upper")
par(mfrow=c(1,1))De estos tres gráficos podemos extraer bastante información. El superior nos muestra como se distribuyen las variables unas en función de las otras. Tanto mirando la primera linea (donde están directamente representados los pares de puntos y una línea de tendencia) como la primera columna (donde aparece el índice de correlación lineal) podemos ver como nuestra variable objetivo no presenta una correlación fuerte con ninguna variable cuantitativa. Si se aprecia ya una tendencia fuerte entre otras variables, pero nada sorprendente ya que todas nuestras variables son medidas de edad y empleo, que son supercategorías que están además relacionadas entre sí. Estas relaciones entre variables podría conducirnos a eliminarlas, sin embargo, para ello emplearemos métodos destinados a ese fin. Esta etapa es fundamentalmente descriptiva y clave para garantizar la consistencia en la construcción del modelo.
En la pareja de gráficas inferiores vemos primero la correlación entre todas las variables. Aquí se ven los grupos relacionados con la edad y el empleo fuertemente relacionados no sólo internamente sino también entre sí. En segundo lugar vemos como según vamos avanzando en la primera fila el color de las variables va perdiendo fuerza, corroborando este gráfico la información extraida de la V de Cramer (las variables están ordenadas).
3 Modelos predictivos
Toca por fin elavorar los modelos predictivos sobre nuestros datos, ya perfectamente preparados. La idea es elaborar múltiples modelos que vayan ganando complejidad e ir viendo que variables son las que tienen mayor impacto, para acabar teniendo un modelo lo más competitivo posible.
3.1 Regresión Lineal
En primer lugar vamos a enfocarnos en la variable contínua, lo primero que hacemos es cargar los datos (sin las transformaciones) y elaborar un primer modelo con todas las variables que nos sirva como tanteo.
3.1.1 Modelo Base
Elaboro los primeros modelos con todas las variables. Para evitar la multicolinealidad pruebo con las generadas por mí en el proceso de depuración y exploración y con las originales.
#Obtengo la partición
set.seed(123456)
names(input)[names(input) == "Age_0-4_Ptge"] <- "Age_0_4_Ptge"
todo<-data.frame(input,varObjCont)
trainIndex <- createDataPartition(todo$varObjCont, p=0.8, list=FALSE)
data_train <- todo[trainIndex,]
data_test <- todo[-trainIndex,]
#Construyo un modelo preliminar con todas las variables
data_train <- todo[trainIndex,-c(34,35)]
data_test <- todo[-trainIndex,-c(34,35)]
modelo1<-lm(varObjCont~.,data=data_train)
summary(modelo1)##
## Call:
## lm(formula = varObjCont ~ ., data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.854 -5.823 -0.339 5.533 52.757
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.897e+01 7.964e+00 2.382 0.017259 *
## CodigoProvincia 2.140e-02 9.702e-03 2.206 0.027426 *
## CCAAAragón -1.114e+01 6.965e-01 -15.999 < 2e-16 ***
## CCAAAsturias -5.164e+00 1.364e+00 -3.785 0.000155 ***
## CCAABaleares -1.027e+01 1.459e+00 -7.037 2.17e-12 ***
## CCAACanarias -1.517e+01 1.327e+00 -11.431 < 2e-16 ***
## CCAACantabria -1.790e+01 1.247e+00 -14.351 < 2e-16 ***
## CCAACastillaLeón -1.966e+01 5.747e-01 -34.212 < 2e-16 ***
## CCAACastillaMancha -1.199e+01 6.394e-01 -18.754 < 2e-16 ***
## CCAACataluña -4.330e+01 6.428e-01 -67.371 < 2e-16 ***
## CCAACeuta -2.209e+01 1.003e+01 -2.202 0.027691 *
## CCAAComValenciana -3.154e+01 6.908e-01 -45.663 < 2e-16 ***
## CCAAExtremadura -3.666e+00 7.325e-01 -5.004 5.76e-07 ***
## CCAAGalicia -3.118e+01 8.062e-01 -38.670 < 2e-16 ***
## CCAAMadrid -1.651e+01 1.024e+00 -16.116 < 2e-16 ***
## CCAAMelilla -1.842e+01 1.003e+01 -1.836 0.066407 .
## CCAAMurcia -1.773e+01 1.779e+00 -9.970 < 2e-16 ***
## CCAANavarra -1.348e+01 8.745e-01 -15.411 < 2e-16 ***
## CCAAPaísVasco -2.267e+01 8.757e-01 -25.892 < 2e-16 ***
## CCAARioja -1.602e+01 1.068e+00 -15.000 < 2e-16 ***
## Population -9.015e-05 1.558e-04 -0.579 0.562759
## TotalCensus 4.730e-04 1.945e-04 2.432 0.015040 *
## Age_0_4_Ptge -3.556e-01 1.241e-01 -2.864 0.004193 **
## Age_under19_Ptge 1.173e-01 8.738e-02 1.343 0.179459
## Age_19_65_pct 2.984e-01 7.676e-02 3.887 0.000102 ***
## Age_over65_pct 7.389e-02 7.461e-02 0.990 0.322036
## WomanPopulationPtge 1.629e-01 3.696e-02 4.408 1.06e-05 ***
## ForeignersPtge -1.013e-01 2.825e-02 -3.585 0.000340 ***
## SameComAutonPtge -7.710e-03 2.153e-02 -0.358 0.720299
## SameComAutonDiffProvPtge -1.263e-03 3.199e-02 -0.039 0.968516
## DifComAutonPtge 2.528e-01 2.772e-02 9.120 < 2e-16 ***
## UnemployLess25_Ptge 8.431e-02 1.874e-02 4.500 6.92e-06 ***
## Unemploy25_40_Ptge 1.516e-02 1.049e-02 1.444 0.148681
## UnemployMore40_Ptge 1.645e-02 9.885e-03 1.664 0.096146 .
## AgricultureUnemploymentPtge 1.145e-01 1.688e-02 6.786 1.26e-11 ***
## IndustryUnemploymentPtge 4.655e-02 1.519e-02 3.065 0.002189 **
## ConstructionUnemploymentPtge 5.959e-02 1.363e-02 4.373 1.25e-05 ***
## ServicesUnemploymentPtge 3.285e-02 8.874e-03 3.702 0.000216 ***
## totalEmpresas -3.284e-03 2.112e-03 -1.555 0.120045
## Industria -1.908e-03 1.535e-02 -0.124 0.901036
## Construccion -3.014e-02 9.730e-03 -3.098 0.001956 **
## Hosteleria -3.753e-03 4.561e-03 -0.823 0.410617
## Servicios -2.559e-03 6.020e-03 -0.425 0.670759
## ActividadPpalHosteleria 2.282e+00 5.433e-01 4.199 2.71e-05 ***
## ActividadPpalOtro -1.016e+00 6.947e-01 -1.462 0.143767
## inmuebles 2.953e-04 1.886e-04 1.566 0.117370
## Pob2010 2.584e-05 1.552e-04 0.166 0.867791
## SUPERFICIE 4.304e-05 2.917e-05 1.476 0.140078
## DensidadBaja -7.591e-01 4.884e-01 -1.554 0.120167
## PobChange_pct -1.030e-02 1.570e-02 -0.656 0.511538
## PersonasInmueble -2.355e-01 3.806e-01 -0.619 0.536034
## Explotaciones 1.963e-03 1.675e-03 1.172 0.241348
## prop_missings -7.786e+00 2.426e+00 -3.210 0.001336 **
## aleatorio -6.894e-01 4.314e-01 -1.598 0.110081
## aleatorio2 1.706e-01 4.305e-01 0.396 0.691941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.969 on 6442 degrees of freedom
## Multiple R-squared: 0.6346, Adjusted R-squared: 0.6316
## F-statistic: 207.2 on 54 and 6442 DF, p-value: < 2.2e-16Rsq(modelo1,"varObjCont",data_test) ## [1] 0.6335737# Pruebo con las variables generadas por mí
data_train <- todo[trainIndex, -c(22:25, 2)]
data_test <- todo[-trainIndex, -c(22:25, 2)]
modelo2<-lm(varObjCont~.,data=data_train)
summary(modelo2)##
## Call:
## lm(formula = varObjCont ~ ., data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.085 -9.764 -0.586 8.828 59.558
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.049e+01 1.224e+01 -2.492 0.01274 *
## CodigoProvincia 1.219e-02 1.183e-02 1.030 0.30306
## Population -1.082e-04 2.078e-04 -0.520 0.60276
## TotalCensus 2.960e-05 2.577e-04 0.115 0.90856
## Age_0_4_Ptge -5.066e-01 1.663e-01 -3.047 0.00232 **
## Age_under19_Ptge 4.623e-01 1.165e-01 3.967 7.37e-05 ***
## Age_19_65_pct 4.164e-01 1.032e-01 4.036 5.51e-05 ***
## Age_over65_pct 1.533e-01 1.003e-01 1.529 0.12628
## WomanPopulationPtge 9.212e-02 4.899e-02 1.880 0.06009 .
## ForeignersPtge -1.801e-01 3.733e-02 -4.823 1.44e-06 ***
## SameComAutonPtge 7.846e-02 2.837e-02 2.766 0.00569 **
## SameComAutonDiffProvPtge -5.054e-01 3.865e-02 -13.076 < 2e-16 ***
## DifComAutonPtge 4.649e-01 3.646e-02 12.752 < 2e-16 ***
## UnemployLess25_Ptge 2.960e-01 2.477e-02 11.946 < 2e-16 ***
## Unemploy25_40_Ptge 3.833e-02 1.405e-02 2.729 0.00638 **
## UnemployMore40_Ptge -2.670e-03 1.324e-02 -0.202 0.84024
## AgricultureUnemploymentPtge 2.118e-01 2.227e-02 9.509 < 2e-16 ***
## IndustryUnemploymentPtge -1.006e-01 2.002e-02 -5.024 5.19e-07 ***
## ConstructionUnemploymentPtge 7.826e-02 1.814e-02 4.313 1.63e-05 ***
## ServicesUnemploymentPtge 1.241e-02 1.182e-02 1.050 0.29385
## totalEmpresas -6.189e-03 2.496e-03 -2.480 0.01317 *
## ActividadPpalHosteleria 4.941e+00 7.116e-01 6.944 4.17e-12 ***
## ActividadPpalOtro 1.279e+01 6.113e+00 2.093 0.03638 *
## inmuebles 4.285e-04 2.467e-04 1.737 0.08239 .
## Pob2010 4.358e-05 2.067e-04 0.211 0.83298
## SUPERFICIE 2.412e-04 3.859e-05 6.250 4.37e-10 ***
## DensidadBaja 1.334e+00 6.504e-01 2.051 0.04029 *
## PobChange_pct 6.298e-03 2.100e-02 0.300 0.76430
## PersonasInmueble 4.608e-01 5.057e-01 0.911 0.36220
## Explotaciones 1.460e-03 2.213e-03 0.660 0.50946
## ActPpal.bin1 9.581e+00 6.116e+00 1.566 0.11729
## CCAA.groupednorteños 3.222e-01 5.024e-01 0.641 0.52130
## CCAA.groupedpijeños -1.372e+01 5.572e-01 -24.624 < 2e-16 ***
## CCAA.groupedsureños 2.740e+00 5.409e-01 5.066 4.17e-07 ***
## prop_missings -7.114e+00 3.208e+00 -2.217 0.02663 *
## aleatorio -3.434e-01 5.793e-01 -0.593 0.55335
## aleatorio2 -3.242e-03 5.782e-01 -0.006 0.99553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.41 on 6460 degrees of freedom
## Multiple R-squared: 0.3374, Adjusted R-squared: 0.3337
## F-statistic: 91.39 on 36 and 6460 DF, p-value: < 2.2e-16Rsq(modelo2,"varObjCont",data_test) ## [1] 0.3419729# Nos fijamos en la importancia de las variables. Podemos sacar un gráfico que muestra lo que se pierde en R2 en train al quitarlas del modelomodelEffectSizes(modelo2)
par(mfrow=c(1,2), mar=c(4.1, 14.1, 4.1, 6.1))
barplot(sort(modelEffectSizes(modelo1)$Effects[-1,4],decreasing =T),horiz=T,las=2,main="Importancia de las variables (R2)")## lm(formula = varObjCont ~ ., data = data_train)
##
## Coefficients
## SSR df pEta-sqr dR-sqr
## (Intercept) 563.8159 1 0.0009 NA
## CodigoProvincia 483.6273 1 0.0008 0.0003
## CCAA 671787.9965 18 0.5120 0.3833
## Population 33.2938 1 0.0001 0.0000
## TotalCensus 587.8858 1 0.0009 0.0003
## Age_0_4_Ptge 815.4237 1 0.0013 0.0005
## Age_under19_Ptge 179.1482 1 0.0003 0.0001
## Age_19_65_pct 1501.6910 1 0.0023 0.0009
## Age_over65_pct 97.4819 1 0.0002 0.0001
## WomanPopulationPtge 1930.8135 1 0.0030 0.0011
## ForeignersPtge 1277.3175 1 0.0020 0.0007
## SameComAutonPtge 12.7434 1 0.0000 0.0000
## SameComAutonDiffProvPtge 0.1548 1 0.0000 0.0000
## DifComAutonPtge 8267.3733 1 0.0127 0.0047
## UnemployLess25_Ptge 2012.4035 1 0.0031 0.0011
## Unemploy25_40_Ptge 207.3490 1 0.0003 0.0001
## UnemployMore40_Ptge 275.2233 1 0.0004 0.0002
## AgricultureUnemploymentPtge 4576.9445 1 0.0071 0.0026
## IndustryUnemploymentPtge 933.4113 1 0.0015 0.0005
## ConstructionUnemploymentPtge 1900.6025 1 0.0030 0.0011
## ServicesUnemploymentPtge 1362.1709 1 0.0021 0.0008
## totalEmpresas 240.2615 1 0.0004 0.0001
## Industria 1.5370 1 0.0000 0.0000
## Construccion 953.9576 1 0.0015 0.0005
## Hosteleria 67.2968 1 0.0001 0.0000
## Servicios 17.9628 1 0.0000 0.0000
## ActividadPpal 4908.2890 2 0.0076 0.0028
## inmuebles 243.7747 1 0.0004 0.0001
## Pob2010 2.7544 1 0.0000 0.0000
## SUPERFICIE 216.4345 1 0.0003 0.0001
## Densidad 240.1032 1 0.0004 0.0001
## PobChange_pct 42.8336 1 0.0001 0.0000
## PersonasInmueble 38.0643 1 0.0001 0.0000
## Explotaciones 136.4574 1 0.0002 0.0001
## prop_missings 1023.8263 1 0.0016 0.0006
## aleatorio 253.8170 1 0.0004 0.0001
## aleatorio2 15.6048 1 0.0000 0.0000
##
## Sum of squared errors (SSE): 640266.4
## Sum of squared total (SST): 1752458.9barplot(sort(modelEffectSizes(modelo2)$Effects[-1,4],decreasing =T),horiz=T,las=2,main="Importancia de las variables (R2)")## lm(formula = varObjCont ~ ., data = data_train)
##
## Coefficients
## SSR df pEta-sqr dR-sqr
## (Intercept) 1115.8805 1 0.0010 NA
## CodigoProvincia 190.6820 1 0.0002 0.0001
## Population 48.6898 1 0.0000 0.0000
## TotalCensus 2.3714 1 0.0000 0.0000
## Age_0_4_Ptge 1668.7240 1 0.0014 0.0010
## Age_under19_Ptge 2828.0957 1 0.0024 0.0016
## Age_19_65_pct 2927.4973 1 0.0025 0.0017
## Age_over65_pct 420.2711 1 0.0004 0.0002
## WomanPopulationPtge 635.5771 1 0.0005 0.0004
## ForeignersPtge 4181.8396 1 0.0036 0.0024
## SameComAutonPtge 1375.2438 1 0.0012 0.0008
## SameComAutonDiffProvPtge 30732.3158 1 0.0258 0.0175
## DifComAutonPtge 29229.2950 1 0.0246 0.0167
## UnemployLess25_Ptge 25652.4178 1 0.0216 0.0146
## Unemploy25_40_Ptge 1338.2407 1 0.0012 0.0008
## UnemployMore40_Ptge 7.3046 1 0.0000 0.0000
## AgricultureUnemploymentPtge 16253.2594 1 0.0138 0.0093
## IndustryUnemploymentPtge 4537.1402 1 0.0039 0.0026
## ConstructionUnemploymentPtge 3343.6338 1 0.0029 0.0019
## ServicesUnemploymentPtge 198.0928 1 0.0002 0.0001
## totalEmpresas 1105.2863 1 0.0010 0.0006
## ActividadPpal 9080.7903 2 0.0078 0.0052
## inmuebles 542.4684 1 0.0005 0.0003
## Pob2010 7.9937 1 0.0000 0.0000
## SUPERFICIE 7021.1868 1 0.0060 0.0040
## Densidad 756.2107 1 0.0007 0.0004
## PobChange_pct 16.1604 1 0.0000 0.0000
## PersonasInmueble 149.2516 1 0.0001 0.0001
## Explotaciones 78.2298 1 0.0001 0.0000
## ActPpal.bin 441.0395 1 0.0004 0.0003
## CCAA.grouped 153925.5690 3 0.1170 0.0878
## prop_missings 883.7963 1 0.0008 0.0005
## aleatorio 63.1577 1 0.0001 0.0000
## aleatorio2 0.0057 1 0.0000 0.0000
##
## Sum of squared errors (SSE): 1161132.8
## Sum of squared total (SST): 1752458.9
par(mfrow=c(1,1))
data_train<- todo[trainIndex, -c(34,35)]
data_test<- todo[-trainIndex, -c(34,35)]Bueno aquí podemos observar que la pérdida de información para la regresión con las variables generadas por mí a lo largo del proceso es excesiva. El \(R^2\) pasa a ser prácticamente la mitad ya que la importancia de CCAA con respecto al resto es altísima y seguramente al avanzar veamos como es prácticamente la variable que explica en mayor medida los resultados. De ahora en adelante el proceso en continuará con las variables originales, iré intentando mejorar el desempeño del modelo y reducir el número de variables implicadas. Debo comentar que los valores de \(R^2\) sobre los datos train y test es muy similar, así que en este sentido y pese a haber implicado una cantidad muy elevada de variables no hemos sobreajustado el modelo.
3.1.2 Modelos intermedios
Pruebo ahora con algunos modelos cuyas configuraciones creo que pueden ser competitivas.
#Pruebo un modelo con menos variables. Sobreescribo el modelo2 anterior
cramer.con<- c("CCAA", "ActividadPpal", "UnemployLess25_Ptge", "SameComAutonDiffProvPtge",
"UnemployMore40_Ptge", "AgricultureUnemploymentPtge", "Age_under19_Ptge",
"Age_over65_pct", "CodigoProvincia", "Age_0_4_Ptge")
frml<- as.formula(paste("varObjCont", paste(cramer.con, collapse="+"), sep="~"))
modelo2<-lm(frml, data=data_train)
summary(modelo2)##
## Call:
## lm(formula = frml, data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.200 -5.966 -0.435 5.834 57.987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.614475 1.376176 41.866 < 2e-16 ***
## CCAAAragón -9.867380 0.676436 -14.587 < 2e-16 ***
## CCAAAsturias -3.817137 1.388471 -2.749 0.00599 **
## CCAABaleares -9.065714 1.472604 -6.156 7.90e-10 ***
## CCAACanarias -15.697205 1.342980 -11.688 < 2e-16 ***
## CCAACantabria -15.951512 1.267625 -12.584 < 2e-16 ***
## CCAACastillaLeón -19.211423 0.565147 -33.994 < 2e-16 ***
## CCAACastillaMancha -9.475128 0.593772 -15.958 < 2e-16 ***
## CCAACataluña -42.922112 0.614327 -69.868 < 2e-16 ***
## CCAACeuta -22.620147 10.254456 -2.206 0.02743 *
## CCAAComValenciana -30.692838 0.668067 -45.943 < 2e-16 ***
## CCAAExtremadura -1.813571 0.736126 -2.464 0.01378 *
## CCAAGalicia -30.591377 0.809948 -37.770 < 2e-16 ***
## CCAAMadrid -15.128887 1.019338 -14.842 < 2e-16 ***
## CCAAMelilla -20.970561 10.258337 -2.044 0.04097 *
## CCAAMurcia -18.799859 1.792784 -10.486 < 2e-16 ***
## CCAANavarra -11.961098 0.852148 -14.036 < 2e-16 ***
## CCAAPaísVasco -21.366642 0.857971 -24.904 < 2e-16 ***
## CCAARioja -14.450775 1.045415 -13.823 < 2e-16 ***
## ActividadPpalHosteleria 2.104048 0.534064 3.940 8.24e-05 ***
## ActividadPpalOtro -0.892443 0.571823 -1.561 0.11864
## UnemployLess25_Ptge 0.096167 0.018022 5.336 9.82e-08 ***
## SameComAutonDiffProvPtge -0.014498 0.032493 -0.446 0.65549
## UnemployMore40_Ptge 0.036224 0.005770 6.278 3.65e-10 ***
## AgricultureUnemploymentPtge 0.069474 0.014638 4.746 2.12e-06 ***
## Age_under19_Ptge -0.096835 0.050322 -1.924 0.05436 .
## Age_over65_pct -0.219831 0.022218 -9.894 < 2e-16 ***
## CodigoProvincia 0.017399 0.009871 1.763 0.07800 .
## Age_0_4_Ptge -0.380630 0.123252 -3.088 0.00202 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.24 on 6468 degrees of freedom
## Multiple R-squared: 0.6133, Adjusted R-squared: 0.6117
## F-statistic: 366.4 on 28 and 6468 DF, p-value: < 2.2e-16Rsq(modelo2,"varObjCont",data_test)## [1] 0.6129429#Pruebo un modelo con (más o menos) la mitad de variables, basánndome en la importancia su variables
impor.con<- c("CCAA", "DifComAutonPtge", "ActividadPpal", "AgricultureUnemploymentPtge",
"UnemployLess25_Ptge", "ConstructionUnemploymentPtge")
frml<- as.formula(paste("varObjCont", paste(impor.con, collapse="+"), sep="~"))
modelo3<-lm(frml,data=data_train)
summary(modelo3)##
## Call:
## lm(formula = frml, data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.112 -5.776 -0.288 5.690 52.765
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.44473 0.71481 67.773 < 2e-16 ***
## CCAAAragón -11.76531 0.64896 -18.130 < 2e-16 ***
## CCAAAsturias -4.55315 1.35343 -3.364 0.000772 ***
## CCAABaleares -11.37601 1.45658 -7.810 6.62e-15 ***
## CCAACanarias -13.98393 1.31393 -10.643 < 2e-16 ***
## CCAACantabria -17.06383 1.23402 -13.828 < 2e-16 ***
## CCAACastillaLeón -20.39346 0.53460 -38.147 < 2e-16 ***
## CCAACastillaMancha -13.06049 0.61680 -21.175 < 2e-16 ***
## CCAACataluña -43.37001 0.59761 -72.573 < 2e-16 ***
## CCAACeuta -25.30136 10.17395 -2.487 0.012912 *
## CCAAComValenciana -31.32625 0.66364 -47.204 < 2e-16 ***
## CCAAExtremadura -3.56423 0.71695 -4.971 6.82e-07 ***
## CCAAGalicia -30.76833 0.77384 -39.761 < 2e-16 ***
## CCAAMadrid -16.62150 1.00256 -16.579 < 2e-16 ***
## CCAAMelilla -23.22358 10.17251 -2.283 0.022464 *
## CCAAMurcia -19.03084 1.77026 -10.750 < 2e-16 ***
## CCAANavarra -13.44646 0.83724 -16.060 < 2e-16 ***
## CCAAPaísVasco -22.65009 0.84680 -26.748 < 2e-16 ***
## CCAARioja -17.06323 1.04175 -16.379 < 2e-16 ***
## DifComAutonPtge 0.27327 0.01795 15.225 < 2e-16 ***
## ActividadPpalHosteleria 2.54639 0.52902 4.813 1.52e-06 ***
## ActividadPpalOtro -0.96821 0.52059 -1.860 0.062952 .
## AgricultureUnemploymentPtge 0.10643 0.01460 7.292 3.43e-13 ***
## UnemployLess25_Ptge 0.10612 0.01748 6.070 1.35e-09 ***
## ConstructionUnemploymentPtge 0.05272 0.01138 4.631 3.70e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.16 on 6472 degrees of freedom
## Multiple R-squared: 0.6188, Adjusted R-squared: 0.6174
## F-statistic: 437.7 on 24 and 6472 DF, p-value: < 2.2e-16Rsq(modelo3,"varObjCont",data_test) ## [1] 0.6230918#Pruebo con una interaccion sobre el anterior
inter1.con<- c("CCAA", "DifComAutonPtge", "ActividadPpal", "AgricultureUnemploymentPtge",
"UnemployLess25_Ptge", "ConstructionUnemploymentPtge", "CCAA:ActividadPpal")
frml<- as.formula(paste("varObjCont", paste(inter1.con, collapse="+"), sep="~"))
modelo4<-lm(frml,data=data_train)
summary(modelo4)##
## Call:
## lm(formula = frml, data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.277 -5.553 -0.245 5.492 51.679
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error t value
## (Intercept) 41.97278 1.64922 25.450
## CCAAAragón -2.96408 3.23078 -0.917
## CCAAAsturias 1.44578 6.02670 0.240
## CCAABaleares -1.67458 2.44043 -0.686
## CCAACanarias -4.30085 3.72552 -1.154
## CCAACantabria -8.18989 4.41706 -1.854
## CCAACastillaLeón -10.51932 2.32229 -4.530
## CCAACastillaMancha -9.51406 2.59754 -3.663
## CCAACataluña -35.09041 1.81056 -19.381
## CCAACeuta -25.11012 10.06585 -2.495
## CCAAComValenciana -26.84497 2.00800 -13.369
## CCAAExtremadura -7.96741 5.28333 -1.508
## CCAAGalicia -21.85057 2.69293 -8.114
## CCAAMadrid -12.46682 2.09441 -5.952
## CCAAMelilla -23.09025 10.06438 -2.294
## CCAAMurcia -13.31971 6.02611 -2.210
## CCAANavarra -2.47639 3.23359 -0.766
## CCAAPaísVasco -10.34221 2.73537 -3.781
## CCAARioja -10.29323 5.28952 -1.946
## DifComAutonPtge 0.25377 0.01806 14.052
## ActividadPpalHosteleria 9.28927 1.71370 5.421
## ActividadPpalOtro 7.36130 1.81019 4.067
## AgricultureUnemploymentPtge 0.09885 0.01458 6.780
## UnemployLess25_Ptge 0.10170 0.01734 5.864
## ConstructionUnemploymentPtge 0.05128 0.01133 4.525
## CCAAAragón:ActividadPpalHosteleria -7.77634 3.48909 -2.229
## CCAAAsturias:ActividadPpalHosteleria -7.20352 6.23227 -1.156
## CCAABaleares:ActividadPpalHosteleria -12.61596 3.31600 -3.805
## CCAACanarias:ActividadPpalHosteleria -10.24630 3.98499 -2.571
## CCAACantabria:ActividadPpalHosteleria -10.36316 4.70154 -2.204
## CCAACastillaLeón:ActividadPpalHosteleria -5.75374 2.50478 -2.297
## CCAACastillaMancha:ActividadPpalHosteleria -4.74443 2.72365 -1.742
## CCAACataluña:ActividadPpalHosteleria -8.15664 2.00743 -4.063
## CCAACeuta:ActividadPpalHosteleria NA NA NA
## CCAAComValenciana:ActividadPpalHosteleria -6.28080 2.19614 -2.860
## CCAAExtremadura:ActividadPpalHosteleria 3.59381 5.37661 0.668
## CCAAGalicia:ActividadPpalHosteleria -9.78888 2.82434 -3.466
## CCAAMadrid:ActividadPpalHosteleria -3.65886 2.63779 -1.387
## CCAAMelilla:ActividadPpalHosteleria NA NA NA
## CCAAMurcia:ActividadPpalHosteleria -5.36389 6.31863 -0.849
## CCAANavarra:ActividadPpalHosteleria -9.36800 3.54165 -2.645
## CCAAPaísVasco:ActividadPpalHosteleria -9.36853 2.94701 -3.179
## CCAARioja:ActividadPpalHosteleria -5.72642 5.83891 -0.981
## CCAAAragón:ActividadPpalOtro -10.51847 3.35723 -3.133
## CCAAAsturias:ActividadPpalOtro -3.88546 6.60659 -0.588
## CCAABaleares:ActividadPpalOtro -17.10590 5.64084 -3.033
## CCAACanarias:ActividadPpalOtro -10.39324 10.74436 -0.967
## CCAACantabria:ActividadPpalOtro -8.61198 4.84419 -1.778
## CCAACastillaLeón:ActividadPpalOtro -11.96851 2.46414 -4.857
## CCAACastillaMancha:ActividadPpalOtro -4.19563 2.74249 -1.530
## CCAACataluña:ActividadPpalOtro -10.58630 2.03982 -5.190
## CCAACeuta:ActividadPpalOtro NA NA NA
## CCAAComValenciana:ActividadPpalOtro -3.04611 2.28846 -1.331
## CCAAExtremadura:ActividadPpalOtro 3.84766 5.39260 0.714
## CCAAGalicia:ActividadPpalOtro -5.48262 3.49935 -1.567
## CCAAMadrid:ActividadPpalOtro -1.68920 2.77357 -0.609
## CCAAMelilla:ActividadPpalOtro NA NA NA
## CCAAMurcia:ActividadPpalOtro -10.05774 9.35140 -1.076
## CCAANavarra:ActividadPpalOtro -13.37857 3.42157 -3.910
## CCAAPaísVasco:ActividadPpalOtro -19.13164 3.05608 -6.260
## CCAARioja:ActividadPpalOtro -8.35738 5.43095 -1.539
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## CCAAAragón 0.358941
## CCAAAsturias 0.810418
## CCAABaleares 0.492622
## CCAACanarias 0.248368
## CCAACantabria 0.063764 .
## CCAACastillaLeón 6.01e-06 ***
## CCAACastillaMancha 0.000252 ***
## CCAACataluña < 2e-16 ***
## CCAACeuta 0.012635 *
## CCAAComValenciana < 2e-16 ***
## CCAAExtremadura 0.131596
## CCAAGalicia 5.82e-16 ***
## CCAAMadrid 2.78e-09 ***
## CCAAMelilla 0.021808 *
## CCAAMurcia 0.027117 *
## CCAANavarra 0.443803
## CCAAPaísVasco 0.000158 ***
## CCAARioja 0.051702 .
## DifComAutonPtge < 2e-16 ***
## ActividadPpalHosteleria 6.16e-08 ***
## ActividadPpalOtro 4.83e-05 ***
## AgricultureUnemploymentPtge 1.31e-11 ***
## UnemployLess25_Ptge 4.76e-09 ***
## ConstructionUnemploymentPtge 6.15e-06 ***
## CCAAAragón:ActividadPpalHosteleria 0.025864 *
## CCAAAsturias:ActividadPpalHosteleria 0.247789
## CCAABaleares:ActividadPpalHosteleria 0.000143 ***
## CCAACanarias:ActividadPpalHosteleria 0.010156 *
## CCAACantabria:ActividadPpalHosteleria 0.027545 *
## CCAACastillaLeón:ActividadPpalHosteleria 0.021645 *
## CCAACastillaMancha:ActividadPpalHosteleria 0.081567 .
## CCAACataluña:ActividadPpalHosteleria 4.90e-05 ***
## CCAACeuta:ActividadPpalHosteleria NA
## CCAAComValenciana:ActividadPpalHosteleria 0.004251 **
## CCAAExtremadura:ActividadPpalHosteleria 0.503892
## CCAAGalicia:ActividadPpalHosteleria 0.000532 ***
## CCAAMadrid:ActividadPpalHosteleria 0.165462
## CCAAMelilla:ActividadPpalHosteleria NA
## CCAAMurcia:ActividadPpalHosteleria 0.395968
## CCAANavarra:ActividadPpalHosteleria 0.008187 **
## CCAAPaísVasco:ActividadPpalHosteleria 0.001485 **
## CCAARioja:ActividadPpalHosteleria 0.326761
## CCAAAragón:ActividadPpalOtro 0.001738 **
## CCAAAsturias:ActividadPpalOtro 0.556473
## CCAABaleares:ActividadPpalOtro 0.002435 **
## CCAACanarias:ActividadPpalOtro 0.333420
## CCAACantabria:ActividadPpalOtro 0.075485 .
## CCAACastillaLeón:ActividadPpalOtro 1.22e-06 ***
## CCAACastillaMancha:ActividadPpalOtro 0.126100
## CCAACataluña:ActividadPpalOtro 2.17e-07 ***
## CCAACeuta:ActividadPpalOtro NA
## CCAAComValenciana:ActividadPpalOtro 0.183212
## CCAAExtremadura:ActividadPpalOtro 0.475557
## CCAAGalicia:ActividadPpalOtro 0.117222
## CCAAMadrid:ActividadPpalOtro 0.542524
## CCAAMelilla:ActividadPpalOtro NA
## CCAAMurcia:ActividadPpalOtro 0.282176
## CCAANavarra:ActividadPpalOtro 9.32e-05 ***
## CCAAPaísVasco:ActividadPpalOtro 4.09e-10 ***
## CCAARioja:ActividadPpalOtro 0.123892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.05 on 6440 degrees of freedom
## Multiple R-squared: 0.629, Adjusted R-squared: 0.6258
## F-statistic: 195 on 56 and 6440 DF, p-value: < 2.2e-16Rsq(modelo4,"varObjCont",data_test)## [1] 0.6281302El segundo (que reemplaza al modelo fallido con las variables imputadas) lo he definido mirando la importancia de las variables en el gráfico de la V de Cramer y seleccionando las 10 variables más relevantes, el tercero siguiendo los gráficos de relevancia de las variables según \(R^2\) para el primer modelo y seleccionando las 5 más importantes y el cuarto definiendo una interacción adicional entre las dos variables cualitativas (CCAA y ActividadPpal) involucradas en el tercer modelo. Mirando el summary de cada modelo podemos ver que en el segundo hay varias variables poco significativas, lo que nos indica que uno no debe sobrestimar el poder descriptivo de la V de Cramer. En el tercero prácticamente todas tienen una significación alta, salvo alguna de las ‘dummies’. En el último la interacción ha introducido muchas variables y vuelve a aumentar el número de variables poco significativas. Si miramos el valor de \(R^2\) podemos ver que la variación entre los modelos se mueve en una franja inferior a 30 milésimas.
3.1.3 Validación cruzada
Hago validación cruzada sobre los modelos para seleccionar el más competitivo.
# Hago validación repetida para ver qué modelo es mejor
modelo1VC <- train(formula(modelo1),
data = todo,method = "lm",
trControl = trainControl(method="repeatedcv", number=5, repeats=20, returnResamp="all")
)
modelo2VC <- train(formula(modelo2),
data = todo,method = "lm",
trControl = trainControl(method="repeatedcv", number=5, repeats=20, returnResamp="all")
)
modelo3VC <- train(formula(modelo3),
data = todo,method = "lm",
trControl = trainControl(method="repeatedcv", number=5, repeats=20, returnResamp="all")
)
modelo4VC <- train(formula(modelo4),
data = todo,method = "lm",
trControl = trainControl(method="repeatedcv", number=5, repeats=20, returnResamp="all")
)
results<-data.frame(rbind(modelo1VC$resample,modelo2VC$resample,modelo3VC$resample,modelo4VC$resample),modelo=c(rep(1,100),rep(2,100),rep(3,100),rep(4,100)))
comparativa<- data.frame(row.names = c("modelo1", "modelo2", "modelo3", "modelo4"))
comparativa$mean_R2<- aggregate(Rsquared~modelo, data = results, mean)$Rsquared
comparativa$std_R2<- aggregate(Rsquared~modelo, data = results, sd)$Rsquared
comparativa$variables<- c(length(coef(modelo1)), length(coef(modelo2)), length(coef(modelo3)),
length(coef(modelo4)))
comparativa$test_R2<- c(Rsq(modelo1,"varObjCont",data_test), Rsq(modelo2,"varObjCont",data_test),
Rsq(modelo3,"varObjCont",data_test), Rsq(modelo4,"varObjCont",data_test))
# Elijo el modelo 3 por el bajo número de variables y la escasa diferencia con el R2 de los demás| mean_R2 | std_R2 | variables | test_R2 | |
|---|---|---|---|---|
| modelo1 | 0,629 | 0,014 | 55 | 0,634 |
| modelo2 | 0,611 | 0,014 | 29 | 0,613 |
| modelo3 | 0,617 | 0,013 | 25 | 0,623 |
| modelo4 | 0,625 | 0,013 | 61 | 0,628 |
## lm(formula = frml, data = data_train)
##
## Coefficients
## SSR df pEta-sqr dR-sqr
## (Intercept) 474126.578 1 0.4151 NA
## CCAA 827203.600 18 0.5532 0.4720
## DifComAutonPtge 23928.201 1 0.0346 0.0137
## ActividadPpal 11867.188 2 0.0175 0.0068
## AgricultureUnemploymentPtge 5488.217 1 0.0081 0.0031
## UnemployLess25_Ptge 3803.520 1 0.0057 0.0022
## ConstructionUnemploymentPtge 2214.045 1 0.0033 0.0013
##
## Sum of squared errors (SSE): 668059.8
## Sum of squared total (SST): 1752458.9
La validación cruzada sobre estos modelos construidos a mano nos permite ver que los resultados se mueven ya en una franja bastante estrecha. Pese a que el mejor en cuanto a \(R^2\) es el que hemos construído con todas las variables su número de variables es muy elevado con respecto al resto. Para mi gusto el más competitivo es el tercer modelo, construído eligiendo las 6 variables más significativas del modelo que las contenía todas, ya que con tan solo 25 permite un ajuste prácticamente igual de bueno que el del resto.
3.1.4 Transformación de las variables
Llegado este punto vamos a buscar en las variables las transformaciones fundamentales que puedan mejorar la predicción de las variables objetivos. Como las transformaciones no serán comunes a ambas variables a partir de aqui separaremos los datos que emplearemos para cada la predicción de cada variable objetivo.
input <- input[, -c(34,35)]
input_cont<-cbind(input,Transf_Auto(Filter(is.numeric, input[,-c(35,36)]),varObjCont))
input_bin<-cbind(input,Transf_Auto(Filter(is.numeric, input[,-c(35,36)]),varObjBin))
saveRDS(data.frame(input, varObjBin, varObjCont),"todoNoTrans")
saveRDS(data.frame(input_bin,varObjBin),"todo_bin")
saveRDS(data.frame(input_cont,varObjCont),"todo_cont")
saveRDS(modelo3, 'modeloManual')
saveRDS(modelo1, 'modeloFull')3.1.5 Selección de Variables
Ahora voy a emplear los métodos de selección de variables estudiados en clase, aunque usaré de forma predominante el Stepwise, probaremos también con Lasso y un un método consistente en la incorporación al azar de una submuestra de nuestros datos sobre la que aplicar el método Stepwise (aquí se ha comentado debido al tiempo de ejecución, en el código aparece sin comentar y lo ejecuté para obtener las fórmulas de los modelos aleatorios). Para elegir el mejor modelo usaremos los criterios de información también vistos: \[ \text{Akaike: }n\ln\left(\frac{SSE}{n}\right) + 2 \cdot p\\ \text{Schwarz: }n\ln\left(\frac{SSE}{n}\right) + p \cdot \ln(n)\\ \text{Siendo }n \text{ el número de observaciones y }p \text{ el número de parámetros} \]
Con ellos determinaré las variables que escogeré en cada modelo. En este caso la introduciré también interacciones (todas las posibles) entre las variables y las transformaciones de las variables.
## Modelos básicos ----
datos<-readRDS("todo_cont")
set.seed(123456)
trainIndex <- createDataPartition(datos$varObjCont, p=0.8, list=FALSE)
data_train <- datos[trainIndex,]
data_test <- datos[-trainIndex,]
null<-lm(varObjCont~1, data=data_train) #Modelo minimo
full<-readRDS("modeloFull")
modeloManual<-readRDS("modeloManual")
# Seleccion de variables "clásica"
modeloStepAIC<-step(null, scope=list(lower=null, upper=full), direction="both")
modeloBackAIC<-step(full, scope=list(lower=null, upper=full), direction="backward")
modeloStepBIC<-step(null, scope=list(lower=null, upper=full), direction="both",k=log(nrow(data_train)))
modeloBackBIC<-step(full, scope=list(lower=null, upper=full), direction="backward",k=log(nrow(data_train)))
## Modelos con interacciones ----
formInt<-formulaInteracciones(datos[,c(1:36,68)],37)#Set original la objetivo está en la columna 68
fullInt<-lm(formInt, data=data_train)
modeloStepAIC_int<-step(null, scope=list(lower=null, upper=fullInt), direction="both")
summary(modeloStepAIC_int) # Fallido
modeloStepBIC_int<-step(null, scope=list(lower=null, upper=fullInt), direction="both",k=log(nrow(data_train)))
## Modelos con transformadas ----
fullT<-lm(varObjCont~., data=data_train)
modeloStepAIC_trans<-step(null, scope=list(lower=null, upper=fullT), direction="both")
modeloStepBIC_trans<-step(null, scope=list(lower=null, upper=fullT), direction="both",k=log(nrow(data_train)))
## Modelos con todas ----
formIntT<-formulaInteracciones(datos,68)
fullIntT<-lm(formIntT, data=data_train)
modeloStepAIC_transInt<-step(null, scope=list(lower=null, upper=fullIntT), direction="both")
modeloStepBIC_transInt<-step(null, scope=list(lower=null, upper=fullIntT), direction="both", k=log(nrow(data_train)))
# # Seleccion aleatoria ----
# rep<-100
# prop<-0.7
# modelosGenerados<-c()
# for (i in 1:rep){
# set.seed(12345+i)
# subsample<-data_train[sample(1:nrow(data_train),prop*nrow(data_train),replace = T),]
# full<-lm(formIntT,data=subsample)
# null<-lm(varObjCont~1,data=subsample)
# modeloAux<-step(null,scope=list(lower=null,upper=full),direction="both",trace=0,k=log(nrow(subsample)))
# modelosGenerados<-c(modelosGenerados,paste(sort(unlist(strsplit(as.character(formula(modeloAux))[3]," [+] "))),
# collapse = "+"))
# }
# 2 únicos modelos que aparezcan más de una vez
modeloAleat1.frml<- as.formula("varObjCont ~ Age_0_4_Ptge+CCAA+CCAA:raiz4DifComAutonPtge+CodigoProvincia+ForeignersPtge+inmuebles+logxConstructionUnemploymentPtge+logxinmuebles+logxServicesUnemploymentPtge+raiz4DifComAutonPtge+sqrtxAge_19_65_pct+sqrtxAgricultureUnemploymentPtge+sqrxCodigoProvincia")
modeloAleat2.frml<- as.formula("varObjCont ~ ActividadPpal+ActividadPpal:AgricultureUnemploymentPtge+ActividadPpal:ForeignersPtge+Age_19_65_pct+Age_over65_pct+Age_under19_Ptge+AgricultureUnemploymentPtge+CCAA+CCAA:raiz4DifComAutonPtge+ForeignersPtge+inmuebles+logxinmuebles+raiz4DifComAutonPtge+sqrtxAgricultureUnemploymentPtge")
modeloAleat1<-lm(modeloAleat1.frml, data=data_train)
modeloAleat2<-lm(modeloAleat2.frml, data=data_train)
# LASSO ----
# Lo hacemos sin interacciones pues, de lo contrario, puede coger interacciones y no las variables que las forman
y<- as.double(as.matrix(data_train[, 68]))
x<- model.matrix(varObjCont~., data=data_train[,c(1:37,68)])[,-1]#no cambiar el -1
set.seed(1712)
cv.lasso <- cv.glmnet(x,y,nfolds=5)
plot(cv.lasso)
betas<-coef(cv.lasso, s=cv.lasso$lambda.1se)Una vez definidos los modelos y realizada la criba de variables voy a usar ‘crossvalidation’ para compararlos de una forma justa.
## Validación cruzada ----
total<-c()
# Elimino el modelo que ha fallado
modelos<-sapply(list(modeloManual,modeloStepAIC,modeloStepBIC,modeloStepBIC_int,modeloStepAIC_trans,
modeloStepBIC_trans, modeloStepAIC_transInt, modeloStepBIC_transInt, modeloAleat1,
modeloAleat2), formula)
model.names = c("modeloManual", "modeloStepAIC","modeloStepBIC","modeloStepBIC_int","modeloStepAIC_trans",
"modeloStepBIC_trans", "modeloStepAIC_transInt", "modeloStepBIC_transInt", "modeloAleat1",
"modeloAleat2")
for (i in 1:length(modelos)){
set.seed(1712)
vcr<-train(modelos[[i]], data = data_train,
method = "lm",
trControl = trainControl(method="repeatedcv", number=5, repeats=20,
returnResamp="all")
)
total<-rbind(total,cbind(vcr$resample[,1:2],modelo=rep(model.names[i],
nrow(vcr$resample))))
}
set.seed(1712)
lassovcr <- train(varObjCont ~ ., data = data_train,
method = "glmnet",
tuneGrid=expand.grid(.alpha=1,.lambda=cv.lasso$lambda.1se),
trControl = trainControl(method="repeatedcv", number=5, repeats=20,
returnResamp="all")
)
total<-rbind(total,cbind(lassovcr$resample[,1:2],modelo=rep("Lasso",
nrow(vcr$resample))))| mean_R2 | std_R2 | variables | test_R2 | |
|---|---|---|---|---|
| modeloManual | 0,616 | 0,016 | 25 | 0,623 |
| Lasso | 0,624 | 0,014 | 0 | 0,629 |
| modeloStepBIC | 0,628 | 0,015 | 31 | 0,632 |
| modeloStepAIC | 0,630 | 0,015 | 37 | 0,632 |
| modeloStepBIC_trans | 0,637 | 0,015 | 32 | 0,640 |
| modeloStepBIC_int | 0,642 | 0,015 | 54 | 0,641 |
| modeloStepBIC_transInt | 0,649 | 0,015 | 50 | 0,642 |
| modeloStepAIC_trans | 0,638 | 0,015 | 48 | 0,643 |
| modeloAleat1 | 0,647 | 0,015 | 48 | 0,644 |
| modeloAleat2 | 0,644 | 0,015 | 52 | 0,645 |
| modeloStepAIC_transInt | 0,662 | 0,015 | 186 | 0,660 |
Por algún motivo el Lasso aparece como si usará 0 variables, lo bueno es que el gráfico podemos ver que el modelo óptimo según este método sería de entre 35 y 50 variables.
En el gráfico de cajas podemos ver como el mejor modelo es el que emplea tanto variables transformadas como interacciones y ha seleccionado las variables mediante el criterio Akaike. No obstante este modelo tiene 200 variables, como podemos ver en la tabla, luego ni se aproxima a ser el mejor modelo. Los demás modelos se mueven en el entorno de 30 - 60 variables, y en esencia existen dos opciones posibles:
- El modelo Aleat1 que emplea 48 variables y tiene un \(R^2\) sobre la partición de test de 0.644
- El modelo BICtrans que emplea 32 variables incluyendo algunas transformadas y tiene un \(R^2\) sobre la partición de test de 0.640.
Desde mi punto de vista la mejor opción de las dos es el modelo con menos variables, ya que la mejora no me parece suficientemente significativa, y menos mirando las desviaciones típicas de los modelos. Damos por ganador a: BICtrans.
3.2 Regresión Logística
Los modelos que involucraban las variables transformadas se dilataban excesivamente en el tiempo así que he prescindido de ellos. Esencialmente he repetido la parte manual del proceso realizado para la regresión lineal, centrándome en la compprensión de la selección del punto de corte, ya que este es el problema sobre el que vamos a trabajar en el módulo de Machine Learning.
3.2.1 Definición de modelos
# REGRESIÓN LOGÍSTICA ----
#Cargo los datos depurados (excluyo las variables transformadas)
todo<-readRDS("todo_bin")
todo<- todo[,c(1:36,68)]
#veo el reparto original. Compruebo que la variable objetivo tome valor 1 para el evento y 0 para el no evento
freq(todo$varObjBin) #ese ha de ser el error de referencia## n % val%
## 0 3076 37.9 37.9
## 1 5043 62.1 62.1#Hago la partici?n
set.seed(123456)
trainIndex <- createDataPartition(todo$varObjBin, p=0.8, list=FALSE)
data_train <- todo[trainIndex,]
data_test <- todo[-trainIndex,]
#pruebo un primer modelo sin las transformadas
modeloInicial<-glm(varObjBin~.,data=data_train[],family=binomial)
summary(modeloInicial)##
## Call:
## glm(formula = varObjBin ~ ., family = binomial, data = data_train[])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8954 -0.2238 0.3553 0.5291 3.6565
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.457e+00 2.139e+00 2.083 0.037226 *
## CodigoProvincia 1.578e-03 3.165e-03 0.498 0.618153
## CCAAAragón 1.783e+00 1.734e-01 10.286 < 2e-16 ***
## CCAAAsturias 9.469e-01 2.891e-01 3.275 0.001057 **
## CCAABaleares 7.873e-01 3.292e-01 2.391 0.016785 *
## CCAACanarias 1.543e+00 3.007e-01 5.130 2.90e-07 ***
## CCAACantabria 3.787e+00 4.533e-01 8.356 < 2e-16 ***
## CCAACastillaLeón 3.386e+00 1.578e-01 21.451 < 2e-16 ***
## CCAACastillaMancha 2.188e+00 1.598e-01 13.694 < 2e-16 ***
## CCAACataluña -5.486e+00 7.249e-01 -7.568 3.80e-14 ***
## CCAACeuta 1.331e+01 3.247e+02 0.041 0.967311
## CCAAComValenciana 3.460e+00 2.238e-01 15.458 < 2e-16 ***
## CCAAExtremadura 9.240e-01 1.743e-01 5.300 1.16e-07 ***
## CCAAGalicia 4.090e+00 3.403e-01 12.016 < 2e-16 ***
## CCAAMadrid 2.789e+00 3.014e-01 9.255 < 2e-16 ***
## CCAAMelilla 1.266e+01 3.247e+02 0.039 0.968910
## CCAAMurcia 2.834e+00 5.573e-01 5.084 3.69e-07 ***
## CCAANavarra 2.602e-01 2.144e-01 1.214 0.224817
## CCAAPaísVasco -2.584e+00 4.388e-01 -5.890 3.86e-09 ***
## CCAARioja 2.867e+00 2.833e-01 10.120 < 2e-16 ***
## Population -2.244e-05 4.812e-05 -0.466 0.640965
## TotalCensus -5.790e-05 5.710e-05 -1.014 0.310607
## Age_0_4_Ptge 4.182e-02 3.689e-02 1.133 0.257063
## Age_under19_Ptge -2.777e-02 2.346e-02 -1.184 0.236505
## Age_19_65_pct -5.240e-02 2.003e-02 -2.616 0.008902 **
## Age_over65_pct -1.308e-02 1.928e-02 -0.679 0.497316
## WomanPopulationPtge -2.923e-03 1.067e-02 -0.274 0.784227
## ForeignersPtge 2.971e-02 8.903e-03 3.337 0.000846 ***
## SameComAutonPtge -1.088e-02 6.955e-03 -1.564 0.117765
## SameComAutonDiffProvPtge 2.376e-02 1.171e-02 2.030 0.042394 *
## DifComAutonPtge -2.846e-02 8.425e-03 -3.378 0.000729 ***
## UnemployLess25_Ptge -9.393e-03 5.272e-03 -1.782 0.074824 .
## Unemploy25_40_Ptge -4.014e-04 3.230e-03 -0.124 0.901118
## UnemployMore40_Ptge -1.696e-03 3.053e-03 -0.555 0.578578
## AgricultureUnemploymentPtge -1.003e-02 4.880e-03 -2.056 0.039758 *
## IndustryUnemploymentPtge -7.894e-04 4.635e-03 -0.170 0.864745
## ConstructionUnemploymentPtge -1.839e-03 4.118e-03 -0.447 0.655191
## ServicesUnemploymentPtge -9.645e-04 2.746e-03 -0.351 0.725413
## totalEmpresas 8.817e-04 7.091e-04 1.243 0.213728
## Industria 1.143e-03 4.806e-03 0.238 0.811964
## Construccion 2.620e-03 2.942e-03 0.891 0.373171
## Hosteleria -8.362e-04 1.430e-03 -0.585 0.558821
## Servicios 4.891e-04 1.941e-03 0.252 0.801076
## ActividadPpalHosteleria -8.774e-01 2.042e-01 -4.296 1.74e-05 ***
## ActividadPpalOtro -6.024e-01 2.412e-01 -2.497 0.012523 *
## inmuebles 3.036e-06 5.671e-05 0.054 0.957314
## Pob2010 7.704e-05 4.470e-05 1.724 0.084783 .
## SUPERFICIE -1.607e-05 7.799e-06 -2.060 0.039386 *
## DensidadBaja 1.105e-01 1.454e-01 0.760 0.447485
## PobChange_pct 4.671e-03 4.412e-03 1.059 0.289791
## PersonasInmueble 5.306e-02 1.184e-01 0.448 0.653998
## Explotaciones -2.784e-04 4.706e-04 -0.592 0.554102
## prop_missings 3.011e+00 7.347e-01 4.098 4.17e-05 ***
## aleatorio -5.543e-02 1.283e-01 -0.432 0.665633
## aleatorio2 -5.389e-02 1.268e-01 -0.425 0.670783
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8620.2 on 6495 degrees of freedom
## Residual deviance: 4792.0 on 6441 degrees of freedom
## AIC: 4902
##
## Number of Fisher Scoring iterations: 11# Miro las variables más importantes y la curva ROC
par(mfrow=c(1,2), mar=c(4.1, 12, 4.1, 2.1))
impVariablesLog(modeloInicial,"varObjBin")
curvaRoc<- roc(data_test$varObjBin, predict(modeloInicial,data_test,type = "response"), direction="<")
plot(curvaRoc)
par(mfrow=c(1,1))
#Miro el gráfico V de Cramer para ver las variables más importantes
par(mar=c(12, 4.1, 4.1, 2.1))
graficoVcramer(todo, todo$varObjBin, "V de Cramer para Derecha") 
cramer.bin<- c("Age_under19_Ptge", "Age_over65_pct","PersonasInmueble",
"Age_0_4_Ptge", "CodigoProvincia", "Age_19_65_pct",
"PobChange_pct", "UnemployMore40_Ptge", "CCAA", "ActividadPpal")
frml<- as.formula(paste("varObjBin", paste(cramer.bin, collapse="+"), sep="~"))
modelo2<-glm(frml,data=data_train,family=binomial)
summary(modelo2)##
## Call:
## glm(formula = frml, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7955 -0.2268 0.3808 0.5130 3.6225
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.613777 1.906565 1.371 0.170395
## Age_under19_Ptge -0.017619 0.022482 -0.784 0.433215
## Age_over65_pct -0.007085 0.018554 -0.382 0.702565
## PersonasInmueble 0.038885 0.106230 0.366 0.714331
## Age_0_4_Ptge 0.064170 0.036281 1.769 0.076947 .
## CodigoProvincia 0.001533 0.003111 0.493 0.622145
## Age_19_65_pct -0.040523 0.019303 -2.099 0.035796 *
## PobChange_pct 0.004555 0.004300 1.059 0.289425
## UnemployMore40_Ptge -0.002099 0.001603 -1.309 0.190381
## CCAAAragón 1.780339 0.157444 11.308 < 2e-16 ***
## CCAAAsturias 0.834891 0.272659 3.062 0.002198 **
## CCAABaleares 0.979563 0.311217 3.148 0.001647 **
## CCAACanarias 1.848876 0.283141 6.530 6.58e-11 ***
## CCAACantabria 3.542812 0.442040 8.015 1.10e-15 ***
## CCAACastillaLeón 3.364670 0.147962 22.740 < 2e-16 ***
## CCAACastillaMancha 2.005234 0.137232 14.612 < 2e-16 ***
## CCAACataluña -5.205583 0.719214 -7.238 4.56e-13 ***
## CCAACeuta 13.528740 324.743741 0.042 0.966770
## CCAAComValenciana 3.515049 0.210178 16.724 < 2e-16 ***
## CCAAExtremadura 0.529013 0.159449 3.318 0.000907 ***
## CCAAGalicia 4.101512 0.331430 12.375 < 2e-16 ***
## CCAAMadrid 2.841224 0.282652 10.052 < 2e-16 ***
## CCAAMelilla 13.378031 324.743752 0.041 0.967140
## CCAAMurcia 3.195654 0.540141 5.916 3.29e-09 ***
## CCAANavarra 0.056227 0.187237 0.300 0.763951
## CCAAPaísVasco -2.601900 0.429632 -6.056 1.39e-09 ***
## CCAARioja 2.765279 0.261701 10.567 < 2e-16 ***
## ActividadPpalHosteleria -1.101589 0.192982 -5.708 1.14e-08 ***
## ActividadPpalOtro -1.024593 0.210448 -4.869 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8620.2 on 6495 degrees of freedom
## Residual deviance: 4891.9 on 6467 degrees of freedom
## AIC: 4949.9
##
## Number of Fisher Scoring iterations: 11par(mar=c(4.1, 12, 4.1, 2.1))
impVariablesLog(modelo2,"varObjBin")
# Selecciono las 10 variables más influyentes en el modelo Inicial
main.bin<- c("CCAA", "ActividadPpal", "prop_missings", "ForeignersPtge",
"DifComAutonPtge","PersonasInmueble", "Age_19_65_pct",
"Construccion", "Population", "SameComAutonDiffProvPtge")
frml<- as.formula(paste("varObjBin", paste(main.bin, collapse="+"), sep="~"))
modelo3<-glm(frml,data=data_train,family=binomial)
summary(modelo3)##
## Call:
## glm(formula = frml, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0040 -0.2206 0.3623 0.5270 3.6831
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.225e+00 4.557e-01 2.688 0.007187 **
## CCAAAragón 1.936e+00 1.512e-01 12.800 < 2e-16 ***
## CCAAAsturias 1.041e+00 2.748e-01 3.787 0.000152 ***
## CCAABaleares 8.948e-01 3.207e-01 2.790 0.005275 **
## CCAACanarias 1.770e+00 2.844e-01 6.224 4.85e-10 ***
## CCAACantabria 3.971e+00 4.459e-01 8.906 < 2e-16 ***
## CCAACastillaLeón 3.530e+00 1.449e-01 24.366 < 2e-16 ***
## CCAACastillaMancha 2.239e+00 1.546e-01 14.484 < 2e-16 ***
## CCAACataluña -5.292e+00 7.216e-01 -7.334 2.23e-13 ***
## CCAACeuta 1.367e+01 3.247e+02 0.042 0.966421
## CCAAComValenciana 3.581e+00 2.138e-01 16.754 < 2e-16 ***
## CCAAExtremadura 8.487e-01 1.661e-01 5.110 3.22e-07 ***
## CCAAGalicia 4.212e+00 3.283e-01 12.829 < 2e-16 ***
## CCAAMadrid 2.991e+00 2.941e-01 10.170 < 2e-16 ***
## CCAAMelilla 1.321e+01 3.247e+02 0.041 0.967545
## CCAAMurcia 2.894e+00 5.513e-01 5.249 1.53e-07 ***
## CCAANavarra 4.416e-01 2.004e-01 2.204 0.027509 *
## CCAAPaísVasco -2.355e+00 4.309e-01 -5.466 4.61e-08 ***
## CCAARioja 3.016e+00 2.717e-01 11.101 < 2e-16 ***
## ActividadPpalHosteleria -9.452e-01 1.982e-01 -4.769 1.85e-06 ***
## ActividadPpalOtro -5.942e-01 2.285e-01 -2.601 0.009305 **
## prop_missings 2.882e+00 6.343e-01 4.543 5.54e-06 ***
## ForeignersPtge 3.706e-02 6.547e-03 5.661 1.51e-08 ***
## DifComAutonPtge -1.563e-02 4.870e-03 -3.209 0.001330 **
## PersonasInmueble -6.074e-03 9.779e-02 -0.062 0.950469
## Age_19_65_pct -3.446e-02 6.603e-03 -5.219 1.80e-07 ***
## Construccion 5.414e-03 2.630e-03 2.058 0.039575 *
## Population 7.701e-06 3.112e-05 0.248 0.804508
## SameComAutonDiffProvPtge 2.248e-02 1.158e-02 1.940 0.052325 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8620.2 on 6495 degrees of freedom
## Residual deviance: 4826.4 on 6467 degrees of freedom
## AIC: 4884.4
##
## Number of Fisher Scoring iterations: 11# Pruebo alguna interaccion sobre el 2
mod4.bin<- c("Age_under19_Ptge", "Age_over65_pct","PersonasInmueble",
"Age_0_4_Ptge", "CodigoProvincia", "Age_19_65_pct",
"PobChange_pct", "UnemployMore40_Ptge", "CCAA", "ActividadPpal",
"CCAA:ActividadPpal")
frml<- as.formula(paste("varObjBin", paste(mod4.bin, collapse="+"), sep="~"))
modelo4<-glm(frml,data=data_train,family=binomial)
summary(modelo4)##
## Call:
## glm(formula = frml, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6423 -0.1401 0.3557 0.5639 3.5359
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) 4.052e+00 1.934e+00 2.095
## Age_under19_Ptge -2.372e-02 2.246e-02 -1.056
## Age_over65_pct -7.871e-03 1.855e-02 -0.424
## PersonasInmueble 1.200e-02 1.083e-01 0.111
## Age_0_4_Ptge 6.237e-02 3.606e-02 1.730
## CodigoProvincia -1.501e-04 3.204e-03 -0.047
## Age_19_65_pct -4.279e-02 1.932e-02 -2.215
## PobChange_pct 3.818e-03 4.269e-03 0.894
## UnemployMore40_Ptge -2.212e-03 1.598e-03 -1.384
## CCAAAragón 1.275e+00 1.115e+00 1.144
## CCAAAsturias -1.685e+00 1.282e+00 -1.314
## CCAABaleares -9.586e-01 5.347e-01 -1.793
## CCAACanarias -3.118e-01 8.224e-01 -0.379
## CCAACantabria 5.878e-01 1.159e+00 0.507
## CCAACastillaLeón 1.659e+01 6.328e+02 0.026
## CCAACastillaMancha 2.078e+00 1.090e+00 1.907
## CCAACataluña -1.858e+01 2.980e+02 -0.062
## CCAACeuta 1.887e+01 3.956e+03 0.005
## CCAAComValenciana 1.655e+01 4.755e+02 0.035
## CCAAExtremadura 1.663e+01 3.956e+03 0.004
## CCAAGalicia 1.145e+00 8.362e-01 1.369
## CCAAMadrid 1.727e+00 6.374e-01 2.709
## CCAAMelilla 1.874e+01 3.956e+03 0.005
## CCAAMurcia 1.655e+01 2.797e+03 0.006
## CCAANavarra -1.469e+00 6.491e-01 -2.263
## CCAAPaísVasco -1.860e+01 9.320e+02 -0.020
## CCAARioja 1.652e+01 1.978e+03 0.008
## ActividadPpalHosteleria -2.379e+00 3.935e-01 -6.046
## ActividadPpalOtro -1.948e+00 4.198e-01 -4.640
## CCAAAragón:ActividadPpalHosteleria 7.644e-01 1.147e+00 0.667
## CCAAAsturias:ActividadPpalHosteleria 2.929e+00 1.317e+00 2.224
## CCAABaleares:ActividadPpalHosteleria 2.848e+00 7.042e-01 4.044
## CCAACanarias:ActividadPpalHosteleria 2.487e+00 8.740e-01 2.846
## CCAACantabria:ActividadPpalHosteleria 3.471e+00 1.308e+00 2.654
## CCAACastillaLeón:ActividadPpalHosteleria -1.361e+01 6.328e+02 -0.022
## CCAACastillaMancha:ActividadPpalHosteleria 6.897e-01 1.111e+00 0.621
## CCAACataluña:ActividadPpalHosteleria 1.466e+01 2.980e+02 0.049
## CCAACeuta:ActividadPpalHosteleria NA NA NA
## CCAAComValenciana:ActividadPpalHosteleria -1.240e+01 4.755e+02 -0.026
## CCAAExtremadura:ActividadPpalHosteleria -1.618e+01 3.956e+03 -0.004
## CCAAGalicia:ActividadPpalHosteleria 3.474e+00 9.268e-01 3.748
## CCAAMadrid:ActividadPpalHosteleria 1.393e+00 7.683e-01 1.813
## CCAAMelilla:ActividadPpalHosteleria NA NA NA
## CCAAMurcia:ActividadPpalHosteleria -1.321e+01 2.797e+03 -0.005
## CCAANavarra:ActividadPpalHosteleria 1.053e+00 7.636e-01 1.379
## CCAAPaísVasco:ActividadPpalHosteleria 1.524e+01 9.320e+02 0.016
## CCAARioja:ActividadPpalHosteleria -1.225e+01 1.978e+03 -0.006
## CCAAAragón:ActividadPpalOtro 2.898e-01 1.132e+00 0.256
## CCAAAsturias:ActividadPpalOtro 1.767e+00 1.405e+00 1.258
## CCAABaleares:ActividadPpalOtro 3.342e+00 1.254e+00 2.666
## CCAACanarias:ActividadPpalOtro 1.899e+01 3.956e+03 0.005
## CCAACantabria:ActividadPpalOtro 2.695e+00 1.385e+00 1.946
## CCAACastillaLeón:ActividadPpalOtro -1.337e+01 6.328e+02 -0.021
## CCAACastillaMancha:ActividadPpalOtro -5.475e-01 1.109e+00 -0.494
## CCAACataluña:ActividadPpalOtro 1.352e+01 2.980e+02 0.045
## CCAACeuta:ActividadPpalOtro NA NA NA
## CCAAComValenciana:ActividadPpalOtro -1.373e+01 4.755e+02 -0.029
## CCAAExtremadura:ActividadPpalOtro -1.617e+01 3.956e+03 -0.004
## CCAAGalicia:ActividadPpalOtro 1.740e+00 1.134e+00 1.535
## CCAAMadrid:ActividadPpalOtro 9.615e-01 8.004e-01 1.201
## CCAAMelilla:ActividadPpalOtro NA NA NA
## CCAAMurcia:ActividadPpalOtro 2.142e+00 3.953e+03 0.001
## CCAANavarra:ActividadPpalOtro 1.599e+00 6.937e-01 2.304
## CCAAPaísVasco:ActividadPpalOtro 1.675e+01 9.320e+02 0.018
## CCAARioja:ActividadPpalOtro -1.415e+01 1.978e+03 -0.007
## Pr(>|z|)
## (Intercept) 0.036196 *
## Age_under19_Ptge 0.290921
## Age_over65_pct 0.671383
## PersonasInmueble 0.911787
## Age_0_4_Ptge 0.083711 .
## CodigoProvincia 0.962628
## Age_19_65_pct 0.026748 *
## PobChange_pct 0.371209
## UnemployMore40_Ptge 0.166297
## CCAAAragón 0.252604
## CCAAAsturias 0.188804
## CCAABaleares 0.073041 .
## CCAACanarias 0.704561
## CCAACantabria 0.611947
## CCAACastillaLeón 0.979088
## CCAACastillaMancha 0.056506 .
## CCAACataluña 0.950281
## CCAACeuta 0.996195
## CCAAComValenciana 0.972232
## CCAAExtremadura 0.996646
## CCAAGalicia 0.171110
## CCAAMadrid 0.006742 **
## CCAAMelilla 0.996220
## CCAAMurcia 0.995280
## CCAANavarra 0.023608 *
## CCAAPaísVasco 0.984078
## CCAARioja 0.993336
## ActividadPpalHosteleria 1.48e-09 ***
## ActividadPpalOtro 3.48e-06 ***
## CCAAAragón:ActividadPpalHosteleria 0.505079
## CCAAAsturias:ActividadPpalHosteleria 0.026178 *
## CCAABaleares:ActividadPpalHosteleria 5.26e-05 ***
## CCAACanarias:ActividadPpalHosteleria 0.004434 **
## CCAACantabria:ActividadPpalHosteleria 0.007964 **
## CCAACastillaLeón:ActividadPpalHosteleria 0.982838
## CCAACastillaMancha:ActividadPpalHosteleria 0.534597
## CCAACataluña:ActividadPpalHosteleria 0.960781
## CCAACeuta:ActividadPpalHosteleria NA
## CCAAComValenciana:ActividadPpalHosteleria 0.979192
## CCAAExtremadura:ActividadPpalHosteleria 0.996737
## CCAAGalicia:ActividadPpalHosteleria 0.000178 ***
## CCAAMadrid:ActividadPpalHosteleria 0.069869 .
## CCAAMelilla:ActividadPpalHosteleria NA
## CCAAMurcia:ActividadPpalHosteleria 0.996232
## CCAANavarra:ActividadPpalHosteleria 0.167771
## CCAAPaísVasco:ActividadPpalHosteleria 0.986952
## CCAARioja:ActividadPpalHosteleria 0.995057
## CCAAAragón:ActividadPpalOtro 0.797978
## CCAAAsturias:ActividadPpalOtro 0.208353
## CCAABaleares:ActividadPpalOtro 0.007678 **
## CCAACanarias:ActividadPpalOtro 0.996169
## CCAACantabria:ActividadPpalOtro 0.051603 .
## CCAACastillaLeón:ActividadPpalOtro 0.983140
## CCAACastillaMancha:ActividadPpalOtro 0.621537
## CCAACataluña:ActividadPpalOtro 0.963814
## CCAACeuta:ActividadPpalOtro NA
## CCAAComValenciana:ActividadPpalOtro 0.976968
## CCAAExtremadura:ActividadPpalOtro 0.996740
## CCAAGalicia:ActividadPpalOtro 0.124847
## CCAAMadrid:ActividadPpalOtro 0.229649
## CCAAMelilla:ActividadPpalOtro NA
## CCAAMurcia:ActividadPpalOtro 0.999568
## CCAANavarra:ActividadPpalOtro 0.021197 *
## CCAAPaísVasco:ActividadPpalOtro 0.985658
## CCAARioja:ActividadPpalOtro 0.994293
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8620.2 on 6495 degrees of freedom
## Residual deviance: 4788.5 on 6435 degrees of freedom
## AIC: 4910.5
##
## Number of Fisher Scoring iterations: 16# Pruebo uno con las variables más importantes del 3
impVariablesLog(modelo3,"varObjBin")
mod5.bin<- c("CCAA", "ActividadPpal", "prop_missings", "ForeignersPtge",
"DifComAutonPtge", "Age_19_65_pct")
frml<- as.formula(paste("varObjBin", paste(mod5.bin, collapse="+"), sep="~"))
modelo5<-glm(frml,data=data_train,family=binomial)
summary(modelo5)##
## Call:
## glm(formula = frml, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8812 -0.2221 0.3649 0.5249 3.6615
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.390873 0.443406 3.137 0.001708 **
## CCAAAragón 1.923127 0.149485 12.865 < 2e-16 ***
## CCAAAsturias 0.944902 0.268311 3.522 0.000429 ***
## CCAABaleares 0.816499 0.315013 2.592 0.009543 **
## CCAACanarias 1.670786 0.280761 5.951 2.67e-09 ***
## CCAACantabria 3.899332 0.441709 8.828 < 2e-16 ***
## CCAACastillaLeón 3.519518 0.142444 24.708 < 2e-16 ***
## CCAACastillaMancha 2.188758 0.149746 14.616 < 2e-16 ***
## CCAACataluña -5.219328 0.720768 -7.241 4.44e-13 ***
## CCAACeuta 13.493874 324.743736 0.042 0.966856
## CCAAComValenciana 3.539079 0.211308 16.748 < 2e-16 ***
## CCAAExtremadura 0.766186 0.160548 4.772 1.82e-06 ***
## CCAAGalicia 4.233436 0.327453 12.928 < 2e-16 ***
## CCAAMadrid 2.904824 0.287411 10.107 < 2e-16 ***
## CCAAMelilla 12.983947 324.743742 0.040 0.968107
## CCAAMurcia 2.795576 0.548465 5.097 3.45e-07 ***
## CCAANavarra 0.341248 0.188697 1.808 0.070538 .
## CCAAPaísVasco -2.362357 0.429324 -5.503 3.74e-08 ***
## CCAARioja 2.913225 0.264077 11.032 < 2e-16 ***
## ActividadPpalHosteleria -0.975342 0.197862 -4.929 8.25e-07 ***
## ActividadPpalOtro -0.788825 0.206962 -3.811 0.000138 ***
## prop_missings 2.840473 0.601012 4.726 2.29e-06 ***
## ForeignersPtge 0.037080 0.006531 5.678 1.36e-08 ***
## DifComAutonPtge -0.015741 0.004737 -3.323 0.000890 ***
## Age_19_65_pct -0.032207 0.006306 -5.107 3.27e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8620.2 on 6495 degrees of freedom
## Residual deviance: 4836.0 on 6471 degrees of freedom
## AIC: 4886
##
## Number of Fisher Scoring iterations: 11#Validacion cruzada repetida para elegir entre todos
auxVarObj<-todo$varObjBin
todo$varObjBin<-make.names(todo$varObjBin) #formateo la variable objetivo para que funcione el codigo
total<-c()
modelos<-sapply(list(modeloInicial,modelo2,modelo3,modelo4,modelo5),formula)
for (i in 1:length(modelos)){
set.seed(1712)
vcr<-train(as.formula(modelos[[i]]), data = todo,
method = "glm", family="binomial",metric = "ROC",
trControl = trainControl(method="repeatedcv", number=5, repeats=20,
summaryFunction=twoClassSummary,
classProbs=TRUE,returnResamp="all")
)
total<-rbind(total,data.frame(roc=vcr$resample[,1],modelo=rep(paste("Modelo",i),
nrow(vcr$resample))))
}Vemos el valor del área bajo la curva ROC
| mean_ROC | std_ROC | variables | test_pseudoR2 | |
|---|---|---|---|---|
| Modelo 4 | 0,895 | 0,008 | 65 | 0,426 |
| Modelo 2 | 0,892 | 0,008 | 29 | 0,445 |
| Modelo 5 | 0,897 | 0,007 | 25 | 0,454 |
| Modelo 3 | 0,897 | 0,007 | 29 | 0,455 |
| Modelo 1 | 0,895 | 0,007 | 55 | 0,459 |
Lo primero que debo hacer es celebrar los resultados ya que tanto el \(pseudoR^2\) como el área bajo la curva ROC dan valores considerablemente altos, lo que es una muy buena señal. Claramente los modelos que compiten pos ser el mejor son el modelo1 generado con todas las variables, el modelo3 generado con las variables que presentaban mayor significancia en el modelo1 y el modelo5 generado reduciendo aún más las variables, escogiendo las más importantes del modelo3. Aunque el tanto el modelo1 como el modelo3 presentan un \(pseudoR^2\) algo superior (muy ligeramente), el modelo5 tiene mayor área bajo la curva roc y el menor número de variables, así que va a ser nuestro elegido.
3.2.2 Determinación del punto de corte
Para elegir el mejor punto de corte lo primero es mirar la distribución de las probabilidades sobre el set de datos test. 
Ya vemos una región que parece interesante, podemos ver algunos valores de precisión, sensibilidad y especificidad de el modelo tomando esos valores:
- Si elegimos 0.55: 0.8391867, 0.8998016, 0.7398374, 0.8500469, 0.8183453
- Si elegimos 0.65: 0.8250154, 0.8541667, 0.7772358, 0.8627255, 0.7648
Y viendo que los resultados son bastante buenos, podemos hacer ir probando los valores dentro de una región (un poco más conservadora) en la que creemos que se encuentra el punto óptimo.
## Mejor punto de corte ----
## generamos una rejilla de puntos de corte
posiblesCortes<-seq(0.25,0.75,0.005)
rejilla<-data.frame(t(rbind(posiblesCortes,sapply(posiblesCortes,function(x) sensEspCorte(modelo5,data_test,"varObjBin",x,"1")))))
rejilla$Youden<-rejilla$Sensitivity+rejilla$Specificity-1
rejilla$posiblesCortes[which.max(rejilla$Youden)]## [1] 0.61rejilla$posiblesCortes[which.max(rejilla$Accuracy)]## [1] 0.43rejilla$posiblesCortes[which.max(rejilla$Youden + rejilla$Accuracy)]## [1] 0.435| Accuracy | Sensitivity | Specificity | Pos Pred Value | Neg Pred Value | |
|---|---|---|---|---|---|
| 0.43 | 0,845 | 0,923 | 0,719 | 0,843 | 0,850 |
| 0.435 | 0,845 | 0,922 | 0,720 | 0,844 | 0,849 |
| 0.61 | 0,837 | 0,882 | 0,763 | 0,859 | 0,798 |
Comparando los resultados para cada punto de corte podemos ver que el valor 0.61 baja demasiado la sensibilidad mientras que entre el 0.43 y el 0.435 las diferencias son poco sensibles. Como diría Aristóteles la virtud está en el centro así que creo que la mejor opción es 0.435, que además obedece al valor donde se maximiza la suma de la precisión con el coeficiente índice de Youden que es un criterio que me he inventado yo (la gráfica ha sido una fuerte sugestión).
Ahora con este valor vamos a ver los coeficientes y variables que participan del modelo ganador. Por supuesto aquí no nos vamos a llevar ninguna sorpresa. La Comunidad Autónoma sigue predominando con excesiva fuerza y el resto parecen prácticamente de apoyo.
Hubiera sido interesante hacer un modelo solo con la Comunidad Autónoma y usarlo como floor, posiblemente fuera un modelo bastante competitivo (y no es otra cosa que ese señor en un bar gritando: Los andaluces son unos rojos, los valencianos son unos fachas, …’ (a ver quién va ahora a decirle que no)).
# Vemos las variables m?s importantes del modelo ganador
par(mar=c(4.1, 12, 4.1, 2.1))
impVariablesLog(modelo5,"varObjBin") 
# Vemos los coeficientes del modelo ganador
coef(modelo5)## (Intercept) CCAAAragón CCAAAsturias
## 1.39087329 1.92312684 0.94490232
## CCAABaleares CCAACanarias CCAACantabria
## 0.81649928 1.67078613 3.89933178
## CCAACastillaLeón CCAACastillaMancha CCAACataluña
## 3.51951801 2.18875785 -5.21932772
## CCAACeuta CCAAComValenciana CCAAExtremadura
## 13.49387404 3.53907865 0.76618650
## CCAAGalicia CCAAMadrid CCAAMelilla
## 4.23343569 2.90482387 12.98394676
## CCAAMurcia CCAANavarra CCAAPaísVasco
## 2.79557607 0.34124772 -2.36235739
## CCAARioja ActividadPpalHosteleria ActividadPpalOtro
## 2.91322493 -0.97534186 -0.78882495
## prop_missings ForeignersPtge DifComAutonPtge
## 2.84047291 0.03708012 -0.01574055
## Age_19_65_pct
## -0.03220663# Evaluamos la estabilidad del modelo a partir de las diferencias en train y test:
train.test<- as.data.frame(rbind(sensEspCorte(modelo5,data_train,"varObjBin",0.495,"1"),
sensEspCorte(modelo5,data_test,"varObjBin",0.495,"1")),
row.names=c("train", "test"))
train.test$auc<- c(roc(data_train$varObjBin, predict(modelo5,data_train,type = "response"),
direction="<")$auc,
roc(data_test$varObjBin, predict(modelo5,data_test,type = "response"),
direction="<")$auc)
train.test$pseudoR2<- c(pseudoR2(modelo5,data_train,"varObjBin"),
pseudoR2(modelo5,data_test,"varObjBin"))
kable(train.test , caption = "Comparación entre train y test",
row.names = TRUE,
digits = 3,
format.args = list(decimal.mark = ",")
)| Accuracy | Sensitivity | Specificity | Pos Pred Value | Neg Pred Value | auc | pseudoR2 | |
|---|---|---|---|---|---|---|---|
| train | 0,841 | 0,915 | 0,719 | 0,842 | 0,837 | 0,898 | 0,439 |
| test | 0,842 | 0,910 | 0,732 | 0,848 | 0,832 | 0,903 | 0,454 |
Sorprendentemente la mayoría de valores mejoran para los datos test así que podemos concluir ampliamente satisfechos que nuestra regresión acierta más del 90% de las veces que región gobernará la derecha y al menos un 70% en cual lo hará la izquierda, lo cual es los tiempos que corren no es baladí.
4 Muchas gracias por todo
Aún tengo muchos problemas con R y pierdo mucho tiempo con temas de librerías y errores varios, así que el trabajo no es todo lo bueno que debería. Además me relaje al saber que iba a profundizar en la regresión logística en el módulo de machine learning. Mis más sinceras disculpas por haber leído tan tarde lo del .pdf, y entregar en .html, siempre me salto la parte clave de los enunciados. En definitiva, muchas gracias por su tiempo, ha sido un trabajo muy estimulante.