REGRESIONES NO LINEALES

Los modelos lineales tienen la ventaja de ser fácilmente interpretables, sin embargo, pueden tener limitaciones importantes en capacidad predictiva. Esto se debe a que, la asunción de linealidad, es con frecuencia una aproximación demasiado simple para describir las relaciones reales entre variables. A continuación, se describen métodos que permiten relajar la condición de linealidad intentando mantener al mismo tiempo una interpretabilidad alta.

Regresión polinómica: Consigue añadir curvatura al modelo introduciendo nuevos predictores que se obtienen al elevar todos o algunos de los predictores originales a distintas potencias.

Step functions: Se divide el rango del predictor en K subintervalos de forma que, en cada uno, se emplean únicamente las observaciones que pertenecen a la región para ajustar el modelo.

Regression splines: Se trata de una extensión de la regresión polinómica y de las step functions que consigue una mayor flexibilidad. Consiste en dividir el rango del predictor X en K subintervalos. Para cada una de las nuevas regiones se ajusta una función polinómica, introduciendo una serie de restricciones que hacen que los extremos de cada función se aproximen a los de las funciones de las regiones colindantes.

Smoothing splines: El concepto es similar a regression splines pero consigue la aproximación de los extremos de las funciones colindantes de forma distinta.

Local regression: Se asemeja a regression splines y smoothing splines en cuanto a que también se realizan ajustes por regiones, pero en este método las regiones solapan las unas con las otras.

Generalized additive models: Es el resultado de extender los métodos anteriores para emplear múltiples predictores.

La forma más sencilla de incorporar flexibilidad a un modelo lineal es introduciendo nuevos predictores obtenidos al elevar a distintas potencias el predictor original.

Partiendo del modelo lineal

\[ y_i = \beta_0 + \beta_1x_i + \epsilon_i \]

Se obtiene un modelo polinómico de grado d a partir de la ecuación

\[ y_i = \beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i+ \epsilon_i \]

Los modelos polinómicos se pueden ajustar mediante regresión lineal por mínimos cuadrados ya que, aunque generan modelos no lineales, su ecuación no deja de ser una ecuación lineal con predictores $ x, x^2, x^3, …, x^d $ Por esta misma razón, las funciones polinómicas pueden emplearse en regresión logística para predecir respuestas binarias. Solo es necesario realizar una transformación logit.

\[ P(y_i>Y|x_i=X) = \frac{exp(\beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i)}{1 + exp(\beta_0 + \beta_1x_i + \beta_2x^2_i + \beta_3x^3_i + ... + \beta_dx^d_i)} \]

En el libro Introduction to Statistical Learning desaconsejan el uso de modelos polinómicos con grado mayor de 3 o 4 debido a un exceso de flexibilidad (overfitting), principalmente en los extremos del predictor X. La selección del grado de polinomio óptimo puede hacerse mediante cross validation.

El set de datos Wage del paquete ISRL contiene información sobre 3000 trabajadores. Entre las 12 variables registradas se encuentra el salario (wage) y la edad (age). Dada la relación no lineal existente entre estas dos variables, se recurre a un modelo polinómico de grado 4 que permita predecir el salario en función de la edad.

La función lm() permite ajustar modelos lineales por mínimos cuadrados. En el argumento formula se especifica la variable dependiente y los predictores, que en este caso son wage ~ age + age^2 + age^3 + age^4. lm() tiene flexibilidad a la hora de interpretar el argumento fórmula, por lo que existen diferentes formas de obtener el mismo ajuste. Todas las que se muestran a continuación son equivalentes:

lm(wage ~ age + I(age^2) + I(age^3) + I(age^4), data = Wage)

lm(wage ~ cbind(age + age^2 + age^3 + age^4), data = Wage)

lm(wage ~ poly(age, 4, raw = TRUE), data = Wage)

