census_income.Rmd

---
title: "Analizando el dataset Census Income"
author: "Marta Gómez y Braulio Vargas"
date: "3 de junio de 2016"
lang: es
output: 
  pdf_document:
    fig_caption: yes
    includes:
      in_header: structure.tex
    number_sections: yes
    toc: yes
    highlight: pygments
bibliography: bibliografia.bib
csl: ieee-with-url.csl
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	message = FALSE,
	warning = FALSE,
	cache = TRUE
)
```

# Definición del problema a resolver y enfoque elegido

## Definición del problema

La base de datos escogida se denomina _Census Income_, aunque también es conocida como _Adult_ [@Lichman:2013]. Se puede descargar [aquí](http://archive.ics.uci.edu/ml/datasets/Census+Income "aqui"). Esta base de datos tiene 48842 instancias, cada una con 14 atributos. El problema es que también tiene valores perdidos.

El problema consiste en predecir cuando una persona ingresa más de 50.000 dólares al año, basándonos en los datos del censo. Por tanto, estamos ante un problema de _clasificación binaria_.

Los atributos de los que disponemos son:

* ___age___: atributo numérico que expresa la edad de la persona.

* ___workclass___: atributo categórico que expresa el tipo de empleo que tiene la persona. Los valores que puede tomar son:

    * _Private_: persona contratada en una empresa privada.
    
    * _Self-emp-not-inc_: autónomo.
    
    * _Self-emp-inc_: persona que tiene una empresa grande.
    
    * _Federal-gov_: funcionario del gobierno federal.
    
    * _Local-gov_: funcionario del gobierno local.
    
    * _State-gov_: funcionario del gobierto estatal.
    
    * _Without-pay_: persona en paro.
    
    * _Never-worked_: persona que nunca ha trabajado.
    
* ___fnlwgt (final weight)___: variable numérica que representa un peso, cada peso corresponde a una característica socio-económica de la población, por tanto, personas con características demográficas parecidas deben tener un peso parecido. [@fnlwgt]

* ___education___: variable categórica que representa el nivel de estudios de la persona. Los valores que puede tomar son: _Bachelors_, _Some-college_, _11th_, _HS-grad_, _Prof-school_, _Assoc-acdm_, _Assoc-voc_, _9th_, _7th-8th_, _12th_, _Masters_, _1st-4th_, _10th_, _Doctorate_, _5th-6th_, _Preschool_.

* ___education-num___: variable númerica que representa a la anterior.

* ___marital-status___: variable categórica que expresa el estado civil de la persona. Los valores que puede tomar son: 

    * _Married-civ-spouse_: la persona está casada con un civil.
    
    * _Divorced_: la persona está divorciada.
    
    * _Never-married_: la persona nunca ha estado casada.
    
    * _Separated_: la persona está separada.
    
    * _Widowed_: la persona es viuda.
    
    * _Married-spouse-absent_: la persona aparece como casada en el registro, pero no se encuentra a ninguna pareja [@maritalstatus].
    
    * _Married-AF-spouse_: la persona está casada con alguien de las fuerzas armadas.
    
* ___occupation___: variable categórica que describe el tipo de profesión que tiene la persona. Puede tomar los valores _Tech-support_, _Craft-repair_, _Other-service_, _Sales_, _Exec-managerial_, _Prof-specialty_, _Handlers-cleaners_, _Machine-op-inspct_, _Adm-clerical_, _Farming-fishing_, _Transport-moving_, _Priv-house-serv_, _Protective-serv_, _Armed-Forces_. 

* ___relationship___: variable que contiene el valor que puede tomar la relación de una persona con respecto a otra dentro de una familia. Contiene solo un valor por instancia del dato. Estos valores pueden ser _Wife_ (Esposa), _Own-child_ (hijo propio), _Husband_ (marido), _Not-in-family_ (sin familia), _Other-relative_ (otro tipo de familiar) y _Unmarried_ (soltero).

* ___race___: valor que describe la raza de la persona. Puede ser _White_ (Blanco), _Asian-Pac-Islander_ (asíatico o de las islas del pacífico), _Amer-Indian-Eskimo_ (indio americano o esquimal), _Other_ (otro) y _Black_ (negro).

* ___sex___: sexo de la persona, _Female_ (mujer) o _Male_ (hombre).

* ___capital-gain___: registro de la ganancia de capital de la persona.

* ___capital-loss___: registro de la pérdida de capital de la persona.

* ___hours-per-week___: corresponde a las horas de trabajo a la semana.

* ___country___: país de pertenencia. En este caso, hay una mayor variedad de valores, siendo los siguientes: _United-States_, _Cambodia_, _England_, _Puerto-Rico_, _Canada_, _Germany_, _Outlying-US(Guam-USVI-etc)_, _India_, _Japan_, _Greece_, _South_, _China_, _Cuba_, _Iran_, _Honduras_, _Philippines_, _Italy_, _Poland_, _Jamaica_, _Vietnam_, _Mexico_, _Portugal_, _Ireland_, _France_, _Dominican-Republic_, _Laos_, _Ecuador_, _Taiwan_, _Haiti_, _Columbia_, _Hungary_, _Guatemala_, _Nicaragua_, _Scotland_, _Thailand_, _Yugoslavia_, _El-Salvador_, _Trinadad&Tobago_, _Peru_, _Hong_, _Holand-Netherlands_,  o lo que es lo mismo, Estados-Unidos, Colombia, Inglaterra, Puerto-Rico, Canada, Alemania, Estados supervisados de los Estados Unidos, India, Japón, Grecia, Sudáfrica, China, Cuba, Irán, Honduras, Filipinas, Italia, Polonia, Jamaica, Vietnam, México, Portugal, Irlanda, Francia, República Dominicana, Laos, Ecuador, Taiwn, Haití, Colombia, Hungría, Guatemala, Nicaragua, Escocia, Tailandia, Yugoslavia, El Salvador y Trinidad-Tobago.

* ___income___: corresponde al valor que queremos predecir, puede tomar los valores "<=50K" o ">50K"

## Enfoque elegido

El enfoque que se usará en este trabajo consiste en realizar una serie de modelos _**M'**_ usando distintos algoritmos. Cada uno de ellos tendrá un error de train y test diferentes, una complejidad distinta y comportamientos distintos a la hora de generalizar. Tras analizar cada uno de los distintos modelos, ver su complejidad, su idoneidad para el conjunto de datos, el error de test... se escogerá un modelo _**M**_ que será el "mejor" del conjunto de modelos desarrollados. A lo largo del trabajo se desarrollarán distintos modelos usando los siguientes algoritmos:

  - **Modelo lineal**: se usará como modelo lineal la *regresión logística*.  
  - **Modelo Random Forest**
  - **Modelos de Base Radial**: en este caso, usaremos dos modelos de base radial diferentes:
    - **Support Vector Machine**
    - **KNN**
  - **Red Neuornal**

# Codificación de los datos de entrada para hacerlos útiles a los algoritmos

## Valores perdidos

Respecto a los valores perdidos, podríamos tomar tres enfoques diferentes: sustituirlos por la media del resto de valores para ese dato predictor o predecirlo con un modelo, eliminar todas las filas que contengan algún valor perdido o no usar las variables que contengan valores perdidos. Para tomar esta decisión, primero hemos de estudiar el porcentaje de valores perdidos que presentan los datos y dónde se pierden estos datos.

Para ello, en primer lugar leemos los datos. 

```{r}
leer_datos <- function(fichero = "./Data/adult.data") {
    adult.train <- read.csv(fichero, header=FALSE, col.names = c("age","workclass", 
            "fnlwgt","education","education-num","marital-status","occupation","relationship",
            "race","sex","capital-gain","capital-loss","hours-per-week","country","income"), 
            na.strings = c(" ?", "?", ""), stringsAsFactors = F)
}