lm(wage ~ poly(age, 4), data = Wage

Con la función poly() se puede generar directamente un polinomio, evitando tener que escribir toda la fórmula. Cuando se especifica que su argumento raw = TRUE, devuelve una matriz formada por el valor de la variable original elevada a cada una de las potencias del polinomio. Si el argumento raw = FALSE se devuelve una matriz con polinomios ortogonales, en la que cada columna es una combinación lineal de las otras. Es importante tenerlo en cuenta al emplear la función poly() en el ajuste de modelos, ya que, aunque no influye en el valor de las predicciones que se obtienen (las curvas obtenidas son iguales) sí que cambia el valor estimado de los coeficientes.

library(pacman)
p_load("ISLR", "ggplot2", "xfun", "cluster")
data("Wage")
Wage$age
##    [1] 18 24 45 43 50 54 44 30 41 52 45 34 35 39 54 51 37 50 56 37 38 40 75 40
##   [25] 38 49 43 34 57 18 55 51 33 34 36 56 70 25 32 27 28 27 43 50 39 52 35 57
##   [49] 25 33 57 71 43 23 30 22 59 28 61 34 43 54 69 41 48 49 42 37 55 21 58 31
##   [73] 25 32 40 44 60 23 63 44 47 61 55 24 42 25 34 53 53 70 47 46 33 34 22 74
##   [97] 40 45 43 33 62 37 54 34 50 46 41 63 38 35 29 66 37 39 42 51 55 51 38 49
##  [121] 42 43 38 59 57 25 49 41 38 61 49 52 43 60 46 21 61 32 58 35 26 32 37 22
##  [145] 51 44 35 60 40 35 35 47 43 33 60 38 53 55 57 64 43 35 54 45 58 48 46 46
##  [169] 55 51 49 34 53 40 50 37 39 52 50 48 47 27 39 44 37 52 26 39 25 31 58 30
##  [193] 27 40 55 35 48 29 25 40 27 44 49 22 45 33 63 49 39 25 29 37 35 46 58 39
##  [217] 41 29 37 62 27 37 56 36 23 46 26 43 36 54 62 40 53 48 23 34 43 55 59 60
##  [241] 49 45 31 34 26 22 40 41 36 38 48 47 47 35 49 40 32 57 56 73 29 49 30 29
##  [265] 36 39 48 45 50 42 58 48 48 45 41 27 44 32 46 52 22 53 34 38 42 25 29 28
##  [289] 46 55 56 41 61 47 52 27 33 58 27 52 45 22 25 25 41 54 23 36 31 52 46 45
##  [313] 40 33 42 41 25 24 24 51 41 42 30 55 23 60 40 41 80 53 40 35 55 24 39 30
##  [337] 68 54 27 34 53 48 59 30 45 39 45 43 30 57 53 42 23 58 43 33 56 37 41 33
##  [361] 50 43 53 37 20 26 29 57 47 54 38 23 31 43 33 48 33 52 30 44 38 51 47 28
##  [385] 52 23 65 27 33 49 27 48 38 58 40 46 53 25 42 40 35 50 34 48 49 41 37 58
##  [409] 28 33 50 38 48 30 49 44 50 58 45 44 35 59 45 51 61 45 34 33 45 50 50 39
##  [433] 26 49 20 31 49 56 61 56 51 42 45 30 48 33 50 31 48 35 47 45 52 46 47 49
##  [457] 36 51 46 37 23 43 34 58 44 61 40 37 33 38 41 30 51 39 56 44 51 35 27 47
##  [481] 51 41 63 61 46 33 35 28 50 41 23 45 63 59 30 30 45 42 57 52 54 30 47 38
##  [505] 47 39 50 47 45 44 52 42 18 56 50 50 53 38 27 56 50 56 50 49 42 50 53 28
##  [529] 45 24 25 38 47 38 33 39 52 46 42 61 32 42 61 27 45 58 25 34 44 40 51 50
##  [553] 64 54 42 31 25 32 51 75 28 28 32 53 36 43 51 35 46 38 43 39 37 59 33 23
##  [577] 37 40 25 27 47 52 23 33 37 34 51 40 50 64 54 50 32 37 19 32 42 26 46 18
##  [601] 43 62 44 40 43 32 48 35 61 47 55 28 52 27 22 20 42 29 31 59 36 32 32 25
##  [625] 43 44 57 49 41 38 58 59 39 60 56 44 70 27 47 41 29 42 51 47 48 39 23 31
##  [649] 59 19 42 35 38 45 58 26 30 54 38 44 52 73 65 28 53 53 52 27 57 32 51 43
##  [673] 50 53 54 24 40 30 45 55 45 44 37 41 44 52 55 43 27 43 46 30 43 48 33 29
##  [697] 24 52 50 24 42 53 49 46 40 41 41 36 34 24 44 33 25 40 48 56 25 44 27 58
##  [721] 52 60 39 51 31 45 43 45 40 46 51 47 52 31 40 26 49 38 50 63 35 29 33 29
##  [745] 62 45 56 58 39 49 24 61 58 47 41 71 47 38 31 56 44 33 36 31 55 48 50 42
##  [769] 25 55 28 30 46 48 42 23 60 34 45 28 33 58 47 24 56 56 23 56 45 58 30 46
##  [793] 48 50 43 37 49 31 54 52 70 30 44 34 38 28 38 39 52 34 33 45 58 38 63 38
##  [817] 31 53 21 38 18 28 33 32 29 28 42 28 36 47 53 48 36 55 37 68 65 37 33 19
##  [841] 50 22 40 45 36 37 24 66 35 33 44 48 47 41 55 74 33 62 56 29 56 47 45 44
##  [865] 25 45 40 56 56 47 39 38 48 44 42 55 61 39 30 49 33 31 35 55 36 37 33 48
##  [889] 47 60 28 28 48 29 40 32 33 40 42 36 39 44 43 44 29 49 47 46 41 52 43 38
##  [913] 42 56 39 29 33 33 41 44 45 58 29 42 68 47 30 39 42 49 36 44 54 42 39 30
##  [937] 33 28 38 30 25 52 26 54 62 57 32 60 40 49 54 48 35 46 42 58 30 27 55 50
##  [961] 55 40 32 27 49 50 50 44 57 72 23 31 36 63 48 46 55 44 44 32 60 42 48 36
##  [985] 31 54 29 46 35 50 25 33 34 51 33 59 62 28 40 52 29 34 25 66 43 45 35 43
## [1009] 49 60 46 55 32 53 57 51 59 45 56 30 21 60 60 48 39 32 51 47 44 46 28 24
## [1033] 37 31 72 50 36 27 27 41 20 49 25 41 53 45 32 47 26 31 32 50 46 47 55 49
## [1057] 25 60 62 63 31 29 51 28 36 28 38 35 46 37 41 59 51 26 34 32 23 50 25 43
## [1081] 36 50 34 34 55 39 40 32 63 38 67 39 31 41 47 61 61 21 22 53 47 36 36 42
## [1105] 54 27 36 35 60 42 41 35 34 22 47 59 29 32 49 30 33 51 37 46 63 41 40 46
## [1129] 45 48 28 65 50 53 48 30 37 51 55 35 33 41 34 33 49 62 31 43 43 43 48 48
## [1153] 50 29 52 46 40 43 54 38 32 28 33 44 41 31 48 50 47 31 40 54 24 34 65 28
## [1177] 30 45 35 46 47 71 39 41 33 76 20 34 26 51 25 62 37 37 32 55 36 39 46 48
## [1201] 34 39 50 54 60 48 47 20 60 37 46 54 51 32 35 46 38 62 47 28 61 26 49 28
## [1225] 61 42 31 50 44 46 54 40 38 43 43 30 58 60 45 47 47 59 36 39 40 42 29 38
## [1249] 58 18 35 30 48 27 40 34 54 20 50 42 45 59 44 67 26 40 42 22 43 34 44 53
## [1273] 54 55 30 48 48 30 38 25 59 38 40 36 60 40 28 34 39 26 28 27 34 38 51 58
## [1297] 42 32 27 51 50 25 60 62 53 39 23 47 32 57 40 57 31 28 58 47 41 53 27 19
## [1321] 47 44 32 23 31 42 50 52 20 33 26 40 29 35 34 37 37 22 39 42 44 45 42 52
## [1345] 27 46 59 52 19 64 38 45 53 28 39 52 41 44 48 47 31 28 45 56 47 41 49 28
## [1369] 19 48 47 57 39 45 52 51 46 45 40 48 47 45 49 57 77 48 56 39 44 49 58 46
## [1393] 39 26 61 62 47 40 35 61 38 18 45 27 22 29 50 25 35 41 38 49 36 66 23 28
## [1417] 42 20 37 58 45 64 45 52 42 32 40 55 46 55 53 33 45 54 32 42 35 54 38 31
## [1441] 41 45 49 50 40 29 25 32 51 29 32 42 59 31 29 52 35 28 50 49 39 26 31 50
## [1465] 62 35 22 55 63 63 39 28 43 58 38 49 62 36 54 26 46 58 22 58 63 67 44 51
## [1489] 46 43 40 31 41 21 51 43 34 45 35 33 46 38 44 50 55 49 47 40 21 60 56 61
## [1513] 44 43 29 38 26 28 38 45 39 33 44 48 45 56 40 71 30 45 41 40 59 46 38 26
## [1537] 53 46 42 28 30 34 46 22 37 27 52 43 52 54 47 21 42 72 38 53 35 60 39 45
## [1561] 61 64 39 46 25 80 51 58 52 45 41 42 31 49 47 39 53 37 52 41 39 63 23 62
## [1585] 43 29 43 34 38 25 64 49 46 30 33 29 32 54 40 24 30 49 47 25 44 41 54 59
## [1609] 25 63 49 35 24 66 57 47 44 56 34 23 30 38 43 58 36 48 35 22 50 58 20 44
## [1633] 50 55 39 32 50 24 53 35 25 24 57 43 33 40 39 32 48 29 35 41 52 50 35 27
## [1657] 57 48 41 59 55 54 34 43 43 62 45 40 33 45 21 54 43 54 35 23 61 51 43 20
## [1681] 51 65 37 24 36 45 33 36 32 45 51 33 47 37 41 23 57 61 47 29 37 57 39 59
## [1705] 31 26 51 38 51 50 32 38 49 40 60 30 41 32 26 41 53 35 50 43 50 19 42 53
## [1729] 22 49 33 31 43 57 40 48 46 58 50 38 33 52 48 53 37 43 54 41 57 54 25 35
## [1753] 29 36 43 41 52 48 51 55 37 48 42 43 32 33 48 20 43 58 38 23 54 66 47 35
## [1777] 39 35 40 41 40 50 35 43 26 23 61 24 66 50 22 62 55 41 54 29 39 71 22 48
## [1801] 35 46 51 55 36 57 21 60 41 56 32 32 36 28 60 41 55 36 46 62 44 36 41 22
## [1825] 51 57 55 36 43 38 27 32 44 56 38 36 28 39 42 51 55 39 42 33 23 40 40 51
## [1849] 76 57 27 26 25 27 33 51 22 47 34 47 31 46 47 49 40 51 31 58 48 54 56 61
## [1873] 32 60 42 22 46 53 43 32 30 50 40 55 42 35 33 33 38 41 42 58 45 47 30 34
## [1897] 43 58 44 28 46 36 37 47 56 38 41 62 40 44 56 51 41 52 27 52 51 51 28 49
## [1921] 61 39 23 53 55 48 50 31 53 49 49 43 20 46 48 19 56 44 23 42 32 32 43 57
## [1945] 26 67 31 49 29 41 35 39 33 53 43 36 38 58 45 34 30 21 39 28 54 24 57 48
## [1969] 31 47 44 30 41 26 32 39 42 24 44 30 58 33 56 28 41 22 37 18 37 43 49 36
## [1993] 47 52 34 62 44 48 51 35 42 34 50 50 55 56 50 44 30 34 50 51 53 40 34 36
## [2017] 44 54 36 62 33 44 80 45 58 51 60 33 28 37 56 59 56 34 37 36 45 36 23 56
## [2041] 33 48 28 37 56 22 60 30 50 38 55 26 39 42 71 33 49 45 41 59 33 47 34 23
## [2065] 32 30 54 46 52 50 57 46 30 32 53 48 31 50 41 47 37 24 53 39 47 36 34 45
## [2089] 60 43 20 27 35 50 56 19 35 23 45 54 43 42 51 49 48 50 37 38 58 41 51 51
## [2113] 44 43 34 43 59 61 40 49 55 54 41 29 66 61 34 41 48 26 54 44 60 57 53 54
## [2137] 37 62 40 39 55 57 45 41 28 58 47 60 26 52 31 63 42 30 52 38 29 42 57 46
## [2161] 21 34 39 26 42 48 54 80 37 31 41 37 43 54 39 45 21 59 29 29 33 45 53 63
## [2185] 56 21 39 30 42 40 47 43 42 50 39 50 60 57 38 49 51 33 59 41 33 28 49 51
## [2209] 36 26 22 47 59 41 32 45 39 64 36 61 34 52 33 36 34 39 59 49 40 49 25 59
## [2233] 39 40 51 63 53 25 53 23 36 19 47 43 40 38 63 45 32 37 37 41 45 22 58 44
## [2257] 39 40 33 35 63 45 31 41 66 45 47 35 52 62 47 40 41 48 55 64 29 49 32 41
## [2281] 58 41 34 43 23 40 46 35 42 47 45 34 57 32 55 22 51 26 27 34 48 23 41 49
## [2305] 47 35 50 55 69 45 40 47 22 42 47 48 28 49 51 58 32 55 49 46 44 25 73 18
## [2329] 52 28 51 44 44 29 48 62 40 24 48 54 69 48 73 30 41 48 63 29 40 52 54 33
## [2353] 42 31 46 44 45 40 31 19 43 56 55 59 49 42 27 34 47 35 33 35 40 53 47 45
## [2377] 34 25 56 54 41 46 44 44 46 38 36 27 33 56 52 33 32 39 22 53 36 37 40 51
## [2401] 40 56 50 49 32 21 48 48 32 30 30 33 56 46 49 24 48 37 53 56 34 56 48 47
## [2425] 52 40 25 35 40 22 50 48 40 45 30 46 34 32 53 32 39 51 40 66 40 52 36 32
## [2449] 45 37 47 49 53 66 33 58 41 54 26 24 62 63 26 45 23 31 44 43 48 48 55 45
## [2473] 48 63 57 69 31 38 43 19 44 31 49 29 26 37 30 35 40 44 65 54 56 48 37 30
## [2497] 46 37 58 52 45 33 51 50 42 46 43 47 47 44 36 39 29 40 20 32 36 31 28 41
## [2521] 54 42 19 59 64 31 37 30 24 26 44 39 48 53 42 44 25 46 52 33 30 38 29 57
## [2545] 50 23 46 61 70 59 57 25 48 52 44 40 43 51 43 42 28 54 34 50 36 36 71 31
## [2569] 43 20 23 42 41 27 29 41 41 22 39 50 49 42 38 55 47 52 56 30 56 46 43 35
## [2593] 58 69 68 52 46 39 66 30 30 43 49 54 48 37 46 53 31 18 49 39 40 37 24 35
## [2617] 27 44 46 40 39 42 37 58 40 27 31 44 38 34 41 20 28 55 38 33 43 50 55 24
## [2641] 41 37 42 25 45 40 43 50 42 38 44 19 44 26 54 32 49 35 34 23 60 54 31 47
## [2665] 48 48 35 55 35 37 27 64 52 47 26 41 22 76 46 39 39 23 49 44 20 43 42 46
## [2689] 33 62 38 43 46 18 63 43 44 53 58 23 50 35 38 36 32 56 41 40 26 26 47 25
## [2713] 58 57 30 44 29 51 37 67 34 50 38 55 27 39 54 55 49 42 39 32 39 50 47 40
## [2737] 51 28 53 54 51 33 55 44 30 40 26 49 43 43 44 32 34 36 51 42 60 36 47 60
## [2761] 22 41 56 44 49 45 25 52 30 29 23 48 44 62 33 51 36 23 51 47 34 34 43 43
## [2785] 40 59 53 32 48 37 20 32 40 48 43 56 55 38 52 40 29 45 39 50 50 57 30 56
## [2809] 43 43 67 50 49 50 59 51 39 54 55 45 54 23 74 52 25 45 37 55 25 66 47 36
## [2833] 63 44 37 32 43 37 26 39 45 40 51 52 43 52 50 54 28 28 26 45 34 73 31 29
## [2857] 40 38 42 50 51 40 25 47 27 53 67 40 44 33 29 32 27 40 36 48 56 24 53 34
## [2881] 61 38 55 70 24 22 38 29 30 38 49 71 32 40 32 33 40 25 36 46 59 49 41 42
## [2905] 59 27 30 54 31 34 47 40 29 58 41 56 52 53 37 30 35 40 54 22 37 39 62 47
## [2929] 51 47 28 65 42 58 32 34 49 39 44 38 48 47 40 27 42 30 31 29 53 37 48 56
## [2953] 44 63 25 31 51 25 50 34 20 30 40 34 57 42 30 37 50 26 59 29 22 54 46 51
## [2977] 35 49 53 61 40 52 40 56 39 30 58 33 51 32 50 26 35 31 31 44 30 27 27 55

Primera aproximación a los polinomios

poly(Wage$age, degree = 4, raw = TRUE, simple=TRUE)[1:5,]
##       1    2      3       4
## [1,] 18  324   5832  104976
## [2,] 24  576  13824  331776
## [3,] 45 2025  91125 4100625
## [4,] 43 1849  79507 3418801
## [5,] 50 2500 125000 6250000
poly(Wage$age, degree = 4, raw = FALSE, simple=TRUE)[1:5,]
##                  1            2             3           4
## [1,] -0.0386247992  0.055908727 -0.0717405794  0.08672985
## [2,] -0.0291326034  0.026298066 -0.0145499511 -0.00259928
## [3,]  0.0040900817 -0.014506548 -0.0001331835  0.01448009
## [4,]  0.0009260164 -0.014831404  0.0045136682  0.01265751
## [5,]  0.0120002448 -0.009815846 -0.0111366263  0.01021146

Cálculo del modelo polinómico de grado 4

modelo_poli4 <- lm(wage ~poly(age, 4), data = Wage )
summary(modelo_poli4)
## 
## Call:
## lm(formula = wage ~ poly(age, 4), data = Wage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -98.707 -24.626  -4.993  15.217 203.693 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    111.7036     0.7287 153.283  < 2e-16 ***
## poly(age, 4)1  447.0679    39.9148  11.201  < 2e-16 ***
## poly(age, 4)2 -478.3158    39.9148 -11.983  < 2e-16 ***
## poly(age, 4)3  125.5217    39.9148   3.145  0.00168 ** 
## poly(age, 4)4  -77.9112    39.9148  -1.952  0.05104 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.91 on 2995 degrees of freedom
## Multiple R-squared:  0.08626,    Adjusted R-squared:  0.08504 
## F-statistic: 70.69 on 4 and 2995 DF,  p-value: < 2.2e-16

El p-value obtenido para el estadístico F es muy bajo (< 2.2e-16), lo que indica que al menos uno de los predictores introducidos en el modelo está relacionado con la variable respuesta wage. Los p-values individuales de cada predictor son todos muy bajos a excepción de age^4, lo que apunta a que un polinomio de grado 3 es suficiente para modelar el salario en función de la edad. Acorde al R2, el modelo es capaz de explicar el 8.6% de la variabilidad observada en wage (es un % muy bajo).

Si se quiere generar una representación gráfica del modelo ajustado

\[ wage = 111.70 + 447.07age - 478.32age^2 + 125.52age^3 - 77.91age^4) \]