adult.train <- leer_datos()
adult.test <- leer_datos(fichero = "./Data/adult.test")
```

Para evitar problemas con las funciones de $R$, forzaremos a que todos los datos categóricos, tenga el mismo $factor$. Para ello, realizaremos lo siguiente:

```{r}
# Añadimos una columna a los datos para indicar 
# cuales pertenecen a datos de train y cuales a
# datos de test

adult.test$trainTest = rep(1,nrow(adult.test))
adult.train$trainTest = rep(0,nrow(adult.train))

# Reconstruimos el conjunto de datos al completo,
# uniendo los datos de train y test

fullSet <- rbind(adult.test,adult.train)

# Cada una de las variables categóricas, pasarán 
# de ser cadenas de caracteres a tener un valor 
# numérico o un factor, con lo que 
# tanto datos de train como datos de test, 
# obtendrán el mismo factor

fullSet$workclass = as.factor(fullSet$workclass)
fullSet$country = as.factor(fullSet$country)
fullSet$education = as.factor(fullSet$education)
fullSet$marital.status = as.factor(fullSet$marital.status)
fullSet$sex = as.factor(fullSet$sex)
fullSet$relationship = as.factor(fullSet$relationship)
fullSet$occupation = as.factor(fullSet$occupation)
fullSet$income = as.factor(fullSet$income)
fullSet$race = as.factor(fullSet$race)

# Reconstruimos los datos de train y test originales

adult.train = data.frame(fullSet[fullSet$trainTest == 0,])
adult.test = data.frame(fullSet[fullSet$trainTest == 1,])

# Y eliminamos la columna auxiliar

adult.test$trainTest = NULL
adult.train$trainTest = NULL
```


Una vez leídos, calculamos el número de variables perdidas en cada columna:

```{r}
apply(X=adult.train, MARGIN=2, FUN=function(columna) length(is.na(columna)[is.na(columna)==T]))
```

Sólo tienen valores perdidos los atributos _workclass_, _occupation_ y _country_.

También, vamos a comprobar el número de instancias que contienen algún dato perdido. Como ya hemos visto que los datos perdidos sólo se encuentran en las columnas *workclass*, *occupation* y *country*, nos fijaremos sólo en esas. Para ello, vamos a realizar lo siguiente:

```{r}
getRowsNA <- function(datos = adult.train) {
    aux = is.na(datos)*1
    rowsMissingValues = apply(X=aux, MARGIN=1,
             FUN = function(fila) sum(fila))
}

rowsMissingValues.train = getRowsNA()
length(rowsMissingValues.train[rowsMissingValues.train > 0])
```

En total, sólo hay 2399 filas con valores perdidos. Al tener en total 32561 datos de entrenamiento, perder 2399 no supone una gran diferencia, sólo perdemos un $7,37$% de los datos. Por tanto, lo más sencillo es eliminar las filas que contengan algún dato perdido. Este mismo planteamiento es usado con los datos de test, donde tenemos 16282 filas en total, de las cuales 1222 presentan datos perdidos, un $7,505$% de los datos.

```{r}
adult.train.clean = adult.train[rowsMissingValues.train == 0,]
adult.test.clean = adult.test[getRowsNA(datos = adult.test) == 0,]
fullSet.clean = data.frame(fullSet[getRowsNA(datos = fullSet) == 0,])
```

## Asignación de valores numéricos a los valores categóricos

Al cargar el fichero de datos hemos usado la opción `stringsAsFactors`. Esta función sirve para convertir todas las cadenas de carácteres en factores. Los factores, tal y como se explica en [@factores], asignan a cada posible valor que puede tomar una variable categórica un número. Los factores representan una forma de almacenamiento muy eficiente, ya que las cadenas de carácteres correspondientes sólo se almacenan una vez y los valores se guardan como enteros.

Por ejemplo, los números asociados a la variable `workclass` corresponden con el orden mostrado por la función `levels`. Así, al primer elemento se le asigna el número uno, al siguiente el dos, etc.

```{r}
levels(adult.train.clean$workclass)
```

Para ver el vector numérico almacenado por `R` en vez de el vector con las cadenas de caracteres usamos la función `c`:

```{r}
head(c(adult.train.clean$workclass))
head(adult.train.clean$workclass)
```

Como vemos, los números asignados a los primeros elementos coinciden con el orden mostrado anteriormente.

## Eliminando la variable `education.num`

Aunque cada problema tenga su propio método para escoger y normalizar variables, hemos considerado que la variable `education.num` no nos va a servir, ya que es únicamente una asignación numérica de los valores categóricos de la variable `education`. Al haber almacenado la variable `education` como factor, ya tiene una asignación numérica hecha, y por tanto, eliminamos la variable `education.num`:

```{r}
adult.train.clean[,"education.num"] = NULL
adult.test.clean[,"education.num"] = NULL
adult.train[,"education.num"] = NULL
adult.test[,"education.num"] = NULL
```

# Ajustando un modelo lineal

## Idoneidad del modelo lineal
Antes de pasar a modelos más complejos, vamos a comprobar cómo se comporta un modelo lineal. Esto es así, ya que la complejidad de un modelo lineal es bastante menor que a la de los otros modelos que se usarán a continuación, y si es capaz de ajustar bastante mejor los datos que el resto de modelos, se escogería este.

Para ello, nos vamos a basar en un modelo que utiliza la regresión logística para este problema. Este modelo, nos devuelve la probabilidad de pertenecer a una clase u otra. Usando el concepto de \textit{maximum likelihood} o máxima probabilidad, asigna a los datos la clase a la que es más probable que pertenezca, de acuerdo a la probabilidad que calcula el modelo. 

Estas probabilidades se obtienen a partir de la siguiente expresión:
$$ P(y|x) = \theta(yw^Tx)$$
donde $\theta$ es la función logística, $y$ es la etiqueta del modelo, $w$ es nuestro vector de pesos asociado, y $x$ los datos predictores. Esta función es la que usaremos para crear un modelo de estimación del error, y trataremos de minimizar ese error, que lo calcularemos de la siguiente forma [@lfd]:
$$E_{in} = \frac{1}{N}\sum_{n=1}^N\ln\left(\frac{1}{\theta(y_nw^Tx_n)} \right) = \frac{1}{N}\sum_{n=1}^N\ln\left(1 + e^{-y_nw^Tx_n}\right)$$

## Dimensión del modelo lineal

La dimensión del modelo lineal la mediremos en función de la dimensión de **Vapnik–Chervonenkis** o dimensión *VC*, denotada como $d_{VC}$ que expresa la complejidad del modelo, y el máximo valor de $N$, donde $N$ es el número de datos disponibles, que el modelo puede separar. En el caso de los modelos lineales, como es el caso, $d_{VC} = d + 1$, donde $d$ es el número de variables predictoras que tenemos.

Con esto vemos, que la dimensión $d_{VC}$ es basntante competitiva y puede comportarse bien en el problema.

## Parámetros del modelo

En el modelo de regresión logística, el principal parámetro a configurar es el vector de pesos $w$ que utiliza la regresión. Este vector de pesos se puede inicializar a un valor, fijo o un valor aleatorio, que de algún modo marca el punto de inicio del algoritmo de regresión. Este parámetro se irá actualizando en la regresión logística conforme se vayan comprobando nuevos puntos.

```{r}
predictorGLM <- function(model, pintar = T){
  ypred = predict(model, adult.test.clean, type="response")
  ypred[ypred <= 0.5] = ">50K"
  ypred[ypred > 0.5] = "<=50K"

  if(pintar){
    print("Matriz de confusión datos de test")
    print(table(predict=ypred, truth=(adult.test.clean$income)))
  }
  cat("Eout = ",mean((ypred != adult.test.clean$income)*1))
  ypred
}
trainingIndex=which(fullSet.clean$trainTest==0)
fullSet.clean$trainTest = NULL
fullSet.clean$education.num = NULL
set.seed(1)
glmModel = glm(income ~ ., data = fullSet.clean, 
    subset = trainingIndex, family = binomial(logit))