se tiene que representar la curva que describe su ecuación. Para conseguirlo, es necesario predecir el valor de la variable respuesta en puntos interpolados dentro del rango del predictor. Cuantos más puntos se interpolen, mejor será la representación. R permite obtener predicciones (junto con el intervalo de confianza) mediante la función predict().

# INTERPOLACIÓN DE PUNTOS DENTRO DEL RANGO DEL PREDICTOR
# ------------------------------------------------------
limites <- range(Wage$age)
nuevos_puntos <- seq(from = limites[1], to = limites[2], by = 1)
nuevos_puntos <- data.frame(age = nuevos_puntos) 

# PREDICCIÓN DE LA VARIABLE RESPUESTA Y DEL ERROR ESTÁNDAR
# --------------------------------------------------------
predicciones <- predict(modelo_poli4, newdata = nuevos_puntos, se.fit = TRUE, level = 0.95)

# CÁLCULO DEL INTERVALO DE CONFIANZA SUPERIOR E INFERIOR 95%
# ----------------------------------------------------------
intervalo_conf <- data.frame(inferior = predicciones$fit -
                                        1.96*predicciones$se.fit,
                             superior = predicciones$fit +
                                        1.96*predicciones$se.fit)

attach(Wage)
plot(x = age, y = wage, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: wage ~ age")
points(x = nuevos_puntos$age, predicciones$fit, col = "red", pch = 20)
points(x = nuevos_puntos$age, intervalo_conf$inferior, col = "blue", pch = 4)
points(x = nuevos_puntos$age, intervalo_conf$superior, col = "blue", pch = 4)

El ejemplo anterior muestra que, para conseguir el efecto visual de curva continua, es necesario interpolar muchos datos. La función lines() facilita el proceso uniendo los puntos con segmentos, lo que mejora la representación gráfica de curvas

attach(Wage)
## The following objects are masked from Wage (pos = 3):
## 
##     age, education, health, health_ins, jobclass, logwage, maritl,
##     race, region, wage, year
plot(x = age, y = wage, pch = 20, col = "darkgrey")
title("Polinomio de grado 4: wage ~ age")
lines(x = nuevos_puntos$age, predicciones$fit, col = "red", lwd = 2)
lines(x = nuevos_puntos$age, intervalo_conf$inferior, col = "blue", lwd = 2)
lines(x = nuevos_puntos$age, intervalo_conf$superior, col = "blue", lwd = 2)

El paquete gráfico ggplot2 permite obtener representaciones de modelos de forma rápida. Solo con especificar la fórmula del modelo, ejecuta automáticamente todos los cálculos necesarios (ajuste, predicción…).

ggplot(data = Wage, aes(x = age, y = wage)) +
  geom_point(color = "grey30", alpha = 0.3) + 
  geom_smooth(method = "lm", formula = y ~ poly(x, 4), color = "red") +
  labs(title = "Polinomio de grado 4: wage ~ age") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

Cuando se realiza regresión polinómica, se debe decidir qué grado de polinomio emplear. Cuanto mayor sea el polinomio más flexibilidad tendrá el modelo pero, a su vez, más riesgo de overfitting. Acorde al principio de parsimonia, el grado óptimo es el grado más bajo que permita explicar la relación entre ambas variables. Para identificarlo se puede recurrir a dos estrategias distintas: contraste de hipótesis o cross-validation.

Comparación de modelos por contraste de hipótesis ANOVA

Identificar el modelo polinómico más simple que permite explicar la relación entre variables equivale a identificar el grado de polinomio a partir del cual ya no hay una mejora significativa del ajuste. Cuando se comparan dos modelos anidados (el modelo de menor tamaño está formado por un subset de predictores del modelo mayor), se puede saber si el modelo mayor aporta una mejora sustancial estudiando si los coeficientes de regresión de los predictores adicionales son distintos a cero. El test estadístico empleado para hacerlo es el ANOVA.

\[ Modelo_{menor}: \ \ y = \beta_0 + \beta_1x_1 +...+ \beta_kx_k \] \[ Modelo_{mayor}: \ \ y = \beta_0 + \beta_1x_1 +...+ \beta_kx_k + \beta_{k+1}x_{k+1} + ... + \beta_{p}x_{p} \]

La hipótesis a contrastar es que todos los coeficientes de regresión de los predictores adicionales son igual a cero, frente a la hipótesis alternativa de que al menos uno es distinto.

\[ H_0: \beta_{k+1}= ... = \beta_{p} \]

El estadístico empleado es:

\[ F = \frac{(SEE_{Modelo_{menor}} - SEE_{Modelo_{mayor}})/(p-k)}{SEE_{Modelo_{mayor}}/(n-p-1)} \] Dado que un polinomio de orden n siempre va a estar anidado a uno de orden n+1, se pueden comparar modelos polinómicos dentro un rango de grados haciendo comparaciones secuenciales.

Se ajustan los modelos polinómicos de grado 1 a 5

modelo_1 <- lm(wage ~ age, data = Wage)
modelo_2 <- lm(wage ~ poly(age, 2), data = Wage)
modelo_3 <- lm(wage ~ poly(age, 3), data = Wage)
modelo_4 <- lm(wage ~ poly(age, 4), data = Wage)
modelo_5 <- lm(wage ~ poly(age, 5), data = Wage)

anova(modelo_1, modelo_2, modelo_3, modelo_4, modelo_5)
## Analysis of Variance Table
## 
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
##   Res.Df     RSS Df Sum of Sq        F    Pr(>F)    
## 1   2998 5022216                                    
## 2   2997 4793430  1    228786 143.5931 < 2.2e-16 ***
## 3   2996 4777674  1     15756   9.8888  0.001679 ** 
## 4   2995 4771604  1      6070   3.8098  0.051046 .  
## 5   2994 4770322  1      1283   0.8050  0.369682    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-value de la comparación entre el modelo lineal (modelo_1) y el cuadrático (modelo_2) es prácticamente cero (< 2.2e-16), indicando que el modelo lineal no es suficiente. La comparación entre el modelo cuadrático y el cúbico también ha resultado en un p-value muy pequeño (0.001679), evidencia suficiente de que el modelo cúbico es superior. El p-value obtenido al comparar el modelo cúbico con el de grado 4 está ligeramente por encima del 0.05 y el de comparar grado 4 con grado 5 es muy alto (0.37). En base a los resultados del ANOVA, el modelo cúbico es el mejor. Polinomios superiores no aportan una mejora significativa y polinomios inferiores pierden mucha capacidad de ajuste.

Es importante recordar que las comparaciones por ANOVA pueden hacerse entre cualquier par de modelos, siempre y cuando estos sean anidados

Combinación de Métodos (Asignación)

Relación entre Ingreso y Servicio de salud

Ahora, pasaremos a combinar métodos para determinar las personas que tienen o no servicio de salud, para esto utilizaremos K-Means.

En Estados Unidos el servicio de salud pública es inexistente, por lo que es necesario que cada persona pague por su propio servicio de salud, esto puede (en teoría) ser un problema para aquellas personas que tengan un bajo ingreso económico. Para demostrar lo anterior, analizaremos la variable cualitativa “health_ins” (Seguro de Salud) utilizando K-Means para determinar si las personas que ganan más dinero son las que tienen servicio de salud.

Partiendo de la hipótesis nula: “El salario no afecta en si la persona cuenta o no con servicio de salud”. Veremos al finalizar el análisis si es correcta o incorrecta.

Viendo los datos gráficamente

ggplot(Wage, aes(wage, age)) + geom_point(aes (col=health_ins), size=4   )

Por el momento podemos observar, respecto al salario y su edad los grupos de personas que cuentan con seguro médico y los que no se encuentran muy mezclados entre sí, aún así, podemos observar cómo a medida de que el salario aumenta, va disminuyendo la cantidad de personas que no cuentan con un seguro de vida.

K-Means

Ahora vamos a utilizar la función K-Means para construir los clusters

set.seed(95)
wageCluster <- kmeans(data.frame(wage, age), center=2, nstart = 20 )
wageCluster
## K-means clustering with 2 clusters of sizes 2261, 739
## 
## Cluster means:
##        wage      age
## 1  94.05367 41.44449
## 2 165.70429 45.38295
## 
## Clustering vector:
##    [1] 1 1 2 2 1 1 2 1 1 1 1 1 1 2 2 1 1 2 1 1 2 1 1 2 1 2 1 1 1 1 2 1 1 1 1 1 1
##   [38] 1 1 1 1 1 1 1 1 2 2 2 1 1 2 1 2 1 1 1 2 2 1 1 2 1 1 1 1 1 2 2 1 1 1 1 1 1
##   [75] 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1
##  [112] 1 1 1 1 1 1 2 2 1 1 1 2 2 1 1 2 1 1 1 1 2 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 2
##  [149] 2 1 1 2 1 1 2 2 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 2 2 1 2 1 2 1 1 1 1 1 2
##  [186] 1 1 1 1 1 2 1 1 1 2 2 1 2 1 1 1 1 1 1 1 1 2 2 2 1 2 1 2 1 1 1 2 1 1 1 1 1
##  [223] 1 2 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2
##  [260] 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 2 1 1 1
##  [297] 1 1 2 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 1
##  [334] 1 2 1 1 1 1 1 1 1 1 1 2 1 1 2 2 2 2 2 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 2 2
##  [371] 1 1 1 2 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [408] 1 1 1 2 1 1 1 2 1 1 1 2 1 1 1 2 2 2 2 1 1 1 2 1 1 2 1 1 1 1 1 2 1 1 1 2 1
##  [445] 1 1 1 2 2 2 1 2 1 1 1 2 1 1 2 2 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1
##  [482] 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 2
##  [519] 1 2 1 2 1 1 1 2 2 2 1 1 1 1 1 2 1 2 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 2 2 1
##  [556] 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 2 1 1 1 1
##  [593] 1 1 1 1 1 1 2 1 1 1 2 2 2 1 1 1 2 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 2 1 1 2 2
##  [630] 2 1 2 2 1 1 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1
##  [667] 2 1 2 1 2 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 2 1 2 1 1 2 1 2 1 1 1 2 2 1 1 1 2
##  [704] 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 2 1 2 2 2 2 2 1 1 1 1 2
##  [741] 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 2 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 2 1 1 1
##  [778] 1 1 1 2 1 2 1 2 1 1 2 2 2 1 2 2 1 1 1 2 2 1 2 1 2 1 2 1 1 2 1 1 2 1 1 1 1
##  [815] 1 1 1 2 1 2 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1
##  [852] 1 1 1 1 2 1 2 1 1 1 2 2 2 1 1 2 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1
##  [889] 1 1 1 1 2 1 1 1 1 1 2 1 1 2 1 1 1 1 2 2 1 1 1 2 1 1 1 1 2 1 1 1 2 1 1 1 1
##  [926] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1 2 2 1 2 1 2 2 1 1 1 1 2
##  [963] 1 1 1 1 1 2 2 1 1 2 1 1 1 1 2 2 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
## [1000] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 1 1 1 1
## [1037] 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1074] 1 1 1 1 2 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
## [1111] 1 2 2 1 1 1 1 1 1 1 1 2 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1
## [1148] 1 2 1 1 1 2 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2 1
## [1185] 2 2 1 2 1 2 1 1 2 1 1 2 1 2 2 1 2 1 1 1 2 1 1 1 1 1 2 1 2 1 2 1 1 1 1 2 2
## [1222] 1 1 1 1 2 1 1 1 2 1 1 1 2 1 1 1 1 2 2 1 1 1 1 1 1 1 2 1 1 2 2 1 1 1 1 2 1
## [1259] 1 2 1 1 1 1 1 2 1 1 2 1 1 2 1 2 1 2 1 1 1 1 1 2 2 2 1 1 1 2 2 1 1 1 1 2 1
## [1296] 1 2 1 1 1 2 1 1 2 1 2 1 2 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2
## [1333] 1 1 1 1 2 1 1 1 1 1 2 2 1 2 1 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1
## [1370] 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1
## [1407] 1 1 1 1 2 1 1 1 1 1 1 1 2 1 2 1 2 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1
## [1444] 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1
## [1481] 2 2 1 2 2 2 1 2 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1 1 2 2 2 1 1 1 1 2 2 1 2 1
## [1518] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 2 2
## [1555] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 1 2 1 1 1 1
## [1592] 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2 2 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1
## [1629] 1 1 1 1 2 1 1 2 2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1 1 1 2 1 1 1
## [1666] 1 1 2 1 2 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 2 1 1 1 1 2 1 2 1 1 1 2 1 1 1 2
## [1703] 1 2 2 1 1 1 1 1 1 2 2 1 2 1 2 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
## [1740] 1 1 2 1 2 1 1 1 2 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1
## [1777] 2 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1
## [1814] 1 2 1 1 1 1 1 1 1 2 1 2 1 1 2 2 1 1 1 1 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1
## [1851] 1 1 1 1 1 1 1 1 2 1 1 1 1 2 2 1 1 1 1 1 2 1 2 2 2 1 1 1 1 1 2 1 2 1 2 1 1
## [1888] 1 1 1 1 2 1 1 1 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2
## [1925] 1 1 1 1 1 2 2 1 1 1 2 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1
## [1962] 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1 2 1 1 2 1 1 1 2 1 1 1 1 2 1 2 1 2 1 2 1 2
## [1999] 1 2 1 1 1 2 1 2 2 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1
## [2036] 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 1 1 1 2 1 1 1 1 1 2 2 1 1
## [2073] 1 1 1 2 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2
## [2110] 1 2 1 1 1 2 2 2 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
## [2147] 2 1 2 1 1 2 1 1 1 2 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 1 2 2
## [2184] 1 2 1 1 2 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [2221] 1 2 1 2 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 1
## [2258] 2 1 1 1 2 1 1 2 1 1 1 1 2 1 2 2 1 1 1 1 1 2 1 2 2 1 2 1 1 2 1 2 1 2 1 1 1
## [2295] 1 1 1 1 1 2 1 1 1 2 1 1 1 2 1 2 2 2 1 1 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1
## [2332] 2 2 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 2 1 1 2 2 1 2 1 1 2 1 1 1 2 1 1 1 1 1 1
## [2369] 2 2 1 2 1 1 1 2 2 1 1 1 1 2 2 1 2 1 1 1 1 2 1 1 1 1 1 1 2 1 2 1 1 1 1 2 1
## [2406] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1
## [2443] 1 1 2 2 1 1 1 2 2 2 1 1 1 1 1 2 1 1 2 2 1 1 1 1 2 1 1 1 1 2 2 1 1 2 1 2 1
## [2480] 1 1 2 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1
## [2517] 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 2 2 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1
## [2554] 1 2 1 2 2 1 1 1 2 2 1 2 2 1 1 2 1 1 1 2 1 1 1 1 1 1 1 2 1 1 2 1 1 1 2 2 1
## [2591] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 1 2 2 1 1 2 1 1 1 1 2 1 2 1 1 2 1 1 1 2 1 1
## [2628] 1 1 2 2 1 2 1 1 1 1 1 1 1 2 1 2 1 2 1 2 2 1 2 1 1 1 1 1 1 2 1 1 1 1 1 2 2
## [2665] 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 2 1 1 2 2 1 2 1 2 1 2 1 1 1 1
## [2702] 1 2 1 2 2 1 1 1 2 1 1 2 1 1 2 1 1 1 2 1 1 2 2 1 2 2 2 2 1 1 1 1 2 1 1 1 1
## [2739] 1 1 1 2 1 1 1 1 1 1 2 1 1 1 2 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2
## [2776] 1 2 1 1 2 1 1 2 1 2 1 1 1 1 2 1 1 1 1 2 1 2 1 1 1 1 2 1 2 2 2 1 1 2 1 1 1
## [2813] 1 1 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 2 2 1 2 2 2 1 2 1 1 1 2 2 2 1 1 1 2 1 1
## [2850] 2 1 1 1 1 1 1 2 2 2 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1
## [2887] 1 1 2 2 1 1 1 2 1 1 1 1 1 2 1 2 1 2 1 1 1 1 1 2 1 2 1 1 1 1 2 1 1 1 2 1 1
## [2924] 1 2 2 1 1 2 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1
## [2961] 1 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1 1 2 1 1 2 1 1 1 2 2 1 2 2 1 1 1 2 2 1
## [2998] 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 1342984 1410685
##  (between_SS / total_SS =  51.0 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Comparando estos clusters con los datos originales

table(wageCluster$cluster, Wage$health_ins)
##    
##     1. Yes 2. No
##   1   1437   824
##   2    646    93

Agrupando los datos en clusters

clusplot(data.frame(wage, age), wageCluster$cluster, color=T, shade=T, lines=0)

Al ver los datos agrupados y divididos en sus respectivos clusters podemos observar que existe mucho ruido entre la transición entre tener y no tener seguro de vida en base al salario, pero se puede notar una clara diferencia entre los salarios muy altos y muy bajos.

Conclusión 1

En conclusión, una vez realizado el análisis utilizando k-means podemos observar que sí existe una relación entre el salario y el tener o no seguro médico, pero este no es un factor determinante, pues si vemos a nuestra primer gráfica, podemos ver que los datos se encuentran muy revueltos entre sí, por lo que sólo nos podemos dar una idea de en que rangos las personas si tienen seguro médico y en cuales no, esto se debe a que en Estados Unidos no existe un seguro gratuito como tal, lo que hace que tener o no seguro medico se vea influido por el sueldo percibido, pero también por otros factores que no son tomados en cuenta en este análisis con k-means.

Relación entre el Ingreso percibido y la Raza

¿Es posible saber tu raza con sólo tu salario?

Como vimos en el análisis anterior, utilizando k-means es difícil determinar a partir del salario si una persona tiene o no seguro médico, entonces, ¿que si podemos saber a partir del salario? En este apartado vamos a análizar con k-means si a partir del ingreso percibido, es posible conocer tu raza; partimos de la hipótesis “No hay relación entre el ingreso percibido y la raza a la que perteneces en EU”.

Sabemos que en Estados Unidos el racismo está muy marcado si no eres de tez blanca, esto debido a la historia que tiene el país, en dónde no hubo un mestizaje al momento de la conquista y este problema ha sido acarreado por generaciones, en la última década se ha intentado reducir el racismo y la separación de las razas en EU, así que a partir de los datos que tenemos, vamos a comprobar si ha surtido efecto o aún hay una diferencia muy marcada entre las razas.

Viendo los datos gráficamente

ggplot(Wage, aes(wage, age)) + geom_point(aes (col=race), size=4   )

K-Means

Ahora vamos a utilizar la función K-Means para construir los clusters, en este caso, nosotros sabemos que son 4 clusters, así que vamos a usar la función con 4 clusters para ver el resultado y al final utilizaremos la ley del codo para conocer si 4 clusters son suficientes.

input <- data.frame(wage, age)

set.seed(53)
raceCluster <- kmeans(input, centers=4, nstart = 20 )
raceCluster
## K-means clustering with 4 clusters of sizes 80, 1084, 571, 1265
## 
## Cluster means:
##        wage      age
## 1 276.90658 46.88750
## 2  76.07952 38.64207
## 3 155.39608 45.21016
## 4 112.06082 44.10277
## 
## Clustering vector:
##    [1] 2 2 4 3 2 4 3 4 4 4 4 2 2 3 3 2 2 3 4 4 3 2 2 3 2 1 2 2 4 2 3 4 2 4 2 2 2
##   [38] 4 4 2 2 4 2 2 2 3 1 3 2 4 3 4 1 2 4 2 3 3 4 4 3 4 2 2 4 4 3 3 4 2 4 2 2 4
##   [75] 2 4 1 2 4 2 4 4 2 2 2 2 2 2 4 4 4 4 2 4 2 4 4 3 3 3 4 4 4 2 4 4 2 2 4 4 4
##  [112] 2 2 4 2 4 4 3 1 4 2 2 3 3 4 2 4 2 4 2 4 4 3 4 3 2 4 4 2 2 2 2 4 2 2 4 3 4
##  [149] 3 2 2 3 2 2 3 3 2 4 4 2 4 4 2 2 3 4 2 2 2 4 2 2 3 3 3 4 3 4 3 4 4 4 4 4 3
##  [186] 4 4 4 4 2 3 4 2 4 3 3 4 3 2 4 2 4 4 2 4 2 1 3 4 2 4 2 3 4 2 4 3 2 4 2 2 2
##  [223] 4 3 2 4 2 4 4 4 4 4 4 3 2 2 1 2 2 2 2 4 2 4 4 2 2 4 4 2 4 3 2 4 2 4 4 2 3
##  [260] 4 2 3 4 2 4 3 4 4 2 2 4 4 4 4 4 2 4 4 4 2 2 4 4 3 4 2 2 3 2 4 4 4 3 2 2 2
##  [297] 2 2 3 2 2 2 2 2 4 3 2 2 2 4 2 4 2 2 3 4 2 2 4 2 4 4 2 2 2 2 3 4 2 4 3 3 2
##  [334] 2 3 2 2 4 2 4 4 2 2 4 3 4 4 1 3 3 3 3 2 2 3 4 4 4 2 2 4 3 3 2 2 2 2 4 3 1
##  [371] 4 2 4 1 4 4 2 4 2 4 4 3 3 3 4 2 2 4 2 2 4 4 4 4 4 4 4 4 4 2 4 4 2 4 4 2 4
##  [408] 4 2 4 3 4 2 2 3 4 4 4 3 4 2 2 4 3 4 3 2 4 2 3 4 2 3 2 2 2 2 2 3 4 4 4 1 4
##  [445] 2 2 2 3 4 3 2 3 4 2 4 1 4 4 3 3 2 4 2 3 4 2 2 4 4 3 4 4 2 2 4 4 4 2 2 3 2
##  [482] 4 2 2 4 2 2 2 4 1 4 4 4 4 4 4 4 2 4 2 4 4 3 1 4 4 2 2 4 4 2 4 2 4 2 1 4 3
##  [519] 2 3 2 4 4 2 4 3 3 3 2 2 2 4 4 3 4 1 2 2 4 2 2 4 3 2 3 2 4 4 2 2 4 4 3 3 4
##  [556] 2 2 4 4 2 2 2 4 1 4 3 4 4 4 2 3 4 4 2 4 2 2 2 2 4 1 4 2 2 4 2 3 3 2 4 2 2
##  [593] 4 4 2 2 4 4 3 2 2 4 3 4 3 4 2 4 4 2 4 2 3 4 2 2 4 2 4 3 3 2 4 4 3 2 4 3 3
##  [630] 3 4 3 3 4 4 3 4 2 4 4 4 4 1 4 4 3 2 4 4 2 3 2 3 4 3 4 4 1 2 4 4 2 4 2 2 2
##  [667] 3 4 4 2 3 2 3 4 2 2 4 2 4 3 3 4 4 4 2 2 3 2 3 4 4 3 2 3 2 2 2 3 3 2 4 4 1
##  [704] 2 3 3 4 2 2 2 4 4 2 4 2 4 4 3 4 2 4 4 4 3 2 4 4 4 3 4 1 3 3 3 3 4 4 4 4 3
##  [741] 2 4 2 2 4 4 4 4 3 2 2 3 2 4 2 4 3 3 2 4 1 4 2 2 4 2 4 3 4 4 4 4 4 4 2 2 4
##  [778] 2 4 2 3 4 3 2 3 4 2 4 4 3 4 3 1 4 4 2 3 3 4 3 2 3 4 3 2 2 3 2 2 3 4 4 2 4
##  [815] 2 4 4 3 4 3 2 4 4 4 2 2 4 2 3 4 1 3 3 4 2 4 2 4 3 2 4 2 2 4 4 2 2 4 2 3 4
##  [852] 4 2 4 2 3 4 4 4 4 2 3 3 3 2 2 3 3 2 3 4 4 4 4 4 4 2 2 4 4 3 2 2 2 4 4 4 4
##  [889] 4 4 4 4 3 4 2 4 2 4 3 4 2 1 4 4 2 4 4 3 4 2 2 3 2 4 4 2 1 4 2 2 3 4 4 2 4
##  [926] 4 2 2 4 3 2 4 2 2 4 2 2 2 2 2 4 2 2 3 4 3 3 2 4 4 4 3 2 3 4 1 3 2 4 2 4 3
##  [963] 4 4 4 4 2 3 3 2 2 3 2 4 2 4 3 3 4 2 2 3 3 2 4 4 4 4 4 4 2 4 4 4 2 3 2 2 2
## [1000] 3 2 2 2 4 2 2 4 4 4 2 4 4 4 4 3 1 3 4 2 2 2 3 4 4 2 4 3 2 2 3 2 4 2 4 4 4
## [1037] 3 2 2 4 2 4 2 2 2 4 2 3 4 4 4 4 4 4 4 3 2 4 2 4 4 4 4 2 4 4 4 4 4 2 4 2 4
## [1074] 4 2 4 2 3 2 2 4 2 3 4 1 2 2 4 4 4 2 4 2 4 2 4 3 2 2 4 2 4 2 4 4 2 4 2 4 2
## [1111] 4 4 1 2 2 4 4 2 4 4 4 4 2 4 4 4 4 1 4 3 2 2 4 4 4 4 2 2 3 2 4 4 4 2 2 2 4
## [1148] 4 3 2 4 2 3 4 4 2 2 2 4 4 4 3 3 3 2 2 4 4 4 2 2 2 2 4 2 4 2 3 3 3 2 4 4 2
## [1185] 3 3 2 3 4 3 2 4 3 2 2 4 2 3 3 4 3 4 4 2 3 4 4 2 4 2 3 2 4 2 3 4 2 2 4 3 3
## [1222] 2 4 4 4 3 4 2 2 1 2 4 2 3 4 2 2 4 4 3 2 2 2 4 4 2 2 3 4 2 3 3 4 4 4 4 3 2
## [1259] 2 3 2 4 2 2 2 3 4 2 3 2 2 1 4 3 4 3 4 4 4 2 4 1 1 3 4 4 4 3 3 2 4 4 4 1 4
## [1296] 2 3 2 4 4 4 4 4 3 4 3 2 3 4 1 2 4 3 2 4 2 2 4 2 2 2 4 2 2 2 1 3 2 2 2 4 3
## [1333] 2 4 2 4 3 2 4 4 4 4 3 3 4 3 4 2 2 4 3 2 3 2 4 3 4 2 4 4 2 4 4 2 3 2 2 4 2
## [1370] 4 4 1 4 2 4 4 2 3 2 4 4 4 4 4 4 3 4 3 4 2 2 4 4 2 4 4 2 4 2 3 3 2 4 4 2 2
## [1407] 4 4 4 4 3 4 4 4 2 4 4 2 3 2 3 4 3 4 4 4 2 3 4 2 4 4 4 4 2 3 4 4 2 4 2 3 4
## [1444] 4 4 2 4 4 2 2 3 3 3 2 4 2 2 4 2 2 2 2 4 2 3 4 2 4 4 2 4 4 4 2 3 3 2 2 2 4
## [1481] 3 3 2 1 3 4 4 3 2 2 2 2 4 2 2 2 4 3 2 3 4 2 2 4 4 3 3 3 2 4 2 4 3 3 4 3 4
## [1518] 2 4 4 2 2 2 4 2 4 4 2 4 3 2 4 2 2 2 2 2 4 4 2 2 4 3 2 2 2 4 2 3 4 4 4 3 3
## [1555] 4 2 2 2 2 2 4 4 2 4 2 4 4 2 4 3 4 4 2 4 2 2 4 2 3 4 3 3 2 2 1 4 3 4 2 2 4
## [1592] 4 2 4 2 4 4 4 2 4 2 3 4 2 3 4 3 3 2 4 4 4 2 4 4 4 4 4 4 2 3 4 2 4 2 2 2 4
## [1629] 4 4 2 4 3 4 4 1 3 3 4 2 2 4 3 2 2 2 2 2 4 2 4 4 3 2 3 4 4 4 4 4 4 4 4 2 4
## [1666] 2 4 3 2 3 2 2 2 4 4 2 4 3 4 2 3 4 4 2 2 4 1 2 4 4 4 3 2 3 4 2 4 3 2 2 4 3
## [1703] 2 3 3 4 4 4 4 2 4 3 3 2 3 4 4 4 4 2 3 4 4 3 2 2 2 2 2 4 4 2 2 4 3 2 4 4 4
## [1740] 4 4 3 2 3 4 4 4 3 2 4 2 2 3 1 3 3 2 2 2 4 4 4 4 4 2 2 2 2 4 2 4 2 4 4 3 2
## [1777] 3 4 2 4 3 2 4 4 4 2 4 2 4 4 2 4 4 2 3 3 4 2 2 4 4 2 4 4 4 4 3 4 2 3 2 2 4
## [1814] 2 3 2 4 4 2 2 4 4 3 4 4 4 4 3 1 2 4 4 2 1 3 4 4 2 2 2 2 4 4 4 2 2 2 4 2 2
## [1851] 2 2 4 4 2 2 2 4 3 4 4 2 2 3 3 2 2 4 2 4 3 4 3 3 3 2 4 4 2 2 3 4 1 2 3 2 4
## [1888] 4 4 4 4 3 2 4 4 3 3 4 4 2 4 4 2 3 4 4 2 4 2 4 4 2 2 3 4 4 2 2 2 2 4 4 2 3
## [1925] 2 2 4 2 4 3 4 4 2 4 3 2 2 4 2 4 4 2 2 3 4 3 3 4 4 4 2 4 4 3 4 2 3 4 4 2 4
## [1962] 2 4 4 4 2 2 4 4 3 3 2 4 4 4 4 4 3 4 4 3 2 4 4 1 2 2 2 2 3 4 3 2 3 4 4 4 3
## [1999] 4 3 2 2 4 3 2 4 4 4 4 4 2 4 4 2 2 4 3 4 3 1 4 4 2 2 2 4 1 2 4 4 4 3 2 2 2
## [2036] 2 3 3 2 4 4 4 2 4 2 2 4 2 4 4 4 4 3 3 4 4 2 2 3 2 4 2 3 2 4 4 2 4 3 3 4 4
## [2073] 2 2 4 3 4 2 4 4 2 2 2 4 3 3 3 4 2 4 2 4 2 4 2 2 4 2 4 2 2 3 3 2 4 4 2 3 3
## [2110] 4 1 4 2 4 3 3 1 4 4 4 4 4 4 4 3 4 4 4 2 2 2 4 3 2 4 4 4 4 2 4 2 4 3 2 4 4
## [2147] 3 4 1 4 4 3 4 4 4 3 4 3 4 4 2 4 2 2 4 2 4 2 2 2 2 3 4 3 3 2 2 4 2 3 2 3 3
## [2184] 4 3 2 2 3 4 2 3 2 4 3 2 2 4 2 2 4 2 2 4 4 4 4 2 4 4 2 2 4 4 4 2 4 2 4 3 2
## [2221] 2 3 2 3 2 2 3 3 3 4 2 2 2 2 4 4 4 4 2 2 2 2 2 4 4 3 4 2 2 4 2 4 2 2 1 2 4
## [2258] 3 2 2 4 3 4 2 1 2 4 4 4 3 2 3 1 4 4 4 2 4 4 2 3 3 4 3 2 4 1 4 3 4 3 4 4 2
## [2295] 4 2 4 2 2 3 2 2 4 3 4 4 2 3 2 3 3 3 2 2 3 2 4 2 4 3 4 2 4 2 3 2 4 2 4 2 2
## [2332] 3 3 2 3 4 4 2 4 3 4 4 2 2 2 4 2 2 3 2 4 1 3 2 3 4 4 3 4 2 4 3 4 4 4 4 2 2
## [2369] 1 1 4 3 4 4 2 3 3 2 4 4 2 1 3 4 3 4 4 4 4 3 4 4 2 4 2 4 1 2 4 2 2 4 4 3 2
## [2406] 2 4 2 4 2 4 2 2 2 4 2 3 4 4 3 4 4 4 2 2 3 4 4 4 2 4 2 4 4 4 4 3 2 3 2 4 4
## [2443] 2 4 3 4 4 2 2 3 3 3 2 2 2 4 2 3 2 2 4 3 4 4 2 4 3 2 2 4 4 1 3 4 4 3 4 1 4
## [2480] 2 4 3 2 2 2 4 2 4 4 2 4 4 3 4 2 4 2 4 2 2 3 2 2 2 2 3 3 2 4 2 2 2 4 4 2 4
## [2517] 2 2 4 4 3 4 2 1 4 2 4 2 4 2 4 1 1 2 4 3 2 2 4 2 4 4 2 4 4 4 3 3 4 2 2 2 4
## [2554] 2 4 2 3 3 2 4 4 3 3 4 3 4 4 4 3 2 2 2 3 4 4 4 4 2 2 2 3 4 4 3 2 2 2 4 3 4
## [2591] 4 4 4 2 2 4 4 4 4 4 4 2 3 3 3 4 3 3 4 2 3 2 2 4 4 4 2 4 2 2 3 4 4 4 3 2 2
## [2628] 2 2 3 3 2 3 4 4 2 2 4 2 2 3 4 3 2 3 2 4 3 2 4 2 2 4 2 2 4 3 2 2 2 4 2 3 3
## [2665] 2 4 4 4 3 3 2 2 2 2 2 4 2 2 4 4 1 2 4 2 2 1 3 4 4 1 4 4 1 2 1 4 3 4 2 2 4
## [2702] 4 3 2 3 1 2 4 2 4 4 2 1 2 2 3 2 4 2 3 2 4 3 3 4 3 3 4 3 2 2 2 4 3 4 2 4 2
## [2739] 4 4 2 3 4 4 4 2 4 4 3 2 2 2 3 3 4 4 4 2 4 3 2 2 4 2 4 4 2 4 3 4 2 4 4 4 3
## [2776] 4 3 2 2 3 4 2 3 4 3 4 4 4 4 3 2 2 2 4 3 4 1 2 4 4 4 4 2 3 3 3 2 2 3 4 4 4
## [2813] 4 2 3 4 2 4 4 4 3 4 2 4 2 2 3 4 2 3 3 4 3 3 1 4 3 4 2 2 3 3 3 2 2 4 3 4 4
## [2850] 3 4 2 2 4 4 2 3 3 3 2 2 4 2 3 2 2 2 4 2 2 2 4 4 4 2 2 4 4 3 4 1 2 4 2 2 2
## [2887] 4 2 4 3 4 4 2 1 2 2 2 2 2 3 4 3 2 4 2 2 2 2 2 3 4 4 2 2 4 4 3 4 2 2 3 4 2
## [2924] 2 3 1 4 4 3 4 4 4 4 4 4 3 2 3 2 2 4 4 2 2 4 2 4 4 3 2 4 3 4 4 2 2 4 2 4 4
## [2961] 2 2 4 2 2 2 2 2 2 2 3 4 2 3 4 4 2 2 2 3 2 2 4 2 4 2 4 3 4 3 4 4 4 4 4 3 4
## [2998] 2 2 2
## 
## Within cluster sum of squares by cluster:
## [1]  20954.07 395481.58 228290.14 302554.02
##  (between_SS / total_SS =  83.1 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Comparando estos clusters con los datos originales

table(raceCluster$cluster, Wage$race)
##    
##     1. White 2. Black 3. Asian 4. Other
##   1       68        5        7        0
##   2      864      139       57       24
##   3      482       26       59        4
##   4     1066      123       67        9

Agrupando los datos en clusters

clusplot(data.frame(wage, age), raceCluster$cluster, color=T, shade=T, lines=0, labels = 4, xlab = "Salario", ylab = "Edad")

### Notación de los clusters: - (1) - Raza Blanca - (2) - Raza Negra - (3) - Raza Asiática - (4) - Otra raza

Como podemos ver en los clusters, hay una clara diferencia en cada uno de ellos, a pesar de que hay un poco de ruido entre las razas Negra, Asiática y Otras. Lo que sí podemos observar es como la Raza Blanca esta perfectamente explicada con los salarios más altos, mientras que los salarios más bajos los tienen las personas de Raza Negra.

Ley del Codo

tot.withinss <- vector(mode="character", length=10)
for (i in 1:10){
  raceCluster <- kmeans(input, center=i, nstart=20)
  tot.withinss[i] <- raceCluster$tot.withinss
}

Graficando el codo obtenido:

plot(1:10, tot.withinss, type="b", pch=19)

Como podemos observar, la ley del codo nos dice que la cantidad de clusters óptima es de 4, por lo que podemos asegurar que nuestro modelo de k-means es óptimo.

Conclusión 2

Como pudimos observar en el segundo análisis con k-means, el salario sí puede servir para saber el tipo de raza a la que perteneces, por lo que rechazamos la hipótesis nula.

Con este análisis podemos corroborar que aún existe esta división de razas dentro del rubro laboral en Estados Unidos, pues tenemos una clara tendencia a que las peronas Blancas son las que más ingresos perciben, mientras que la raza Negra es la que menos ingresos perciben y entre medias, con un salario bajo-medio encontramos a otras razas y en un salario medio-alto a la raza Asiática.

Descargas

Descarga este código haciendo click al siguiente enlace.

xfun::embed_file("A6U.Rmd")

Download A6U.Rmd