glmPred = predictorGLM(glmModel, pintar=T)
```

Como podemos ver, el error que genera la regresión lineal, tomando todas las variables como predictoras, genera un error de aproximadamente el 24%. Como vemos, el error generado por la regresión logística es bastante alto y por ello, vamos a comprobar otros modelos. El principal problema de la regresión logística, es que se ve muy afectado por las clases desbalanceadas como es el caso, donde la clase perteneciente al conjunto `>50K` tiene muchas menos instancias que la clase `<=50K` por lo que el modelo puede ajustar muy bien la clase mayoritaria, mientras que la clase minoritaria prácticamente la obvia. 

```{r}
cat("Error obtenido en la clase <=50K: ", 
    0/length(adult.test.clean$income[adult.test.clean$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    370000/length(adult.test.clean$income[adult.test.clean$income==">50K"]))
```

Como vemos, el modelo no ha fallado ni un solo caso de la clase `<=50K`, pero ha fallado todos los casos `>50K`. Esto se debe a que todas las instancias de test han sido clasificadas como `<=50K`.

```{r}
plot(glmModel, which=c(1))
```

En esta gráfica vemos cómo el error del modelo sigue una tendencia no lineal, esto nos indica que debemos ajustar un modelo más complejo para poder ajustar bien los datos.

# Ajustando un modelo _Random Forest_

## Idoneidad del modelo *Random Forest* para los datos del problema

Un modelo usando un árbol como base, es capaz de clasificar los datos con cierta precisión y tiene la ventaja de que por sí solo, es capaz de explicar bastante bien la clasificación que ha hecho. Pero, utilizar un sólo árbol no ofrece la potencia necesaria para ajustarse lo suficientemente bien a los datos, además de que sufre de tener una alta variabilidad, ya que si particionamos los datos aleatoriamente en subconjuntos, podemos obtener resultados muy diferentes si aprendemos el modelo con un subconjunto u otro, y por ello usaremos el modelo *Random Forest*, que aunque pierda las facilidades de comprensión que tiene un árbol, nos ofrece un poder de ajuste muchísimo mayor, y un poder de generalización mayor.

El porqué de usar *Random Forest* frente a un modelo *Bagging* es que con este último, tenemos que aprender el modelo usando las $p$ variables predictoras, mientras que con *Random Forest* podemos usar un subconjunto de tamaño $m$. *Bagging*, al hacer uso de los $p$ predictores, en el caso de que haya una variable predictora muy fuerte en los datos, todos los árboles que se acaben generando, tendrán esta variable predictora en los primeros nodos del árbol, mientras que al usar varios subconjuntos de los $p$ predictores de tamaño $m$ escogidos de forma aleatoria en *Random Forest*, seremos capaces de detectar otras relaciones más débiles entre las variables predictoras pero que también son capaces de predecir los datos con exactitud.

Respecto a los datos que tenemos (30162 datos de entrenamiento y 15060 datos de test)\footnote{Eliminando las instancias con valores perdidos.} tenemos suficientes datos como para que *Random Forest* trabaje de forma adecuada, y disponemos de 14 datos predictores, ya que la variable `income` es aquella que queremos predecir.

## Parámetros del modelo

El parámetro que utiliza el modelo *Random Forest*, es el parámetro $m$, que indica el número de variables que se escogeran de forma aleatoria para crear un modelo de árbol para ajustar los datos cada vez que el modelo vaya a ajustar un nuevo árbol. Al no ser un modelo paramétrico, no tiene mucho más modificar, ya que los pesos no se pueden establecer desde un principio.

## Hiperparámetros del modelo

En este caso, como hiperparámetro se utiliza el **número de árboles** que generará el modelo para ajustar los datos. En cierto modo, es similar a la regularización, ya que a mayor número de árboles, más sobreajustaremos el modelo, mientras que con un menor número de árboles, menos se sobreajusta el modelo.

## Complejidad del modelo

La dimensión $d_{VC}$ para un modelo de árbol, es infinita, ya que para un conjunto de datos, podemos añadir al árbol tantas variables como queramos para clasificar de forma correcta todos los datos, y si tenemos infinitos datos, podemos añadir infinitas  variables al árbol para separar hasta el último punto del espacio, con lo que $d_{VC} = \infty$. 

Si esto lo aplicamos a *Random Forest* tenemos que la dimensión $d_{VC} = \infty$ para este modelo también, al estar construido sobre árboles.

## Aplicación del modelo *Random Forest*

Para ajustar un modelo *Random Forest*, haremos uso del paquete *randomForest* de *R*. 

En primer lugar, vamos a ajustar un modelo _Random Forest_ con el número de predictores por defecto de `R`: $m = \sqrt{p}$. Siendo $p$ el número de variables predictoras totales. Este número por defecto es el recomendado para problemas de clasificación como éste. En este caso, al tener $p = 13$, el número de predictores a usar será $m = 3$.

```{r}
library(randomForest)

outOfSampleError <- function(model, newdata = adult.test.clean,  printTable = T, ...){
  set.seed(1)
  ypred = predict(object = model, newdata = newdata, ...)
  if(printTable)
    print(table(predict=ypred, truth=newdata$income))
  cat("Eout = ",mean((ypred != newdata$income)*1))
}
set.seed(1)
rf.clean = randomForest(income ~ ., data = adult.train.clean, importance = T)
print(rf.clean)
outOfSampleError(rf.clean)
```

Con $m = 3$ obtenemos un modelo con un $18,18$% de error _Out Of Bag_. Al ajustar el modelo de Random Forest, cada uno de los árboles usa un subconjunto de los datos, y los no usados para calcular dicho árbol se quedan como _Out Of Bag_. Una vez generados los árboles, se predice el valor de cada dato usando los árboles en los que dicho dato ha sido _Out Of Bag_. En el caso de clasificación se hace mediante voto mayoritario. El error de clasificación obtenido para cada dato es válido, porque para predecir el dato se han usado árboles en los que el dato en cuestión no "ha participado" [@islroob].

Respecto a la matriz de confusión generada, vemos que hay muchos más aciertos para la clase `<=50K` que para la clase `>50K` y que la tasa de error en la primera clase es mucho menor que la de la segunda, $0,002$ frente a $0,72$. Esto puede deberse a que el número de instancias para el la clase `<=50K` sea mucho mayor y, por tanto, el modelo tienda a clasificar los datos como `<=50K`:

```{r}
length(adult.train.clean$income[adult.train.clean$income == "<=50K"])
length(adult.train.clean$income[adult.train.clean$income == ">50K"])
```

Como hemos dicho antes, hay una gran diferencia entre la proporción de datos para un clase y otra, esto hace que el modelo tienda a clasificar los datos como `<=50K`. De hecho, si nos fijamos en la matriz de confusión, los fallos de clasificación con `<=50K` son mucho mayores a los de `>50K`: 5446 frente a 54.

Ahora bien, ¿podemos llegar a obtener un error más bajo cambiando éste parámetro? Para ello, usamos la función `tuneRF`, que empieza en el valor por defecto de `R` y, a partir de éste, prueba con otros valores hasta encontrar el óptimo. En cada iteración, incrementa o decrementa el valor de $m$ según un parámetro `stepFactor` y además, si la mejora de error obtenida no es mayor a lo indicado por el parámetro `improve`, se para la búsqueda.

```{r}
set.seed(1)
rf.tuned = tuneRF(x=subset(adult.train.clean, select=-income), 
      y=adult.train.clean$income, doBest=T)
```

Como vemos en la gráfica, con $m=2$ obtenemos el modelo óptimo con menor error, por tanto, en vez de quedarnos con el modelo generado anteriormente con $m=3$, nos quedamos con el mejor obtenido según el error _Out-Of-Bag_, junto con la matriz de confusión para los datos de train, donde vemos que el error para la clase `>50K` es del 72%.

Lo último que podemos preguntarnos es, ¿cómo se comporta el modelo con los datos de test? Para saberlo, generamos la tabla de confusión:

```{r}
outOfSampleError(rf.tuned)
plot(rf.clean$predicted, col = c(3,2), main = "Nº de datos predichos para cada clase por Random Forest")
cat("Error obtenido en la clase <=50K: ", 
    2700/length(adult.test.clean$income[adult.test.clean$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    269900/length(adult.test.clean$income[adult.test.clean$income==">50K"]))
```

El $E_{out}$ obtenido es muy parecido al _Out-Of-Bag_ aunque algo mayor (17,9 en train frente a 18,1 en test), ya que éste es una buen estimador del error fuera de la muestra tal y como hemos compobado.

En esta matriz de confusión se aprecia lo desbalanceadas que están las clases en este problema: en la clase `<=50K` hemos obtenido un error del $0,2$%, mientras que en la clase `>50K` el error ha sido del 72%. 

En los datos de test también se aprecia lo desbalanceadas que están las clases en el problema y esto influye también en los resultados obtenidos:

<!-- De los 15060 datos de test que teníamos en total, el modelo ha acertado 123313 y ha fallado 2747. Al igual que en los datos de entrenamiento, se ve una mayoría de datos con clase `<=50K`. Aunque a la hora de los errores, ha habido más errores clasificando instancias `<=50K`: 34 frente a 2713. -->

```{r}
length(adult.test.clean$income[adult.test.clean$income == "<=50K"])
length(adult.test.clean$income[adult.test.clean$income == ">50K"])
```

Al haber esta diferencia en los datos de entrenamiento y test, el modelo se ve muy afectado.

# Ajustando un modelo de base radial

## Support Vector Machine

### Idoneidad del modelo *Support Vector Machine* para los datos del problema

El uso de los modelos de \textit{Support Vector Machine} o $SVM$ es muy apropiado para el problema, ya que funciona bastante bien con los problemas de clasificación binaria, como es el caso, ya que ofrece una aproximación muy natural al problema, y un buen comportamiento en problemas con grandes dimensiones, gracias al uso del *Kernel*.

Además, los modelos $SVM$, tienen un rendimiento muy bueno frente a conjuntos de datos con clases desbalanceadas, como es nuestro caso, donde la clase con un valor de $income =$ "$<=50K$" representa un subconjunto mucho mayor que la clase con $income =$ "$>50K$"", como podemos ver en [@unbalanced].

### Parámetros del modelo

Los parámetros a ajustar del modelo son:

* *Kernel del SVM*: será el kernel que vamos a usar, es decir, cómo de importantes serán el resto de instancias de los datos para ajustar un dato, en función de la distancia entre los datos. En este caso, se usará un modelo de con un kernel de base radial, que tiene como base la siguiente función: $$w(u,v)=e^{-\gamma\cdot|u-v|^2}$$

* *$\gamma$*: valor que tomará $\gamma$ para la función de base radial. Normalmente, se suele hacer que $\gamma = \frac{1}{d}$, donde $d$ es la dimensión del problema.

* *$\epsilon$*: valor de $\epsilon$ para la función de pérdida de SVM. En este caso valdrá 0.1.

* *Ponderización de la clase*: en caso de que sea necesario, podemos ponderar las clases si queremos para balancear el conjunto de datos, dando más peso a la clase minoritaria que a la mayoritaria. En un principio, las clases tendrán el mismo peso, es decir 1.

### Hiperparámetros del modelo

En este caso, como hiperparámetro del modelo tenemos el parámetro $\lambda$ o *cost*, que controla la regularización del problema. Esta regularización de los datos, se realiza haciendo uso de la formulación de Lagrange.

### Complejidad del modelo

La complejidad de los modelos $SVM$ varía en función del $Kernel$ escogido para la función. En el caso de los modelos de base radial, como es el caso, esta dimensión $d_{VC}$, ya que si nuestro espacio de búsqueda se establece en $\mathbf{R}^d$, dado un punto $x_1$, podemos coger un punto $x_2 \in R^d$ y medir la distancia entre estos dos puntos para ver lo que aporta el modelo. Como esta distancia puede ser infinita, $d_{VC} = \infty$.[@dvcsvm]

### Ajuste del SVM al problema

### Normalización de los datos

Para ver como se comporta *Support Vector Machine* al problema, vamos 

```{r}
library(e1071)
set.seed(1)
svmModel = svm(income ~ ., data = adult.train.clean, kernel = "radial")
outOfSampleError(svmModel)
```

Como vemos, el error fuera de la muestra supone un 15\% de $E_{out}$ aproximadamente con un kernel de base radial, un valor de $\gamma = \frac{1}{14}$, $\epsilon = 0,1$, con el mismo peso para las dos clases, y con el valor *cost = 1*. En la matriz de confusión, vemos como este modelo se comporta mejor que otros modelos frente a la clase minoritaria, reduciendo mucho la tasa de error en esta.

A continuación, vamos a realizar tres nuevos modelos con $SVM$, añadiendo regularziación al modelo. Este valor de regularización tomará los valores de $\lambda = 0,001$, $\lambda = 0,01$, $\lambda = 0,1$ y $\lambda = 0,9$.

```{r}
set.seed(1)
svmModelReg_0001 = svm(income ~ ., data = adult.train.clean, kernel = "radial", cost = 0.001)
outOfSampleError(svmModelReg_0001)
set.seed(1)
svmModelReg_001 = svm(income ~ ., data = adult.train.clean, kernel = "radial", cost = 0.01)
outOfSampleError(svmModelReg_001)
set.seed(1)
svmModelReg_01 = svm(income ~ ., data = adult.train.clean, kernel = "radial", cost = 0.1)
outOfSampleError(svmModelReg_01)
set.seed(1)
svmModelReg_09 = svm(income ~ ., data = adult.train.clean, kernel = "radial", cost = 0.9)
outOfSampleError(svmModelReg_09)
set.seed(1)
svmModelReg_25 = svm(income ~ ., data = adult.train.clean, kernel = "radial", cost = 2.5)
outOfSampleError(svmModelReg_25)
```

Como podemos ver, los resultados cambian dependiendo del valor de $\lambda$ que tome el modelo. Frente al valor de $\lambda$ del primer modelo $SVM$, que era 1, el resto de modelos empeoran al cambiar el valor de la regularización, pasando de una gran deterioro del modelo en valores pequeños de $\lambda$ ($0,001$ y $0,01$) a un deterioro más pequeño con $\lambda = 0,1$. En todos estos casos, estábamos sobreaustando el modelo. Sin embargo, con $\lambda = 0,9$ el modelo es capaz de mejorar un poco al primer modelo y un con un $\lambda$ mayor, podíamos llegar a subajustar el modelo. También podemos ver como si hacemos el valor de $\lambda$ mayor, el modelo no realiza mejoras significativas, y en este caso, incluso empeora un poco.

```{r}
plot(svmModelReg_09$fitted, col = c(3,2), main = "Nº de datos predichos para cada clase por SVM")
cat("Error obtenido en la clase <=50K: ", 
    71400/length(adult.test.clean$income[adult.test.clean$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    151300/length(adult.test.clean$income[adult.test.clean$income==">50K"]))
```

## Vecino más cercano con base radial

### Idoneidad del modelo de _Vecino más cercano_ para los datos del problema

La razón principal para usar el vecino más cercano es que tenemos una gran cantidad de datos y además, tenemos una dimensión $d=15$. Por tanto, tenemos instancias suficientes como para que el modelo del vecino más cercano actúe de forma correcta.

Otro punto a favor para el modelo del _Vecino más cercano_ es que funciona bien con clases desbalanceadas, tal y como se indica en [@unbalanced].

Además, al usar una aproximación con base radial, obtenemos una clasificación mucho mejor dando a cada vecino un peso dependiende de su _distancia de Minkowski_ [@teoriakknn], que es una generalización de la _distancia Euclídea_ y la _distancia de Manhattan_:

$$d(x_i,x_j) = \Bigg(\sum_{s=1}^p | x_{is} - x_{js}|^q \Bigg)^{\frac{1}{q}}$$

Hay varios kernel distintos para actualizar los pesos, como por ejemplo una _gaussiana_ o la inversa de la distancia:

$$Gauss\;kernel =  \frac{1}{\sqrt{2\pi}}\mathrm{exp}\bigg(-\frac{d^2}{2}\bigg) \qquad\ Inversion\;kernel =  \frac{1}{|d|}$$

### Dimensión VC del modelo del Vecino más cercano

Al igual que _Random Forest_, la $d_{VC}$ de éste modelo es $d_{VC} = \infty$. Esto se cumple tanto para la técnica $1-NN$ como para la técnica $k-NN$ y se debe a que el conjunto de hipótesis $\mathcal{H}$ que genera el modelo ajusta con $E_{in} = 0$ cualquier dataset, sin importar el tamaño tal y como se indica en [@lfd].

### Normalización de los datos

En primer lugar, vamos a normalizar las variables numéricas de nuestro modelo: `age`, `fnlwgt`, `capital.gain`, `capital.loss` y `hours.per.week`.

Para ello, he desarrollado una función que normaliza los datos en el intervalo $[0,1]$ aplicando la siguiente fórmula a cada uno de los datos:

$$z_i = \frac{x_i - min(x)}{max(x) - min(x)}$$

donde $x$ representa el vector de datos originales y $z$ el vector de datos normalizados. Este método busca normalizar en base al mínimo (que representaría el 0) y el máximo (que representaría el 1) todos los datos. En `R`, la forma de implementar este método es mediante la función `scale`:

```{r}
normalizar_maxmin <- function(data=adult.train.clean, numbercols = c("age","fnlwgt","capital.gain",
                            "capital.loss","hours.per.week"), data_test = adult.test.clean) {
    # nos quedamos sólo con las columnas numéricas
    cols_numericas = subset(data, select=numbercols)
    colstest_numericas = subset(data_test, select=numbercols)
    # calculamos el máximo y el mínimo de cada columna de los datos de train
    maxs = apply(X=cols_numericas, MARGIN=2, FUN=max)
    mins = apply(X=cols_numericas, MARGIN=2, FUN=min)
    # aplicamos el escalado a los datos de train
    datos_normalizados_numericos = scale(x=cols_numericas, center=mins, scale=maxs)
    # aplicamos los valores de normalización de train sobre los de test
    test_normalizado = as.data.frame(scale(x = colstest_numericas, 
                    center = attr(datos_normalizados_numericos, "scaled:center"),
                    scale = attr(datos_normalizados_numericos, "scaled:scale")))
    # juntamos los valores normalizados con el resto de columnas
    datos_normalizados_numericos = as.data.frame(datos_normalizados_numericos)
    for (c in numbercols) {
        data[c] = datos_normalizados_numericos[c]
        data_test[c] = test_normalizado[c]
    }
    list(data, data_test)
}

norm = normalizar_maxmin()
adult.train.clean.norm = norm[[1]]
adult.test.clean.norm = norm[[2]]
norm = NULL
```

### Aplicación del modelo de _Vecino más cercano con base radial_

Para elegir el $k$ a usar, vamos a entrenar el modelo con _Validación cruzada leave-one-out_. Como `kmax` vamos a poner el doble de la dimensión del dataset, para que la función tenga flexibilidad en cuanto al número de vecinos entre los que escoger. Como kernel, vamos a elegir entre uno gaussiano y otro inversamente proporcional a la distancia, ya que elegir entre alguno más incrementa muchísimo el tiempo de ejecución.

```{r, message=FALSE, warning=FALSE}
library(kknn)
set.seed(1)
best_model_knn <- train.kknn(formula = income ~ ., data = adult.train.clean.norm, 
  kmax = 2*ncol(adult.train), kernel = c("gaussian", "inversion"))
best_model_knn
```

El error obtenido con este modelo ha sido del 16%. A pesar de que este error haya sido con los datos de entrenamiento, es una buena medida debido a que ha sido calculado mediante validación cruzada _leave-one-out_ y la estimación se ha hecho en base a ese dato que se ha dejado fuera y no ha participado en el modelo. Como esto se ha repetido una vez por cada dato que tenemos, obtenemos una buena estimación del error del modelo.

Una vez obtenidos los mejores parámetros, vamos a ajustar el modelo para obtener también un vector con predicciones y poder hacer la matriz de confusión:

```{r}
set.seed(1)
model_knn <- kknn(formula = income ~ ., train = adult.train.clean.norm, 
    test = adult.test.clean.norm, k = best_model_knn$best.parameters$k, 
    kernel = best_model_knn$best.parameters$kernel)
print(table(predict = model_knn$fitted.values, truth = adult.test.clean$income))
cat("Eout = ",mean((model_knn$fitted.values != adult.test.clean$income)*1))
```

El error de test obtenido ha sido del 16%, al igual que el error de validación cruzada. En este caso el modelo ha conseguido corregir las clases desbalanceadas que había, equilibrando el error obtenido en una y en otra:

```{r}
plot(model_knn$fitted.values, col = c(3,2), 
     main = "Nº de datos predichos para cada clase por KNN")
cat("Error obtenido en la clase <=50K: ", 
    94700/length(adult.test.clean$income[adult.test.clean$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    154100/length(adult.test.clean$income[adult.test.clean$income==">50K"]))
```

A pesar de que el error en la clase `>50K` sigue siendo bastante alto, ha disminuido bastante y el error en la clase `<=50K` ha aumentado pero aún sigue siendo bajo. El error en ambas clases se ha equilibrado un poco y por eso, hemos obtenido un mejor error en el modelo que con `RandomForest`.

# Ajustando un modelo de _Red Neuronal_

## Idoneidad del modelo de _Red Neuronal_ para los datos del problema

La _Red Neuronal_ es un modelo muy extendido y con mucha flexibilidad en cuanto a la parametrización del modelo. Es tal la flexibilidad que ofrece, que tiene una gran potencia para ajustar los datos y capacidad para crear un modelo muy competitivo, pero tienen un alto coste a pagar, y es que el sobreajuste a los datos está practicamente servido de entrada. Este sobreajuste se puede mejorar usando regularización con _early stopping_ o _weight decay_.

Uno de los parámetros que utiliza la _Red Neuronal_ es el número de _unidades ocultas_, $m$. Este parámetro controla el número de capas internas de las que dispone nuestra red. Éste parámetro afecta mucho a su comportamiento, ya que tal y como se indica en [@lfd], con éste parámetro suficientemente alto podemos asegurar $E_{in} \approx E_{out}$. A mayor número de capas, las transformaciones que existen entre capa y capa, dan potencia y flexibilidad al modelo para generar una función $g$ que ajuste a $f$, con la contrapartida de que aumentamos la posibilidad de sobreajustar los datos. Además, como también tenemos un número alto de instancias, $N$, incrementamos la posibilidad de obtener un buen rendimiento en el modelo obtenido. La mejor forma de escoger este parámetro es mediante validación cruzada.

Al tener un problema de dimensión media/alta ($d=15$) y tantas instancias ($N = 30162$), conviene usar un modelo lo suficientemente flexible para poder obtener una buena aproximación de forma sencilla. 

Por último, al igual que otros modelos usados en este proyecto, las _Redes Neuronales_ son buenos modelos para clases desbalanceadas, tal y como se indica en [@unbalanced].

### Dimensión VC del modelo

La dimensión $d_{VC}$ para una red neuronal sigmoidal, generalmente, puede ser infinita. Pero, para el sigmoide $tanh(\cdot)$, donde el nodo salida es el $sign(\cdot)$, se puede ver que $$d_{VC} = O(VQ)$$ donde $V$ es el número de capas ocultas y $Q$ el número de pesos. 

## Normalización de los datos

Para normalizar los datos, vamos a usar el mismo resultado usado en el modelo del _Vecino más cercano_. 

## Parámetros e Hiperparámetros del modelo

Para elegir los parámetros, vamos a seguir los consejos dados en [@estadisticosinr]:

* ___Pesos iniciales___: los pesos iniciales serán números aleatorios cercanos a cero. Así, el modelo empezará de forma lineal pero se transformará en uno no lineal conforme los pesos aumenten. Si se empezase con pesos exactamente cero, el algoritmo no sería capaz de moverse, ya que estaríamos haciendo derivadas por cero. Empezar con pesos muy grandes suele llevar a soluciones malas.

* ___Regularización___: al tener tantos pesos, las redes neuronales pueden sobreajustar los datos en un mínimo local. Para evitar esto, como hemos dicho antes, podemos usar _early stopping_ o _weight decay_. En este caso vamos a usar _weight decay_, ya que la función que vamos a utilizar, `nnet`, nos permite ajustar este hiperparámetro de forma sencilla. Debido a la gran cantidad de instancias de las que disponemos, la función `tune.nnet` necesitaba un gran tiempo de cálculo y, en vez de realizar validación cruzada, hemos usado $\lambda = 0,1$ para aplicar regularización, ya que de entre todos los parámetros que hemos probado era el que mejor resultado obtenía.

* ___Número de capas ocultas___: es recomendable tener entre 5 y 100 capas, para que el modelo pueda tener flexibilidad suficiente. La regularización se encarga de dejar a 0 aquellos pesos que lleven a sobreajusar el modelo. El problema es que, debido a que el paquete sólo permite usar una única capa oculta, no vamos a poder configurar este parámetro y obtendremos unos resultados peores que si usáramos un mayor número de capas ocultas. Hay otros paquetes que sí permiten configurar el número de capas ocultas, pero no permiten usar datos categóricos.

* ___Número de unidades en la capa oculta___: lo que sí podemos configurar es el número de unidades en la capa oculta, aunque de forma limitada: si establecemos un número mayor a 15, la función devuelve un error diciendo que el modelo tiene demasiadas unidades.

## Aplicación del modelo

Hemos establecido a 10000 el número máximo de iteraciones, para que el algoritmo pudiese conveger antes de llegar a dicho límite. Esto no es un problema ya que estamos aplicando regularización con _weight decay_. El parámetro `trace = F` sólo sirve para que no se imprima cada paso que da el algoritmo en el PDF.

```{r}
library(nnet)
set.seed(1)
model.nnet = nnet(formula = income ~ ., maxit=10000, 
    data = adult.train.clean.norm, size = 10, decay=0.1,trace=F)
outOfSampleError(model.nnet, adult.test.clean.norm, type="class")
```

A pesar de disponer de una única capa oculta, el modelo ha sido capaz de obtener un resultado muy competitivo: $E_{out} = 15$%. Además, el número de fallos en cada clase está bastante equilibrado, casi tanto como en el modelo KNN con base radial:

```{r}
cat("Error obtenido en la clase <=50K: ", 
    94700/length(adult.test.clean$income[adult.test.clean$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    154100/length(adult.test.clean$income[adult.test.clean$income==">50K"]))
```

Por tanto, la conclusión obtenida sobre éste modelo es, que a pesar de haber estado limitado por la implementación del paquete de `R`, hemos obtenido un resultado muy competitivos y de una alta calidad, tanto en cuanto a error fuera de la muestra como error en cada clase y equilibrio en clases desbalanceadas.

# Comparación de los modelos obtenidos

Para concluir con los resultados obtenidos, vamos a hacer un estudio comparativo de todos los modelos que hemos ajustado en este proyecto. Para ello, vamos a representar la curva ROC de cada uno y realizar una tabla comparativa.


|                   | Regresión logística | Random forest | Support Vector Machine | Vecino más cercano | Red Neuronal |
|:-----------------:|:-------------------:|:-------------:|:----------------------:|:------------------:|:------------:|
|     $E_{out}$     | 24,56%              | 18,1%         |  14,78%                |  16,52%            | 15,69%       |
|  $E_{out >50K}$   | 100%                | $72,95$%      |  40,9%                 |  41,65%            | 37,59%       |
| $E_{out <=50K}$   | 0%                  | $0,24$%       |  6,29%                 |  8,34%             | 8,56%        |

Table: Comparativa de los resultados de los modelos ajustados


En la Tabla 1 vemos una comparativa de los modelos obtenidos. El modelo con menor $E_{out}$ ha sido _Support Vector Machine_ utilizando un $Kernel$ de base radial. Es el modelo con menor tasa de error en la clase minoritaria, ya que a diferencia de otros modelos, no ha "obviado" la clase minoritaria, aunque ha sido superado por la _Red Neuronal_, que lo ha conseguido con un error menor al 40%. Este modelo ha conseguido balancear mejor el error entre ambas clases, con lo que puede generalizar un poco mejor que el modelo $SVM$, con el coste de aumentar un poco el error fuera de la muestra. Ahora bien, el modelo _Red Neuronal_ ha jugado con desventaja, al no poder ajustar uno de los parámetros que más flexibilidad da al modelo: el número de capas ocultas, pero esto se debe a una limitación de la función del lenguaje. A pesar de esta desventaja, ha logrado ser el segundo mejor modelo.

El peor modelo ha sido la _Regresión logística_, ya que ha obviado la clase minoritaria completamente. Éste modelo ha obtenido un 24% de error, ya que ha clasificado de forma correcta todos los ejemplos de la clase `<=50K` y sólo ha fallado en los de la clase `>50K`. Al obviar la clase minoritaria, el error fuera de la muestra no se ve demasiado afectado, pero, puede darse que esta clase minoritaria, fuera de nuestro conjunto de train y test tenga una representación mucho mayor y deje de ser minoritaria, con lo que el modelo estaría fallando enormemente. El modelo _Random Forest_ también ha hecho una clasificación bastante mala, ya que el error obtenido en la clase `>50K` ha sido demasiado alto. Por tanto, al no haber podido trabajar correctamente con datos desbalanceados estos dos modelos han obtenidos los peores resultados.

Por último, el modelo del _Vecino más cercano_ con base radial ha obtenido buenos resultados, pero en comparación con _Support Vector Machine_ o _Red Neuronal_ ha obtenido un error mayor fuera de la muestra, debido a un mayor error en la clase minoritaria. Ha obtenido un error menor en `<=50K` que _Red Neuronal_, 8,34 frente a 8,65, pero esta diferencia no es significativa.

## Curvas ROC de los modelos obtenidos

A continuación, podemos ver las curvas $ROC$ asociadas a cada uno de los modelos estudiados:

```{r}
library(ROCR)

getPerfomance <- function(model, newdata = adult.test.clean, svmPred = F, 
                          nnetOrGlm=F, calculatePred = T,...){
    set.seed(1)
    if(calculatePred) {
        if(!nnetOrGlm && !svmPred)
            preds = predict(object = model, newdata = newdata, ...)[,2]
        else if (!svmPred && nnetOrGlm)
            preds = predict(object = model, newdata = newdata, ...)
        else{
            preds = attributes(predict(model, newdata, decision.values=T))$decision.values
            preds = preds*-1
        }
    } else
        preds = as.numeric(model)
    
    pred = prediction(preds, newdata$income)
    # tpr --> True Positive Rate
    # fpr --> False Positive Rate
    performance(pred, "tpr", "fpr")
}
plot(getPerfomance(model = glmModel, type="response", nnetOrGlm = T),
     col=2,lwd=2,main="Curvas ROC Para los distintos modelos estudiados") 
plot(getPerfomance(model = rf.tuned, type = "prob"),col=3,lwd=2,add=T)
plot(getPerfomance(model = svmModelReg_09, svmPred = T),lwd=2,col=4,add=T)
plot(getPerfomance(model_knn$prob[,2], newdata = adult.test.clean.norm, 
                   calculatePred = F),lwd=2,col=5,add=T)
plot(getPerfomance(model = model.nnet, newdata = adult.test.clean.norm,
  type="raw", nnetOrGlm = T),lwd=2,col=6,add=T)
abline(a=0,b=1,lwd=2,lty=2,col="gray")
legend("bottomright",col=c(2:6),lwd=2,legend=c("Regresión logística",
            "Random Forest","SVM","KNN","Red Neuronal"),bty='n')
```

Como podemos ver, las curvas que más destacan, a pesar de estar todas muy igualadas, son las curvas del modelo $SVM$ y de la red neuronal frente al resto, lo que nos indica, que el ratio de aciertos es mayor en estos que en el resto de modelos.

Respecto a la curva del modelo de regresión logística, la tasa de aciertos de este modelo es "bastante alta" ya que el modelo, al obviar la clase minoritaria, ha conseguido acertar la clase mayoritaria. Al ser clase mayoritaria, la tasa de aciertos es muy alta frente a los otros modelos, aunque esto es "virtual".

# Selección del mejor modelo

En base a los resultados, y el análisis de los errores que obtiene cada modelo, junto con el estudio de la curva ROC para cada modelo, podemos elegir como mejores modelos al modelo $SVM$ y el modelo de \textit{Red Neuronal}, ya que son los que mejores resultados ofrecen.

Respecto a la hora de elegir entre estos dos modelos, podemos elegir el modelo más sencillo, ya que, si seguimos el criterio de la \textit{Navaja de Ockham}, si tenemos dos modelos que se comportan de forma similar, se debe elegir el modelo más "simple". Tal y como hemos visto, la dimensión $d_{VC}$ de los modelos $SVM$ puede llegar a ser infinita, mientras que para los modelos de redes neuronales, esta dimensión se establece como $d_{VC} = O(VQ)$, por lo que, la red neuronal es un modelo un poco más sencillo en este caso que el modelo $SVM$, a pesar de que el error es un poco peor. Otro aspecto a tomar en cuenta en la elección del modelo es la facilidad para poder explicar lo que hace este modelo para predecir los datos, y por ello en este aspecto gana el modelo $SVM$.

Por todo esto, el mejor modelo elegido es \textbf{el modelo de \textit{Red Neuronal}}.

\newpage

# Ajustando un modelo de red neuronal recuperando los datos perdidos

## Recuperando datos perdidos

La idea de esto es recuperar todas las instancias de datos que hemos eliminado al contener valores perdidos. La forma de recuperar estos datos va a consistir en tratar de predecir esos datos con un modelo, y una vez predichos, ver cómo se comporta el nuevo modelo con todo el conjunto de datos al completo frente al que teníamos sin los datos perdidos.

Como hemos dicho antes, las variables que tienen valores perdidos son _Workclass_, _Occupation_ y _Country_. Vamos a hacer tres modelos diferentes de red neuronal, esta vez con regularización _early stopping_ (usando un límite de cien iteraciones).

<!-- ```{r} -->
<!-- set.seed(1) -->
<!-- # normalizamos todos los datos de la clase -->
<!-- norm = normalizar_maxmin(data = adult.train, data_test = adult.test) -->
<!-- adult.train.norm = norm[[1]] -->
<!-- adult.test.norm = norm[[2]] -->
<!-- norm = NULL -->

<!-- # debemos eliminar las columnas con NA para que funcione la función predict -->
<!-- workclass.nnet = nnet(formula = workclass ~ ., maxit=100,  -->
<!--         data = subset(adult.train.clean.norm, select=c(-occupation, -country)), -->
<!--         size = 8, decay=1, trace=F) -->
<!-- pred_workclass = predict(object = workclass.nnet,  -->
<!--             newdata = subset(adult.train.norm[is.na(adult.train$workclass),],  -->
<!--             select=c(-occupation, -country)), type="class") -->
<!-- pred_test_workclass = predict(object = workclass.nnet,  -->
<!--             newdata = subset(adult.test.norm[is.na(adult.test$workclass),], -->
<!--             select = c(-occupation, -country)), type = "class") -->
<!-- adult.train$workclass[is.na(adult.train$workclass)] = pred_workclass -->
<!-- adult.train.norm$workclass[is.na(adult.train.norm$workclass)] = pred_workclass  -->
<!-- adult.test$workclass[is.na(adult.test$workclass)] = pred_test_workclass -->
<!-- adult.test.norm$workclass[is.na(adult.test.norm$workclass)] = pred_test_workclass -->

<!-- occupation.nnet = nnet(formula = occupation ~ ., maxit=100,  -->
<!--         data = subset(adult.train.clean.norm, select=c(-workclass, -country)), -->
<!--         size = 8, decay=1, trace=F) -->
<!-- pred_occupation = predict(object = occupation.nnet,  -->
<!--             newdata = subset(adult.train.norm[is.na(adult.train$occupation),],  -->
<!--             select=c(-workclass, -country)), type="class") -->
<!-- pred_test_occupation = predict(object = occupation.nnet,  -->
<!--             newdata = subset(adult.test.norm[is.na(adult.test$occupation),],  -->
<!--             select=c(-workclass, -country)), type="class") -->
<!-- adult.train$occupation[is.na(adult.train$occupation)] = pred_occupation -->
<!-- adult.train.norm$occupation[is.na(adult.train.norm$occupation)] = pred_occupation -->
<!-- adult.test$occupation[is.na(adult.test$occupation)] = pred_test_occupation -->
<!-- adult.test.norm$occupation[is.na(adult.test.norm$occupation)] = pred_test_occupation -->

<!-- country.nnet = nnet(formula = country ~ ., maxit=100,  -->
<!--         data = subset(adult.train.clean.norm, select=c(-workclass, -occupation)), -->
<!--         size = 8, decay=1, trace=F) -->
<!-- pred_country = predict(object = country.nnet,  -->
<!--             newdata = subset(adult.train.norm[is.na(adult.train$country),],  -->
<!--             select=c(-workclass, -occupation)), type="class") -->
<!-- pred_test_country = predict(object = country.nnet,  -->
<!--             newdata = subset(adult.test.norm[is.na(adult.test$country),],  -->
<!--             select=c(-workclass, -occupation)), type="class") -->
<!-- adult.train$country[is.na(adult.train$country)] = pred_country -->
<!-- adult.train.norm$country[is.na(adult.train.norm$country)] = pred_country -->
<!-- adult.test$country[is.na(adult.test$country)] = pred_test_country -->
<!-- adult.test.norm$country[is.na(adult.test.norm$country)] = pred_test_country -->
<!-- ``` -->

```{r}
set.seed(1)
# normalizamos todos los datos de la clase
norm = normalizar_maxmin(data = adult.train, data_test = adult.test)
adult.train.norm = norm[[1]]
adult.test.norm = norm[[2]]
norm = NULL

# debemos quitar las variables con NA para poder predecir correctamente
predice_na <- function(attr, noselectvars, formula, test_dataset = adult.test.norm,
                       train_dataset = adult.train.norm, dataset = adult.train.clean.norm) {
    # el vector no selectvars siempre tendrá dos variables string
    model <- nnet(formula = formula, maxit = 100,
                  data = subset(dataset, select=c(-which(colnames(dataset) == noselectvars[1]),
                                                  -which(colnames(dataset) == noselectvars[2]))),
                  size = 8, decay = 1, trace = F)
    pred_train <- predict(object = model, newdata = subset(train_dataset[is.na(train_dataset[,attr]),],
                         select=c(-which(colnames(train_dataset) == noselectvars[1]),
                                  -which(colnames(train_dataset) == noselectvars[2]))), type = "class")
    pred_test <- predict(object = model, newdata = subset(test_dataset[is.na(test_dataset[,attr]),],
                         select=-which(colnames(test_dataset) == noselectvars)), type = "class")
    list(pred_train, pred_test)
}

workclass = predice_na(attr = "workclass", noselectvars = c("country", "occupation"), 
                       formula = workclass ~ .)
adult.train$workclass[is.na(adult.train$workclass)] = workclass[[1]]
adult.train.norm$workclass[is.na(adult.train.norm$workclass)] = workclass[[1]]
adult.test$workclass[is.na(adult.test$workclass)] = workclass[[2]]
adult.test.norm$workclass[is.na(adult.test.norm$workclass)] = workclass[[2]]

occupation = predice_na(attr = "occupation", noselectvars = c("country", "workclass"),
                        formula = occupation ~ .)
adult.train$occupation[is.na(adult.train$occupation)] = occupation[[1]]
adult.train.norm$occupation[is.na(adult.train.norm$occupation)] = occupation[[1]]
adult.test$occupation[is.na(adult.test$occupation)] = occupation[[2]]
adult.test.norm$occupation[is.na(adult.test.norm$occupation)] = occupation[[2]]

country = predice_na(attr = "country", noselectvars = c("occupation", "workclass"),
                     formula = country ~ .) 
adult.train$country[is.na(adult.train$country)] = country[[1]]
adult.train.norm$country[is.na(adult.train.norm$country)] = country[[1]]
adult.test$country[is.na(adult.test$country)] = country[[2]]
adult.test.norm$country[is.na(adult.test.norm$country)] = country[[2]]
```

Una vez hecho esto, vemos que ya no tenemos ningún valor perdido en los datos: 

```{r}
apply(X=adult.train, MARGIN=2, FUN=function(columna) length(is.na(columna)[is.na(columna)==T]))
apply(X=adult.test, MARGIN=2, FUN=function(columna) length(is.na(columna)[is.na(columna)==T]))
```

## Ajustando una red neuronal con el dataset completo

Una vez hemos conseguido predecir todos los valores perdidos, vamos a ajustar el mismo modelo de _Red Neuronal_ ajustado anteriormente para ver la diferencia obtenida al usar más datos.

```{r}
library(nnet)
set.seed(1)
model.nnet = nnet(formula = income ~ ., maxit=10000, 
    data = adult.train.norm, size = 10, decay=0.1,trace=F)
outOfSampleError(model.nnet, adult.test.norm, type="class")
```

En este caso, el $E_{out}$ ha mejorado algo respecto del modelo anterior, bajando del 15% de error al 14. Respecto al error obtenido en cada clase, vemos en la matriz de confusión que también se ha mejorado el error en ambas clases:

```{r}
cat("Error obtenido en la clase <=50K: ", 
    97200/length(adult.test$income[adult.test$income=="<=50K"]))
cat("Error obtenido en la clase >50K: ", 
    143800/length(adult.test$income[adult.test$income==">50K"]))
```

En ambas hemos reducido el error en un 1% y se puede comprobar como el error del modelo se reduce al tener más instancias para poder entrenar el modelo y a la vez, el tener más instancias para datos de test nos da una mejor representación del error real que tenemos fuera de la muestra. Aun así, tras haber recuperado los datos perdidos, no existe una mejora sustancial en el modelo, ya que los datos que eliminamos para aprender el modelo inicial sólo suponían un 7.37% del total de los datos que disponíamos, pero aun así, se puede notar cómo el modelo se comporta mejor que antes.

# Referencias