En el Capítulo 1, en el que exploramos brevemente los dominios del aprendizaje automático, comenzamos con la regresión lineal, porque probablemente es algo con lo se ha encontrado en algún momento de su formación matemática. El proceso es bastante intuitivo y más fácil de explicar como primer concepto que algunos otros modelos de aprendizaje automático. Además, muchos ámbitos del análisis de datos se basan en modelos de regresión, desde una empresa que trata de predecir sus beneficios, hasta las fronteras de la ciencia que trata de averiguar nuevos descubrimientos que rigen las leyes del universo. Podemos encontrar la regresión en cualquier escenario en en el que se necesite una predicción contra el tiempo. En este capítulo, examinamos cómo utilizar el modelado de regresión en R en profundidad, pero también exploramos algunas advertencias y escollos a tener en cuenta en el proceso.

La motivación principal de la regresión es construir una ecuación mediante la cual podamos aprender más sobre nuestros datos. No obstante, no existe una regla rígida sobre el tipo de modelo de regresión que debe ajustarse a los datos. Elegir entre una regresión logística, una regresión lineal, o un modelo de regresión multivariada depende del problema y de los datos que se tenga. Se puede ajustar una línea recta a una serie determinada de puntos de datos, pero ¿es siempre el mejor caso? Lo ideal es encontrar un equilibrio entre la simplicidad y el poder explicativo.

Una línea recta ajustada a una serie compleja de datos puede ser sencilla, pero puede no describir el cuadro completo. Por otro lado, tener un conjunto de datos muy simple que es básicamente una línea recta y ajustar un modelo con todo tipo de curvas extravagantes puede puede proporcionar un alto grado de precisión, pero deja muy poco espacio para que se ajusten nuevos puntos de datos se ajusten a ella. Puede que recuerde en su educación matemática de la escuela secundaria sobre tener un par de puntos de datos y ajustar una línea a través de ellos. Este ajuste a los datos es la forma más sencilla de aprendizaje automático y se utiliza a menudo sin darse cuenta de que es un tipo de aprendizaje automático. Aunque el ajuste de una línea a dos puntos de datos es relativamente fácil de aprender, el ajuste de una línea con tres o más puntos de datos se convierte en una tarea más adecuada para que los computadores se encarguen de calcular ese tipo de problemas muy fácilmente. R hace que muchos de estos pasos sean bastante simples de calcular, y este capítulo proporciona una base para evaluar preguntas o cuestiones sobre dónde trazamos la línea entre la complejidad y la precisión del modelo.

Regresión lineal

En el capítulo 1, nos encontramos brevemente con la regresión lineal con un ejemplo del conjunto de datos mtcars. En ese ejemplo, se determinó una regresión lineal. En ese ejemplo, determinamos una relación lineal de la eficiencia del combustible como función del peso del vehículo y vimos que la tendencia era descendente. Extrajimos los coeficientes de una ecuación matemática lineal y nos quitamos el polvo listo. Sin embargo, hay mucho más más que simplemente poner una ecuación en un montón de datos y darlo por terminado. Volvamos a nuestro ejemplo de mtcars (Figura 4-1).

model <- lm(mtcars$mpg ~ mtcars$disp)
mpg <- predict(model); mpg
##        1        2        3        4        5        6        7        8 
## 23.00544 23.00544 25.14862 18.96635 14.76241 20.32645 14.76241 23.55360 
##        9       10       11       12       13       14       15       16 
## 23.79677 22.69220 22.69220 18.23272 18.23272 18.23272 10.14632 10.64090 
##       17       18       19       20       21       22       23       24 
## 11.46520 26.35622 26.47987 26.66946 24.64992 16.49345 17.07046 15.17456 
##       25       26       27       28       29       30       31       32 
## 13.11381 26.34386 24.64168 25.68030 15.13335 23.62366 17.19410 24.61283
b0 <- round(model$coefficients[1],2); b0
## (Intercept) 
##        29.6
b1 <- round(model$coefficients[2],2); b1
## mtcars$disp 
##       -0.04
model <- lm(mtcars$mpg ~ mtcars$disp)
plot(y = mtcars$mpg, x = mtcars$disp, xlab = "Tamaño del motor (Pulgadas cúbicas)",
 ylab = "Eficiencia del combustible (Millas por galón)", main = "Eficiencia del combustible del conjunto de datos mtcars")
abline(a = coef(model[1]), b = coef(model)[2], lty = 2)
*Figura 4-1. Una regresión lineal simple ajustada a los datos*

Figura 4-1. Una regresión lineal simple ajustada a los datos

Repasemos nuestro ejemplo de mtcars (Figura 4-1), donde modelamos la eficiencia del combustible (mpg) en función del tamaño del motor (disp) y luego observe las salidas del modelo con la función de resumen:

summary(model)
## 
## Call:
## lm(formula = mtcars$mpg ~ mtcars$disp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8922 -2.2022 -0.9631  1.6272  7.2305 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
## mtcars$disp -0.041215   0.004712  -8.747 9.38e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.251 on 30 degrees of freedom
## Multiple R-squared:  0.7183, Adjusted R-squared:  0.709 
## F-statistic: 76.51 on 1 and 30 DF,  p-value: 9.38e-10

La función summary() nos da una gran cantidad de información de este modelo lineal. Por lo general, el único número que la gente suele mirar para obtener una evaluación de precisión de referencia es el valor del R-cuadrado múltiple. Cuanto más cerca esté ese valor de \(1\), más preciso será el modelo de regresión lineal. Sin embargo, hay muchos otros términos en este resultado, así que repasemos cada elemento para obtener una comprensión sólida:

Call: Esto muestra la forma de la función que usamos. En este caso, usamos una variable de respuesta, mpg, en función de una variable dependiente, disp, estas son llamadas desde el conjunto de datos mtcars.

Residuals: Los residuos son una medida de la distancia vertical desde cada punto de datos hasta la línea ajustada en nuestro modelo. En este caso, tenemos estadísticas resumidas para todas las distancias verticales de todos nuestros puntos en relación con la línea ajustada. Cuanto menor sea este valor, mejor será el ajuste

Coefficients: Estas son las estimaciones de los coeficientes de nuestra ecuación lineal. Nuestra ecuación en este caso sería \(y = -0.04 x + 29.59\)

  • Std. Error: Con dichos coeficientes vienen las estimaciones de error dadas por Std. Parte del error de la tabla de coeficientes. En realidad, nuestra ecuación sería algo así como \(y = (- 0. 04 ± 0. 005) x + (29. 59 ± 1. 23)\)

  • t-value: Esta es la medida de la diferencia relativa a la variación en nuestros datos. Este valor está vinculado con el p-valor, los cuales se utilizan con mucha más frecuencia.

  • p-value: Los p-valores son evaluaciones estadísticas de importancia. El funcionamiento de los p-valores son más complicados que eso, pero para nuestros propósitos si el p-valor es menor que \(0.05\) significa que podemos tomar el número como estadísticamente significativo. Si el número en cuestión tiene un p-valor mayor que \(0.05\), deberíamos decir que este no es estadísticamente significativo. Los valores con asterisco se definen a continuación.

Residual standard error Esta estimación de error se refiere a la desviación estándar de nuestros datos.

Multiple R-squared: Este es el valor R cuadrado para cuando tenemos múltiples predictores. Esto no es totalmente relevante para nuestro ejemplo lineal, pero cuando agregamos más predictores al modelo, invariablemente nuestro R-cuadrado múltiple aumentará. Esto se debe a que alguna característica que agreguemos al modelo explicará alguna parte de la varianza, sea verdadera o no.

Adjusted R-squared: Para contrarrestar los sesgos introducidos por tener un valor de R cuadrado en constante aumento con más predictores, el R cuadrado ajustado tiende a ser una mejor representación de la precisión de un modelo cuando hay múltiples características.

f-statistic (prueba estadística de Fisher): Finalmente, el test-f es la razón de la varianza explicada por los parámetros en el modelo y la varianza no explicada.

Este simple ejemplo lineal tiene un poder explicativo decente. Hemos determinado una relación entre la eficiencia del combustible y el tamaño del motor. A menudo, aquí es donde los ejemplos de regresión lineal simples agotan su utilidad. Las cosas que más buscamos en este caso específico son la pendiente y la intersección. Si este ejemplo se aplicara a las ventas a lo largo del tiempo, por ejemplo, nuestro resultado de este ejercicio de modelado sería un valor inicial para la intersección y una tasa de crecimiento para el coeficiente.

Regresión multivariada

Suponga que desea construir un modelo más sólido de eficiencia de combustible con más variables integradas. La eficiencia de combustible de un vehículo puede ser un fenómeno complejo con muchos factores que afectan, además del tamaño del motor, por lo que para encontrar todas las características que son responsable del consumo de combustible debemos usar el modelo de regresión multivariada.

Recuerde que nuestro ejemplo de regresión lineal simple se basó en: \[{y} = b+ m_1x_1\] donde los coeficientes son la intersección \(b\), y \(m\) la pendiente, vinculada a la única variable que teníamos en el modelo. Si desea incorporar más factores que contribuyan al modelo, cambie la forma matemática a:

\[ {y} = b + m_{1} x_{1} + m_{2} x_{2} + m_{3} x_{3} + (...) \]

dónde \(x_1, x_2, x_3\), y así sucesivamente, son diferentes características en el modelo, como el peso del vehículo, el tamaño del motor, el número de cilindros, etc. Dado que el nuevo objetivo es encontrar coeficientes para un modelo de la forma \(y = f (x_1, x_2, x_3, ( …))\) Usemos de nuevo la función lm () de R:

lm.wt<-lm(mtcars$mpg ~ mtcars$disp + mtcars$wt )
summary ( lm.wt )
## 
## Call:
## lm(formula = mtcars$mpg ~ mtcars$disp + mtcars$wt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4087 -2.3243 -0.7683  1.7721  6.3484 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
## mtcars$disp -0.01773    0.00919  -1.929  0.06362 .  
## mtcars$wt   -3.35082    1.16413  -2.878  0.00743 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.917 on 29 degrees of freedom
## Multiple R-squared:  0.7809, Adjusted R-squared:  0.7658 
## F-statistic: 51.69 on 2 and 29 DF,  p-value: 2.744e-10

Este código amplía el modelo lineal de antes para incluir el peso del vehículo en el procedimiento de ajuste del modelo. En este caso, lo que se ve es que el R-cuadrado ajustado ha subido ligeramente de 0.709 cuando se ajusta a un modelo del tamaño del motor, a 0.7658 después de incluir el peso del motor en el ajuste. Sin embargo, observe que la relevancia estadística de la función anterior se ha reducido considerablemente. Antes, el p-valor del peso estaba muy por debajo del umbral de 0.05, que es el valor de un p-valor significativo; ahora es 0,06. Esto podría deberse a que la eficiencia de combustible del vehículo es más sensible a los cambios en el peso del vehículo que al tamaño del motor.

Si desea extender este análisis aún más, puede traer otra característica del conjunto de datos y ver cómo el R-cuadrado del modelo y los p-valores de los coeficientes cambian en consecuencia:

lm.cyl <- lm ( mtcars$mpg ~ mtcars$disp + mtcars$wt + mtcars$cyl )
summary( lm.cyl )
## 
## Call:
## lm(formula = mtcars$mpg ~ mtcars$disp + mtcars$wt + mtcars$cyl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4035 -1.4028 -0.4955  1.3387  6.0722 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 41.107678   2.842426  14.462 1.62e-14 ***
## mtcars$disp  0.007473   0.011845   0.631  0.53322    
## mtcars$wt   -3.635677   1.040138  -3.495  0.00160 ** 
## mtcars$cyl  -1.784944   0.607110  -2.940  0.00651 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.595 on 28 degrees of freedom
## Multiple R-squared:  0.8326, Adjusted R-squared:  0.8147 
## F-statistic: 46.42 on 3 and 28 DF,  p-value: 5.399e-11

Este código adopta el mismo enfoque que antes, pero agrega el recuento de cilindros del motor al modelo. Observe que el valor de R-cuadrado ha aumentado una vez más de 0.709 a 0.8147. Sin embargo, la relevancia estadística de la cilindrada (disp) en los datos es básicamente defectuosa, con un p-valor 10 veces el umbral en 0.53322 en lugar de estar más cerca de 0.05. Esto nos dice que la eficiencia del combustible está más relacionada con el conjunto de características combinadas del peso del vehículo y la cantidad de cilindros que con el tamaño del motor. Puede volver a ejecutar este análisis con solo las características estadísticamente relevantes:

lm.cyl.wt <- lm ( mtcars$mpg ~ mtcars$wt + mtcars$cyl  )
summary ( lm.cyl.wt )
## 
## Call:
## lm(formula = mtcars$mpg ~ mtcars$wt + mtcars$cyl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2893 -1.5512 -0.4684  1.5743  6.1004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.6863     1.7150  23.141  < 2e-16 ***
## mtcars$wt    -3.1910     0.7569  -4.216 0.000222 ***
## mtcars$cyl   -1.5078     0.4147  -3.636 0.001064 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.568 on 29 degrees of freedom
## Multiple R-squared:  0.8302, Adjusted R-squared:  0.8185 
## F-statistic: 70.91 on 2 and 29 DF,  p-value: 6.809e-12

Al eliminar la característica estadísticamente irrelevante del modelo, ha conservado más o menos la precisión R-cuadrado en 0.8185 versus 0.8147, mientras mantiene solo las características relevantes para los datos.

Sin embargo, debe tener cuidado al agregar funciones a los datos. En R, puede modelar fácilmente una respuesta a todas las características de los datos usando la función lm () de la siguiente forma:

lm.todos <- lm ( mpg ~ ., data =mtcars )
summary( lm.todos )
## 
## Call:
## lm(formula = mpg ~ ., data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4506 -1.6044 -0.1196  1.2193  4.6271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 12.30337   18.71788   0.657   0.5181  
## cyl         -0.11144    1.04502  -0.107   0.9161  
## disp         0.01334    0.01786   0.747   0.4635  
## hp          -0.02148    0.02177  -0.987   0.3350  
## drat         0.78711    1.63537   0.481   0.6353  
## wt          -3.71530    1.89441  -1.961   0.0633 .
## qsec         0.82104    0.73084   1.123   0.2739  
## vs           0.31776    2.10451   0.151   0.8814  
## am           2.52023    2.05665   1.225   0.2340  
## gear         0.65541    1.49326   0.439   0.6652  
## carb        -0.19942    0.82875  -0.241   0.8122  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.65 on 21 degrees of freedom
## Multiple R-squared:  0.869,  Adjusted R-squared:  0.8066 
## F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07

Esta sintaxis crea un modelo lineal con la variable dependiente mpg siendo modelada contra todo el conjunto de datos, El problema con este enfoque, sin embargo, es que ve muy poco valor estadístico en los coeficientes del modelo. Así mismo, el error estándar para cada uno de los coeficientes es muy alto, por lo que es muy difícil precisar un valor exacto para los coeficientes. En lugar de este enfoque de arriba hacia abajo para ver qué características son las más importantes en el conjunto de datos, es mejor abordarlo de abajo hacia arriba como lo hemos hecho hasta ahora. Aunque el tema de la selección de funciones en sí es un tema muy amplio, (uno que exploramos en profundidad con otros algoritmos de aprendizaje automático) podemos mitigar algunos de estos problemas de dos maneras:

Selección cuidadosa de funciones:
Elija características para agregar al modelo una a la vez y elimine las que sean estadísticamente insignificantes. Hemos logrado esto en los fragmentos de código anteriores agregando un parámetro a la vez y verificando si el p-valor de la salida del modelo para ese parámetro es estadísticamente significativo.

Regularización:
Mantenga todas las características pero reduzca matemáticamente los coeficientes de las menos importantes para minimizar su impacto en el modelo.

Regularización

Regularización puede ser un concepto difícil matemáticamente, pero en principio es bastante sencillo. La idea es que desee incluir tantas características en sus datos como pueda incluir en el modelo final. Cuantas más funciones, mejor podrá explicar todas las complejidades del conjunto de datos. El problema aquí es que el grado en que cada característica explica parte del modelo puede ser bastante diferente después de que se aplica la regularización.

Mediante el uso de la regularización, puede hacer que su modelo sea más conciso y reducir el ruido en el conjunto de datos que podría provenir de características que tienen poco impacto en lo que está tratando de modelar.

Veamos cuál es el modelo lineal para el mtcars el conjunto de datos se vería si incluyéramos todas las características. Tendríamos una ecuación como esta:

\[ {mpg} = 12.3 - 0.11_{cyl} + 0.01_{disp} - 0.02_{hp} + 0.79_{drat} - 3.72_{wt} + 0.82_{qsec} + 0.31_{vs} + 2.42_{am} + 0.66_{gear} - 0.20_{carb}\] Según esta ecuación lineal, la eficiencia del combustible es más sensible al peso del vehículo \((-3.72_{wt})\), dado que este tiene el coeficiente más grande. Sin embargo, la mayoría de estos están todos dentro de un orden de magnitud entre sí. La regularización mantendría todas las características, pero las menos importantes tendrían sus coeficientes reducidos mucho más

Para utilizar esta técnica de regularización, usamos un tipo particular de modelo de regresión, conocido como regresión lasso, como se muestra aquí:

require(lasso2)
lm.lasso<-l1ce( mpg ~ ., data = mtcars ) 
summary(lm.lasso)$coefficients
##                   Value  Std. Error     Z score   Pr(>|Z|)
## (Intercept) 36.01809203 18.92587647  1.90311355 0.05702573
## cyl         -0.86225790  1.12177221 -0.76865686 0.44209704
## disp         0.00000000  0.01912781  0.00000000 1.00000000
## hp          -0.01399880  0.02384398 -0.58709992 0.55713660
## drat         0.05501092  1.78394922  0.03083659 0.97539986
## wt          -2.68868427  2.05683876 -1.30719254 0.19114733
## qsec         0.00000000  0.75361628  0.00000000 1.00000000
## vs           0.00000000  2.31605743  0.00000000 1.00000000
## am           0.44530641  2.14959278  0.20715850 0.83588608
## gear         0.00000000  1.62955841  0.00000000 1.00000000
## carb        -0.09506985  0.91237207 -0.10420075 0.91701004

Este código usa función l1ce() del paquete lasso2 en el conjunto de datos mtcars. Este usa la misma función que queremos que modele la variable de eficiencia de combustible mpg como una función de todas las otras variables en el conjunto de datos. Integrado en la regresión lasso está la técnica de regularización que solo se aplica durante la parte matemática pesada del algoritmo. La parte de regularización de la regresión escala los coeficientes de acuerdo con el impacto real que tienen en el modelo de una manera más estadística. En algunos casos, esto puede provocar que algunas características se reduzcan a un valor tan bajo que se aproximen a cero. Como resultado de este modelo de regresión, ahora tiene una ecuación diferente:

\[ {mpg} = 36.02 - 0.86_{cyl} + 0_{disp} - 0.014_{hp} + 0.06_{drat} - 2.69_{wt} + 0_{qsec} + 0_{vs} + 0.45_{am} + 0_{gear} - 0.095_{carb}\] O, más simplemente:

\[ {mpg} = 36.02 - 0.86_{cyl} - 0.014_{hp} + 0.06_{drat} - 2.69_{wt} + 0.45_{am} - 0.095_{carb}\] La característica más importante antes del cambio a una regresión de lasso era el peso del vehículo wt, peso que no ha variado en cuanto a su importancia relativa. Aunque el coeficiente ha cambiado algo, el hecho de que sea el coeficiente de mayor magnitud sigue siendo el mismo. Lo que ve en términos de características menos útiles que se reducen, en este caso a cero, son características que probablemente pensaría que tienen poco impacto en la eficiencia del combustible para empezar. Tiempo de carrera de un cuarto de milla (qsec), configuración del motor en términos de una forma V o en línea recta (vs), y número de marchas hacia adelante (gear), han sido todos reducidos a cero.

Sin embargo, la variable de desplazamiento mostró una clara relación con la eficiencia del combustible que vimos anteriormente. Que se reduzca a cero no significa que no haya relación entre esa única variable y nuestra respuesta, pero cuando se toma junto con todas las demás variables del conjunto de datos, su importancia es insignificante.

Recuerde, en este caso nos interesa un modelo de todas las características, no necesariamente en la importancia de una sola característica.

Observe en el nuevo modelo de regresión de lasso que algunos de los coeficientes se han eliminado más o menos matemáticamente del modelo. Para refinar aún más el modelo y reducir la cantidad de características en él, puede volver a ejecutar la regresión sin esas características y ver qué cambios:

require(lasso2)
lm.lasso2 <- l1ce(mpg ~ cyl + hp + wt + am + carb , data = mtcars )
summary(lm.lasso2)$coefficients
##                     Value Std. Error     Z score     Pr(>|Z|)
## (Intercept) 31.2819166926 4.51160542  6.93365527 4.100942e-12
## cyl         -0.7864202230 0.86107128 -0.91330444 3.610824e-01
## hp          -0.0009037003 0.02343634 -0.03855979 9.692414e-01
## wt          -1.9248597501 1.38749433 -1.38729198 1.653527e-01
## am           0.0000000000 2.22143917  0.00000000 1.000000e+00
## carb         0.0000000000 0.67947216  0.00000000 1.000000e+00

Con el conjunto de datos reducido que luego se pasa a otra regresión de lasso, puede ver que él tipo de transmisión del automóvil, am, y el número de carburadores, carb ambos se han reducido a cero. Al eliminar estas funciones y volver a ejecutarlas, puede ver si se bajar más:

lm.lasso3 <- l1ce( mpg ~ cyl + hp + wt , data = mtcars )
summary(lm.lasso3)$coefficients
##                  Value Std. Error    Z score  Pr(>|Z|)
## (Intercept) 30.2106931 1.97117597 15.3262284 0.0000000
## cyl         -0.7220771 0.82941877 -0.8705821 0.3839824
## hp           0.0000000 0.01748364  0.0000000 1.0000000
## wt          -1.7568469 1.07478525 -1.6346028 0.1021324

En este caso, los caballos de fuerza del automóvil, hp, se ha retirado. Puede continuar ejecutándose siempre que tenga varias funciones para probar:

lm.lasso4 <- l1ce ( mpg ~ cyl + wt , data = mtcars )
summary( lm.lasso4 )$coefficients
##                  Value Std. Error   Z score  Pr(>|Z|)
## (Intercept) 29.8694933  1.4029760 21.290096 0.0000000
## cyl         -0.6937847  0.5873288 -1.181254 0.2375017
## wt          -1.7052064  1.0720172 -1.590652 0.1116879

El resultado final es un modelo que tiene solo dos características en lugar de las 11 con las que comenzó:

\[ {mpg} = 29.87-0.69_{cyl}-1.70_{wt} \]

Regresión polinómica

La regresión polinómica consiste simplemente en ajustar una función de grado superior a los datos. Anteriormente, hemos visto ajustes a nuestros datos de la siguiente forma

\[ {y} = b + m_{1} x_{1} + m_{2} x_{2} + m_{3} x_{3} + (...) \] La regresión polinómica se diferencia de los casos lineales simples por tener múltiples grados para cada característica del conjunto de datos. Su forma podría representarse de la siguiente manera

\[ {y} = b + m_{1} x_{1}^2 \]

El siguiente ejemplo ayudará a razonar mejor (Figura 4-2):

pop<-data.frame(uspop) 
pop$uspop <- as.numeric(pop$uspop)
pop$year <- seq(from = 1790, to = 1970, by = 10)
plot(y = pop$uspop, x = pop$year, main = "Población de Estados Unidos de 1790 a 1970",xlab = "Año", ylab = "Población")
*Figura 4-2. La población trazada de los Estados Unidos en las décadas de 1790 a 1970*

Figura 4-2. La población trazada de los Estados Unidos en las décadas de 1790 a 1970

Aquí tenemos un conjunto de datos incorporado en R que hemos modificado ligeramente para fines de demostración. Normalmente, uspop es un objeto de serie de tiempo que tiene sus propios criterios de trazado, pero aquí lo hemos ajustado para trazar sólo los puntos de datos. Estos datos son la población de los Estados Unidos en períodos de 10 años desde 1790 hasta 1970. Comencemos por ajustar un modelo lineal a los datos:

lm1<- lm(pop$uspop ~ pop$year) 
summary(lm1)
## 
## Call:
## lm(formula = pop$uspop ~ pop$year)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.569 -14.776  -2.933   9.501  36.345 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.958e+03  1.428e+02  -13.71 1.27e-10 ***
## pop$year     1.079e+00  7.592e-02   14.21 7.29e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.12 on 17 degrees of freedom
## Multiple R-squared:  0.9223, Adjusted R-squared:  0.9178 
## F-statistic: 201.9 on 1 and 17 DF,  p-value: 7.286e-11

Este ajuste lineal simple de los datos parece funcionar bastante bien. Los valores p de las estimaciones son muy bajos, lo que indica una buena significación estadística. Asimismo, los valores de R-cuadrado son muy buenos. Sin embargo, los residuos muestran un grado de variabilidad bastante amplio, que llega hasta una diferencia de 36, como se demuestra en la siguiente figura (Figura 4-3):

plot(y = pop$uspop, x = pop$year, main = "Población de Estados Unidos de 1790 a 1970",
xlab = "Año", ylab = "Población")
abline(a = coef(lm1[1]), b = coef(lm1)[2], lty = 2, col = "red")
*Figura 4-3. Datos de población con un ajuste de modelo lineal*

Figura 4-3. Datos de población con un ajuste de modelo lineal

El ajuste de la línea de puntos del modelo lineal parece estar bien. Se ajusta a algunos de los datos mejor que a otros, pero está bastante claro por los datos que no es exactamente una relación lineal. Además, sabemos por intuición que la población a lo largo del tiempo tiende a tener una forma más exponencial que una línea recta. Lo que se quiere hacer a continuación es ver cómo se compara un modelo de mayor grado con el caso lineal, que es el polinomio de menor grado que se puede ajustar:

lm2 <- lm(pop$uspop ~ poly(pop$year, 2))
summary(lm2)
## 
## Call:
## lm(formula = pop$uspop ~ poly(pop$year, 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5997 -0.7105  0.2669  1.4065  3.9879 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         69.7695     0.6377  109.40  < 2e-16 ***
## poly(pop$year, 2)1 257.5420     2.7798   92.65  < 2e-16 ***
## poly(pop$year, 2)2  73.8974     2.7798   26.58 1.14e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.78 on 16 degrees of freedom
## Multiple R-squared:  0.9983, Adjusted R-squared:  0.9981 
## F-statistic:  4645 on 2 and 16 DF,  p-value: < 2.2e-16
lm3 <- lm(pop$uspop ~ poly(pop$year, 3)) 
summary(lm3)
## 
## Call:
## lm(formula = pop$uspop ~ poly(pop$year, 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2615 -0.7823  0.4491  1.4350  3.6180 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         69.7695     0.6504 107.273  < 2e-16 ***
## poly(pop$year, 3)1 257.5420     2.8350  90.844  < 2e-16 ***
## poly(pop$year, 3)2  73.8974     2.8350  26.066  6.6e-14 ***
## poly(pop$year, 3)3   1.7543     2.8350   0.619    0.545    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.835 on 15 degrees of freedom
## Multiple R-squared:  0.9983, Adjusted R-squared:  0.998 
## F-statistic:  2977 on 3 and 15 DF,  p-value: < 2.2e-16
lm4 <- lm(pop$uspop ~ poly(pop$year, 4))
summary(lm4)
## 
## Call:
## lm(formula = pop$uspop ~ poly(pop$year, 4))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2223 -1.4251  0.5952  1.6344  3.5346 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         69.7695     0.6387 109.234  < 2e-16 ***
## poly(pop$year, 4)1 257.5420     2.7841  92.505  < 2e-16 ***
## poly(pop$year, 4)2  73.8974     2.7841  26.543 2.25e-13 ***
## poly(pop$year, 4)3   1.7543     2.7841   0.630    0.539    
## poly(pop$year, 4)4   3.4701     2.7841   1.246    0.233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.784 on 14 degrees of freedom
## Multiple R-squared:  0.9985, Adjusted R-squared:  0.9981 
## F-statistic:  2316 on 4 and 14 DF,  p-value: < 2.2e-16
lm5 <- lm(pop$uspop ~ poly(pop$year, 5))
summary(lm5)
## 
## Call:
## lm(formula = pop$uspop ~ poly(pop$year, 5))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8151 -0.9071  0.0811  0.7011  3.4618 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         69.7695     0.4763 146.480  < 2e-16 ***
## poly(pop$year, 5)1 257.5420     2.0762 124.046  < 2e-16 ***
## poly(pop$year, 5)2  73.8974     2.0762  35.593 2.41e-14 ***
## poly(pop$year, 5)3   1.7543     2.0762   0.845    0.413    
## poly(pop$year, 5)4   3.4701     2.0762   1.671    0.119    
## poly(pop$year, 5)5   7.2442     2.0762   3.489    0.004 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.076 on 13 degrees of freedom
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9989 
## F-statistic:  3334 on 5 and 13 DF,  p-value: < 2.2e-16
lm6 <- lm(pop$uspop ~ poly(pop$year, 6))
summary(lm6)
## 
## Call:
## lm(formula = pop$uspop ~ poly(pop$year, 6))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9298 -0.6066 -0.0685  0.3543  4.0751 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          69.769      0.469 148.770  < 2e-16 ***
## poly(pop$year, 6)1  257.542      2.044 125.986  < 2e-16 ***
## poly(pop$year, 6)2   73.897      2.044  36.149 1.29e-13 ***
## poly(pop$year, 6)3    1.754      2.044   0.858  0.40761    
## poly(pop$year, 6)4    3.470      2.044   1.698  0.11535    
## poly(pop$year, 6)5    7.244      2.044   3.544  0.00404 ** 
## poly(pop$year, 6)6    2.427      2.044   1.187  0.25808    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.044 on 12 degrees of freedom
## Multiple R-squared:  0.9993, Adjusted R-squared:  0.999 
## F-statistic:  2866 on 6 and 12 DF,  p-value: < 2.2e-16

Este código usa de nuevo la función lm(), pero esta vez con un parámetro extra alrededor de la variable dependiente, la función poly(). Esta función toma los datos y calcula un vector ortogonal, que luego se escala adecuadamente. Por defecto, la función poly() no cambia los valores de los datos, pero puede utilizarla para ver si produce mejores resultados que el ajuste de orden inferior que hizo anteriormente. Recuerde que el ajuste lineal es técnicamente un polinomio, pero de grado \(1\). En una ecuación, aquí está el ajuste del modelo resultante: \[y=b+m_1x_1^ 2 +m_2 x_1 \]

Vamos a caminar lentamente a través de la salida de summary() primero. Mirar la salida de los residuos nos da un poco de alivio: ¡no hay residuos en el rango de 30!. Los residuos más pequeños son siempre mejores en términos de ajuste del modelo. La tabla de coeficientes tiene ahora tres entradas: una para el intercepto, una para el término de primer grado, y ahora una para el término de segundo grado. Cuando usó poly(pop$year, 2), le indicó a R que quería un polinomio de los datos con el grado más alto siendo 2. Volviendo a la tabla de coeficientes, puede ver que todos los p-valores son estadísticamente significativos, lo que también es una buena indicación de que este es un modelo sólido que se ajusta a sus datos (ver Figura 4-4):

plot(y =pop$uspop, x = pop$year, main = "Población de Estados Unidos de 1790 a 1970",xlab ="Año",ylab = "Población") 
pop$lm2.predict = predict(lm2, newdata = pop);lines(sort(pop$year),fitted(lm2)[order(pop$year)], col = "blue",lty = 2)
*Figura 4-4. La población a lo largo del tiempo modelada con un ajuste cuadrático parece ajustarse a los datos mucho mejor que uno lineal; Sin embargo, si desea el modelo más preciso posible, es posible que deba aumentar el grado de polinomio en el que ajusta los datos*

Figura 4-4. La población a lo largo del tiempo modelada con un ajuste cuadrático parece ajustarse a los datos mucho mejor que uno lineal; Sin embargo, si desea el modelo más preciso posible, es posible que deba aumentar el grado de polinomio en el que ajusta los datos

Figura 4-5. Población a lo largo del tiempo con ajuste de múltiples modelos

De la (Figura 4-4), parece bastante obvio que el polinomio de mayor grado (en este caso una ecuación cuadrática) se ajusta mejor a los datos. Es evidente que el uso de polinomios de mayor grado funciona mejor que los de menor grado, ¿verdad? ¿Qué ocurre si se ajusta un polinomio de tercer grado? ¿O algo más alto aún? Apuesto a que si se utiliza un polinomio de sexto grado tendrías un modelo muy preciso. Lo que salta inmediatamente es que el ajuste lineal simple que tenías antes se ajustaba a los datos lo mejor posible, pero el polinomio de segundo grado o más alto (es decir, un cuadrático simple) se ajusta mejor. Una mejor manera de distinguir la diferencia entre los ajustes polinómicos de orden superior es observar los gráficos de los residuos de cada modelo, que puede ver en la (Figura 4-5).

La figura 4-5 muestra el ajuste lineal en comparación con polinomios de grado creciente. Los polinomios son difíciles de separar en términos de qué tan bien encajan, pero todos parecen encajar mejor que el caso lineal. Para comparar modelos que son aproximaciones cercanas visualmente en este nivel, necesitaría sumergirse en mirar gráficos de sus residuos, como se muestra en la Figura 4-6:

par(mfrow = c(2, 3))
plot(resid(lm1),main ="Grado 1",xlab = "Sucesión de años", 
plot(resid(lm2),main = "Grado 2",xlab = "Sucesión de años", 
plot(resid(lm3),main = "Grado 3",xlab = "Sucesión de años",
plot(resid(lm4),main = "Grado 4",xlab="Sucesión de años",
plot(resid(lm5),main="Grado 5",xlab="Sucesión de años", 
plot(resid(lm6), main="Grado 6", xlab="Sucesión de años",ylab ="Residuos"),ylab="Residuos"),ylab="Residuos"),ylab="Residuos"),ylab = "Residuos"),ylab="Residuos")
*Figura 4-6. Una gráfica de residuos de cada uno de los modelos*

Figura 4-6. Una gráfica de residuos de cada uno de los modelos

Recordemos que el residuo es la distancia vertical entre un punto de datos y la línea del modelo ajustado. Un modelo que se ajusta exactamente a los puntos de datos debería tener un gráfico de residuos tan cercano a una línea plana como sea posible. En el caso de tu ajuste lineal, la escala del gráfico de residuos es mucho más grande que el resto, y puede ver que el ajuste lineal tiene una distancia residual bastante mala en su inicio, punto medio y final. Este no es un modelo ideal. Por otro lado, los polinomios de mayor grado parecen hacerlo bastante bien. La escala de sus gráficos de residuos es mucho más agradable, pero el que realmente destaca es el de sexto grado. El gráfico de residuos es prácticamente nulo al principio, y luego se vuelve un poco más propenso a los errores.

Todo esto está muy bien, pero podría ser más fácil clasificar el ajuste del modelo mirando sus residuos numéricamente:

c(sum(abs(resid(lm1))), sum(abs(resid(lm2))), sum(abs(resid(lm3))), sum(abs(resid(lm4))), sum(abs(resid(lm5))), sum(abs(resid(lm6))))
## [1] 272.51432  33.77224  34.54039  36.95125  25.45242  19.59938

Este código suma el valor absoluto del residuo. Si sólo se toma la suma bruta de los residuos, se obtiene una imagen inexacta porque algunos residuos pueden ser negativos. Así, el residuo total para el ajuste lineal es cuantitativamente malo comparado con el resto de los modelos, siendo el polinomio de sexto grado el claro ganador en términos de mejor ajuste a los puntos de datos.

Pero, ¿el mejor ajuste a los puntos de datos es realmente el mejor modelo?. Debemos tener en cuenta los conceptos de sobreajuste y subajuste de los datos. El modelo lineal ajustado a los datos en el caso anterior sería un buen ejemplo de un escenario de subajuste. Está claro que hay alguna estructura en los datos que no se explica con un simple ajuste lineal. Por otro lado, un modelo puede estar sobreajustado si es demasiado específico para los datos presentados y ofrece poco poder explicativo para cualquier dato nuevo que pueda entrar en el sistema. Este es el riesgo que se corre al aumentar el grado de los modelos polinómicos.

Bondad del ajuste con datos: los peligros del sobreajuste

Acabamos de ejecutar un ejemplo para tratar de obtener un modelo que tenga el mejor ajuste a nuestros datos. Este es un buen objetivo, pero hay que tener cuidado de no ir demasiado lejos en el ajuste perfecto a los datos. Hemos visto hasta ahora que el ajuste lineal a una curva de población probablemente no es el mejor modelo para el trabajo. Un ajuste cuadrático o un polinomio cúbico parecen hacerlo mucho mejor en comparación. Sin embargo, ¿merece la pena seguir aumentando el grado de ajuste del modelo? ¿Es la minimización del residuo el único objetivo a la hora de seleccionar el mejor modelo para el trabajo?

Error cuadrático medio

En estadística, el error cuadrático medio (RMSE) es una forma cuantificable de ver cómo se comparan los valores predichos de nuestro modelo con los valores reales. Matemáticamente, el RMSE se da como sigue:

\[RMSE =\sqrt{( \text{predictedvalue} − \text{actualvalue})^ 2}\]

Para evaluar los ajustes polinómicos, puede realizar un análisis de RMSE en cada uno de ellos. A continuación, puede comparar los errores resultantes y seleccionar el que tenga el resultado más bajo. Para ello, necesitas nuevos datos que no estén en tu modelo. Para ello, vamos a utilizar los datos del censo de población de EE.UU. de 1980 a 2010:

uspop.2020<-data.frame(year = c(1980, 1990, 2000, 2010), uspop=c(226.5,249.6, 282.2, 309.3))
uspop.2020.predict <- uspop.2020
pop2 <- data.frame(uspop)
pop2$uspop <- as.numeric(pop$uspop)
pop2$year <- seq(from = 1790, to = 1970, by = 10)

Este código también reinicializa los datos de población antiguos con fines de predicción como medida de limpieza general. A partir de ahí, puedes hacer tu rutina de predicción habitual y luego calcular el RMSE para cada polinomio:

uspop.2020.predict$lm1 <- predict(lm(uspop ~ poly(year, 1), data = pop2), uspop.2020)
uspop.2020.predict$lm2 <- predict(lm(uspop ~ poly(year, 2), data = pop2), uspop.2020)
uspop.2020.predict$lm3 <- predict(lm(uspop ~ poly(year, 3), data = pop2), uspop.2020)
uspop.2020.predict$lm4 <- predict(lm(uspop ~ poly(year, 4), data = pop2), uspop.2020)
uspop.2020.predict$lm5 <- predict(lm(uspop ~ poly(year, 5), data = pop2), uspop.2020)
uspop.2020.predict$lm6 <- predict(lm(uspop ~ poly(year, 6), data = pop2), uspop.2020)
#And, finally, calculate the RMSE:
c(sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm1)^2)),
sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm2)^2)), sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm3)^2)), sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm4)^2)), sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm5)^2)), sqrt(mean((uspop.2020.predict$uspop - uspop.2020.predict$lm6)^2)))
## [1]  75.622445   8.192311   5.070814   9.153189  73.632318 124.429798

A partir de estos resultados, puede ver que el ajuste lineal simple tuvo un RMSE de 75, el polinomio de segundo grado tuvo 8, y el polinomio de tercer grado tuvo 5. Los errores se disparan después del polinomio de tercer grado, lo que es otra indicación de que los modelos se ajustaron demasiado a los datos. En este caso, se seleccionaría el modelo que tiene el menor RMSE para los nuevos datos predichos, eligiendo el polinomio de grado tres.

Simplicidad del modelo y bondad del ajuste

Si recuerda los coeficientes del modelo, cada uno de ellos tiene un p-valor de significación estadística vinculado al procedimiento de ajuste del modelo lm(). Si el p-valor de un coeficiente es inferior a 0.05, es seguro asumir que es estadísticamente importante para su modelo.

Para ayudarle a decidir qué modelo utilizar, identifique dónde está el equilibrio entre la precisión y la complejidad del modelo. Cuanto más complejo sea el modelo, es decir, el grado del polinomio utilizado, más se ajustará a los datos, pero se corre el riesgo de que algunos de los coeficientes sean menos válidos estadísticamente a medida que el modelo se hace más complejo. Para evitarlo, observe primero el R-cuadrado y el número de coeficientes estadísticamente válidos de cada uno de sus modelos:

table((summary(lm1)$coefficients[, 4]) < 0.05)
## 
## TRUE 
##    2
summary(lm1)$r.squared
## [1] 0.9223434

Este ejemplo toma los coeficientes del ajuste lineal simple, lm1, y luego extrae los p-valores vinculados a los coeficientes. A continuación, tabula cuántos de ellos son estadísticamente válidos (si son superiores a 0.05). El resultado del caso lineal simple es que hay dos coeficientes: la pendiente y el intercepto, y que ambos son estadísticamente válidos. El valor de R-cuadrado también confirma que el ajuste es bastante bueno, pero vamos a utilizarlo como referencia para comparar.

En lugar de calcular esto para cada modelo y mirar hacia adelante y hacia atrás en los resultados, puede volcar toda esta información en un marco de datos práctico para facilitar la lectura. Vamos a definir model.order como el grado más alto del ajuste polinómico (esto es simplemente el número que se pasa a la función poly() durante el uso de la función lm() del modelo lineal. A continuación, defina coef.true como el número de coeficientes que son estadísticamente válidos en el modelo. En este caso, usted está mirando sólo los coeficientes relacionados con las variables dependientes y no el intercepto del modelo, que es estadísticamente válido en todos los casos, de ahí que reste el valor de coef.true por 1. A continuación, define un término coef.false como el caso opuesto: cuántos de los coeficientes del modelo sobre las variables dependientes no son estadísticamente significativos. Por último, se define un valor model.rsq, que es la exactitud extraída del modelo R-cuadrado. A continuación, se reúne todo en un marco de datos y se define una última métrica: la bondad. Esta medida compara la relación entre los coeficientes estadísticamente significativos y el orden del modelo:

model.order <- c(1,2,3,4,5,6)
coef.true <- c(
table((summary(lm1)$coefficients[,4])<0.05) - 1
,table((summary(lm2)$coefficients[,4])<0.05) - 1
,table((summary(lm3)$coefficients[,4])<0.05)[2] - 1
,table((summary(lm4)$coefficients[,4])<0.05)[2] - 1
,table((summary(lm5)$coefficients[,4])<0.05)[2] - 1
,table((summary(lm6)$coefficients[,4])<0.05)[2] - 1
)
coef.false <- c(
0
,0
,table((summary(lm3)$coefficients[,4])<0.05)[1]
,table((summary(lm4)$coefficients[,4])<0.05)[1]
,table((summary(lm5)$coefficients[,4])<0.05)[1]
,table((summary(lm6)$coefficients[,4])<0.05)[1]
)
model.rsq <- c(
summary(lm1)$r.squared
,summary(lm2)$r.squared
,summary(lm3)$r.squared
,summary(lm4)$r.squared
,summary(lm5)$r.squared
,summary(lm6)$r.squared
)
model.comparison <- data.frame(model.order, model.rsq, coef.true, coef.false)
model.comparison$goodness <- (model.comparison$coef.true / model.comparison
$model.order)
knitr::kable(model.comparison)
model.order model.rsq coef.true coef.false goodness
1 0.9223434 1 0 1.0000000
2 0.9982808 2 0 1.0000000
3 0.9983235 2 1 0.6666667
4 0.9984910 2 2 0.5000000
5 0.9992208 3 2 0.6000000
6 0.9993027 3 3 0.5000000

El resultado demuestra que, aunque la precisión del R-cuadrado del modelo puede aumentar a medida que el ajuste se vuelve más complejo, la bondad de ese ajuste disminuye con el tiempo porque el número de coeficientes estadísticamente significativos en comparación con el número total de coeficientes tiende a disminuir. Una forma de cuantificar estadísticamente esto es clasificar los elementos asociados que le interesa optimizar con lo siguiente:

model.comparison$rank<-sqrt(model.comparison$goodness^2 + model.comparison$model.rsq^2) 
knitr::kable(model.comparison)
model.order model.rsq coef.true coef.false goodness rank
1 0.9223434 1 0 1.0000000 1.360411
2 0.9982808 2 0 1.0000000 1.412998
3 0.9983235 2 1 0.6666667 1.200456
4 0.9984910 2 2 0.5000000 1.116685
5 0.9992208 3 2 0.6000000 1.165522
6 0.9993027 3 3 0.5000000 1.117410

Lo que podemos decir de este procedimiento es que tenemos un modelo óptimo a elegir que tiene el valor de rango más alto. Un modelo que tiene un R-cuadrado más bajo y un rango más bajo no se ajusta a los datos. Un modelo que tiene un R-cuadrado más alto y un rango más bajo es un modelo sobreajustado a nuestros datos.

Regresión Logística

Hasta ahora hemos considerado los modelos de regresión en términos de tomar algún tipo de datos numéricos a los que queremos ajustar algún tipo de curva para poder utilizarla con de predicción. La regresión lineal toma algún tipo de datos numéricos y hace una ecuación dela forma \[y = mx + b\] La regresión lineal también puede tener múltiples entradas y podríamos tener una ecuación de la forma: \[y= b+m_1x_1+ m_2x_2 +(...).\] Además, estos tipos de modelos de de regresión pueden convertirse en casos no lineales de la forma \[y = b + m_1x_1 + m_2x_1^2 +m_3x_1^3 + (…)\] Todos ellos tienen sus propios casos de uso y dependen totalmente de los datos con los que trabajemos y de la estrategia que sigamos sobre el tipo de precisión que queremos optimizar. Hasta ahora, todos estos métodos han incorporado una entrada numérica y nos han proporcionado una salida numérica. ¿Y si, en cambio, quisiéramos un resultado “sí” o “no” de nuestros datos? ¿Y si tratáramos de hacer algo como determinar si nuestros datos de entrada tuvieran un resultado positivo o negativo?.
En este caso, estaríamos tomando datos numéricos continuos y obteniendo algún tipo de resultado discreto. Esta es la base para el extremo de clasificación de nuestro modelado de regresión.
La regresión logística es un tipo particular de clasificación y relativamente simple para ser utilizada como un ejemplo introductorio. La regresión logística, a diferencia de la regresión lineal, encuentra el punto en el que los datos pasan de un tipo de clasificación a otro, en lugar de tratar de ajustar todos los puntos de datos individuales en sí mismos.

La motivación de la clasificación

Supongamos que se intenta diagnosticar a los pacientes para determinar si tienen un tumor maligno. Veamos el código y el gráfico resultante en la Figura 4-7:
*Figura 4-7. Ajustar una línea de regresión lineal a datos binarios no proporciona un modelo preciso*

Figura 4-7. Ajustar una línea de regresión lineal a datos binarios no proporciona un modelo preciso 

## (Intercept)  tumor.size 
##  0.20380952  0.06095238
## [1] 0.4063492

Este código crea un conjunto de datos de tamaños de tumores de 1 a 20 y clasifica si son malignos, clasificando un tumor benigno o no canceroso como 0, y un tumor maligno o canceroso se clasifica como 1. Un instinto ingenuo podría ser el de aplicar un modelo de regresión sobre estos datos y ver cuál es el resultado de la ecuación. Con este enfoque, usted obtendría una ecuación como la siguiente \[\textrm{Malignidad del tumor} = 0,204+0,061 ×\textrm{Tamaño del tumor} \] El mal ajuste de la R-cuadrado en 0,406 sugiere que podríamos obtener un modelo mas preciso. Además, debería cuestionar la evaluación lógica de lo que significa tener un tumor que es 0,2 maligno cuando se registran en los datos como siendo malignos o no, sin espacio intermedio.
Esto tampoco tendría mucho sentido con el conjunto de datos mtcars si tuviéramos algo modelado contra el tipo de transmisión. ¿Qué sería un 0,2 de transmisión si el 0 fuera manual y el 1 fuera automático? Tenemos que replantearnos este enfoque. En lugar de ajustar una función matemática normal necesitamos ajustar algo llamado límite de decisión a los datos.

El límite de decisión

El límite de decisión es simplemente una línea en la arena de nuestros datos que dice “cualquier cosa en este lado se clasifica como \(X\) y cualquier cosa del otro lado se clasifica como \(Y\)”. La figura 4-8 revisa el gráfico de los tamaños de los tumores y si son malignos. Puede ver claramente que cualquier tumor que tenga un tamaño superior a 5 siempre parece ser maligno:

*Figura 4-8.Trazar un límite de decisión en lugar de una línea de regresión clasifica los datos menores de aproximadamente 4.5 como 0; los datos por encima de ese umbral se clasifican como 1*

Figura 4-8.Trazar un límite de decisión en lugar de una línea de regresión clasifica los datos menores de aproximadamente 4.5 como 0; los datos por encima de ese umbral se clasifican como 1

La regresión logística establece el límite con el que se pueden clasificar los datos. El límite en la Figura 4-8 muestra que cualquier tamaño de tumor superior a 4,5 es maligno mientras que cualquier tamaño inferior es benigno.

La función sigmoide

La forma en que la regresión logística (así como muchos otros tipos de algoritmos de clasificación) se basa en los fundamentos matemáticos de la función sigmoide. La función sigmoide tiene la siguiente forma matemática:

\[ h(x)=\frac{1}{1+e^{-x}}\] La Figura 4-9 muestra cómo se ve la gráfica:

  e<-exp(1)
curve(1/(1 + e^-x), -10, 10, main = "Función sigmoide", xlab = "Entrada",
      ylab = "Probabilidad")
*Figura 4-9. La función sigmoidea es la base de la regresión logística*

Figura 4-9. La función sigmoidea es la base de la regresión logística

Esta función se utiliza en la regresión logística para clasificar los datos. Por sí misma, la función toma algún valor numérico que nos interesa y lo asigna a una probabilidad entre 0 y 1. Podríamos tener la tentación de introducir algunos de los valores de nuestro ejemplo anterior en la función sigmoide y ver cuál es el resultado.
Si lo hiciéramos, como poner \(x = 1\), por ejemplo, obtendríamos \(h(1) = 0,73\), o aproximadamente un 73% de probabilidad de que un tumor sea maligno si nuestra entrada es 1. Sin embargo, nuestro sistema de clasificación es 0 para benigno y 1 para maligno. La entrada de longitud 1 da un resultado de 0,73, que es incorrecto.

  e<-exp(1)
x=1
h<-(1/(1 + e^-x))
h
## [1] 0.7310586

En su lugar, necesitamos pasar un conjunto de parámetros ponderados al regresor logístico. Dado que por el momento sólo tenemos una variable dependiente (teniendo en cuenta que el eje Y de nuestro resultado de clasificación no es una variable de entrada), deberíamos esperar pasar una función a nuestro regresor logístico que tenga una forma similar a la siguiente: \[ g(\textrm{longitud}) = θ_0+ θ_1\textrm{longitud}\] A priori, aún no sabemos cuáles son los pesos. Lo que sí queremos es que sean tales que nuestra función \(g(x)\), al pasarla a nuestra función sigmoidea, nos dé una clasificación que se parezca razonablemente a lo que vemos en nuestros datos:

lengths <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
t1 = -4.5
t2 = 1
g = t1 + t2 * lengths
s = 1/(1 + e^-g)
data.frame(lengths, g, s)
##    lengths    g          s
## 1        1 -3.5 0.02931223
## 2        2 -2.5 0.07585818
## 3        3 -1.5 0.18242552
## 4        4 -0.5 0.37754067
## 5        5  0.5 0.62245933
## 6        6  1.5 0.81757448
## 7        7  2.5 0.92414182
## 8        8  3.5 0.97068777
## 9        9  4.5 0.98901306
## 10      10  5.5 0.99592986

Este fragmento de código toma las longitudes de los tumores de entrada, que van de \(1\) a \(10\), y elige dos pesos de \[\theta_0= -4,5\] y \[\theta_1=1\] En la práctica, habría que experimentar con la elección de los valores de los pesos y viendo cómo reaccionan los resultados, o a través de un algoritmo que le dé la respuesta. El código anterior proporciona la respuesta como resultado final. Se utilizan como pesos para la función \(g(x)\) que se pasa a la sigmoidea. La tabla del código presenta la clasificación resultante de los datos como \(s\). Un tumor de longitud \(1\), cuando se pasa por la función de entrada \(g(x)\) da un resultado de \(-3,5\). Este valor, cuando se pasa a la función sigmoidea, da un resultado que es bastante cercano a cero. Esto significa que un tumor de longitud \(1\) tiene una probabilidad muy baja de ser maligno, como se demuestra en la Figura 4-10:

plot(y = s, x =lengths, pch = 1, main = "Entradas de la función sigmoide y estimaciones de redondeo",
xlab = "Longitud del tumor", ylab = "Probabilidad de tipificación de clase 1")
points(y = round(s), x = lengths, pch = 4)
*Figura 4-10. Para una longitud de entrada dada, se puede ver la estimación de la función sigmoide en círculos, y su valor redondeado en cruces*

Figura 4-10. Para una longitud de entrada dada, se puede ver la estimación de la función sigmoide en círculos, y su valor redondeado en cruces

La Figura 4-10 presenta las probabilidades de que las longitudes de los tumores se clasifiquen como malignas si la probabilidad es de 1.0 y benigna si la probabilidad es de 0.0. El resultado es bastante cercano, pero hay algún error en él. Se obtendría una imagen mucho mejor si simplemente se redondean los valores al número entero más cercano.
El resultado final es una clasificación que se parece exactamente como los datos de partida.
Empezamos con los datos de entrada siendo la longitud del tumor. La salida del tipo de tumor entre benigno, \(y = 0\), y maligno, \(y = 1\), ya nos fue dada. El objetivo era diseñar un modelo que calculara la probabilidad de que un tumor fuera benigno o maligno en función de su longitud. Para ello, partimos de la ecuación \[g (x) = θ_0 + θ_1x\] y, a continuación, encontrar los pesos \(θ_0\) y \(θ_1\) , que ayudaron a sacar valores que, cuando se pasan a través de una función sigmoide, proporcionan valores que parecen correctos para lo que necesitábamos.
Lo que obtenemos al final es un límite de decisión en \(\textrm{longitud} = 4.5\); cualquier valor por encima se clasifica como 1, y cualquier valor por debajo se clasifica como 0. Los mecanismos por los que funcionan los algoritmos de clasificación como la regresión logística para determinar esas ponderaciones de modelización son algo similares en su alcance a cómo se calculan las ponderaciones de regresión lineal simple. se calculan los pesos de la regresión. Sin embargo, dado que el objetivo de este texto es ser de naturaleza introductoria, le remitiré al apéndice estadístico de la regresión lineal. La regresión logística y muchos otros algoritmos de aprendizaje automático funcionan de forma similar.pero adentrarse demasiado en el ámbito de la optimización de algoritmos puede resultar demasiado matemático y perderíamos el enfoque de los algoritmos. perderíamos el foco en la comprensión y aplicación del ecosistema del aprendizaje automático en su conjunto.

Clasificación binaria

Todo lo que hemos hecho hasta ahora en términos de clasificación ha sido sobre datos binarios: el tumor es maligno o benigno. La Figura 4-11 muestra otro ejemplo en el que determinamos las clases basándonos en la distribución de los datos:

plot(iris$Sepal.Length ~ iris$Sepal.Width, main = "Longitud del sépalo vs Ancho del sépalo en iris",
xlab = "Ancho del sépalo", ylab = "Longitud del sépalo")
*Figura 4-11. Puede utilizar la regresión logística para datos más dispersos en lugar de ser discretos en la salida*

Figura 4-11. Puede utilizar la regresión logística para datos más dispersos en lugar de ser discretos en la salida

La figura 4-11 Puede utilizar la regresión logística para los datos que están más repartidos en lugar de ser discretos en la salida. En la Figura 4-11, hay un montón de puntos de datos y lo que parece ser dos clases diferentes de plantas. Parece que hay una agrupación de puntos de datos en la parte inferior del gráfico que parecen estar más separados que los otros. Puedes ajustar un modelo de regresión logística a estos datos y encontrar la ecuación de la línea que hace de límite de decisión. Todos los puntos por debajo de ese umbral se clasificarán como un tipo, y todos los puntos por encima de la línea se clasificarán como de otro tipo. Este ejercicio utiliza un modelo lineal generalizado, dado por la función glm(). Su uso es más flexible que el de la función de modelo lineal estándar lm(), ya que se puede utilizar para propositos de clasificación.

iris.binary<- iris
iris.binary$binary<-as.numeric(iris[, 5]=="setosa")
iris.logistic<-glm(binary~Sepal.Width + Sepal.Length, data = iris.binary,family="binomial")
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
iris.logistic 
## 
## Call:  glm(formula = binary ~ Sepal.Width + Sepal.Length, family = "binomial", 
##     data = iris.binary)
## 
## Coefficients:
##  (Intercept)   Sepal.Width  Sepal.Length  
##        437.2         137.9        -163.4  
## 
## Degrees of Freedom: 149 Total (i.e. Null);  147 Residual
## Null Deviance:       191 
## Residual Deviance: 2.706e-08     AIC: 6

El resultado de este método proporciona algunos coeficientes e interceptos que no parecen totalmente correctos. Se necesita un paso adicional para calcular la pendiente y los interceptos de la límite de decisión de esta manera. Recuerda por un momento cómo utilizaste la función sigmoide \[g(z) = 1/(1 + e^{-z})\] \(z\) es una función con la siguiente forma \[z = θ_0 + θ_1x_1 + θ_2x_2\] Como se desea un valor entre 0 y 1 para la clasificación binaria, la clasificación es 1 cuando tienes tu función sigmoide \(g(z) ≥ 0,5\). Eso sólo es cierto cuando la función \(z\) que le pasas es a su vez mayor que 0: \[0 ≤ θ_0 + θ_{1}x_{1} + θ_{2}x_2\] Puedes reescribir esta ecuación y resolver el valor de nuestro \(x_2\) en consecuencia: \[x2 ≥\frac{−θ_0}{θ_2}+\frac{−θ_1}{θ_2}x1\] Esta ecuación tiene la misma forma que una recta \(y = b + mx\), donde podemos resolver computacionalmente la pendiente y el intercepto para construir la función que determina el límite de decisión límite: \[\text{pendiente} = - θ_1/θ_2 = - 137.9/-163.4=0.843\] \[\text{intercepto} = - b/θ_2 = - 437.2/-163.4=2.675\] Esto se puede calcular directamente desde el modelo logístico:

knitr::kable(head(iris))
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa 
slope.iris <- coef(iris.logistic)[2]/(-coef(iris.logistic)[3])
int.iris<- coef(iris.logistic)[1]/(-coef(iris.logistic)[3])
slope.iris
## Sepal.Width 
##   0.8440957
int.iris
## (Intercept) 
##    2.675511

A continuación, puede trazar esto sobre sus datos y ver cómo las clases se presentan, como se ilustra en la Figura 4-12:

iris.binary$binary[iris.binary$binary == 0] <- 2
plot(Sepal.Length ~ Sepal.Width, data = iris.binary, pch = (binary), main = "Longitud del sépalo vs Ancho del sépalo en iris", xlab = "Ancho del sépalo",
ylab = "Largo del sépalo")
abline(a = int.iris, b = slope.iris)
*Figura 4-12. Dividir una distribución de datos en dos clases*

Figura 4-12. Dividir una distribución de datos en dos clases

Ahora tiene una ecuación que ayuda a determinar cómo separar las especies de flores de iris que tenemos. Todas las flores del conjunto de datos que tengan un valor por debajo de la línea de \(y = 2.68 + 0.844 × \text{(Ancho del sépalo)}\) se clasificarán en consecuencia.

Clasificación multiclase

Si quiere encontrar las divisiones en sus datos que definen múltiples clases y no sólo una clasificación binaria, necesita utilizar la clasificación multiclase. Este enfoque es ligeramente diferente, ya que básicamente está aplicando el mismo esquema de clasificación binaria que ha estado haciendo hasta ahora, pero está comparando la clase que le interesa frente a todo lo demás.

La Figura 4-13 ilustra cómo se vería un ejercicio de clasificación multiclase:

multi <- data.frame(x1 = c(0.03, 0.24, 0.21, 0, 0, 0.23, 0.6, 0.64, 0.86, 0.77), x2 = c(0.07, 0.06, 0.19, 1.15, 0.95, 1,0.81, 0.64, 0.44, 0.74), lab = c(1, 1, 1, 2, 2, 2, 3, 3,3, 3))
plot(x2 ~ x1, pch = lab, cex = 2, data = multi,main = "Clasificación multiclase",xlab = "x", ylab = "y")
*Figura 4-13. Los datos multiclase se pueden separar solo mediante el uso de un enfoque de uno contra a muchos*

Figura 4-13. Los datos multiclase se pueden separar solo mediante el uso de un enfoque de uno contra a muchos

Hay tres clases distintas de datos, y usted quiere encontrar algún tipo de líneas que las dividan en sus propias categorías, de forma muy parecida a como lo hizo en el caso binario. En esencia, esto se reduce a nuestra simple prueba binaria, pero se cambia el grupo con el que se compara. Esto se llama una prueba “uno contra todos” o “uno contra muchos”, en la que se prueban tres casos: triángulos contra resto, círculos contra resto y cruces contra resto, como se muestra en la Figura 4-14:

par(mfrow = c(1, 3))
multi$lab2 <- c(1, 1, 1, 4, 4, 4, 4, 4, 4, 4) 
plot(x2 ~ x1, pch = lab2, cex = 2, data = multi,main = "Clasificación multiclase",xlab = "x", ylab = "y")
multi$lab3<- c(4, 4, 4, 2, 2, 2, 4, 4, 4, 4) 
plot(x2 ~x1, pch = lab3, cex = 2, data = multi,main = "Clasificación multiclase",xlab = "x", ylab = "y")
multi$lab4 <- c(4, 4, 4, 4, 4, 4, 3, 3, 3, 3) 
plot(x2 ~ x1, pch = lab4, cex = 2, data = multi,main="Clasificación multiclase",xlab = "x", ylab = "y")
*Figura 4-14. La clasificación uno vs muchos utiliza el mismo enfoque que la clasificación estándar, pero reclasifica otros grupos en un solo grupo para comparar*

Figura 4-14. La clasificación uno vs muchos utiliza el mismo enfoque que la clasificación estándar, pero reclasifica otros grupos en un solo grupo para comparar

En un enfoque de clasificación uno contra varios, se utiliza un límite de decisión para clasificar los datos de un tipo o clase frente a todos los demás tipos o clases. A continuación, se hace lo mismo con el resto de los tipos o clases de los datos hasta que se dispone de una serie de límites de decisión que se pueden utilizar para tipificar los datos en consecuencia. Así, para el ejemplo de la Figura 4-14, se calculan (en la gráfica de la izquierda) los círculos frente al resto de los datos, y luego se calculan los triángulos frente al resto de los datos (gráfica del medio) y, finalmente, las cruces frente al resto de los datos(gráfica de la derecha). Al dividir un problema de tres clases en tres problemas de dos clases, puede encontrar más fácilmente un único límite de decisión para cada gráfica y luego combinar esos límites de decisión para un modelo final. En este caso, se recurre a la función multinom() de la biblioteca nnet. Se utiliza para pasar un caso multinomial que es básicamente lo mismo que se ha hecho para el caso binario simple, pero con tres valores en lugar de dos. Esta metodología puede aplicarse para más de tres categorías:

require(nnet) 
multi.model<-multinom(lab~x2 + x1, data = multi, trace = F)
multi.model
## Call:
## multinom(formula = lab ~ x2 + x1, data = multi, trace = F)
## 
## Coefficients:
##   (Intercept)       x2        x1
## 2   -12.47452 28.50805 -17.97523
## 3   -19.82927 12.95949  33.39610
## 
## Residual Deviance: 0.0004050319 
## AIC: 12.00041

De nuevo, sin embargo, necesita hacer el cálculo especial para las pendientes e interceptos de los límites de decisión basados en la salida de este modelo. Puede aplicar la misma matemática que antes, pero tiene que aplicarla a cada una de las ecuaciones del modelo. A continuación, puede trazar las líneas límite de decisión, como se ilustra en la Figura 4-15:

multi.int.1<- coef(multi.model)[1]/coef(multi.model)[3] 
multi.slope.1<- -coef(multi.model)[5]/coef(multi.model)[3]
multi.int.2 <- -coef(multi.model)[2]/coef(multi.model)[4] 
multi.slope.2 <- -coef(multi.model)[6]/coef(multi.model)[4]
plot(x2 ~ x1, pch = lab, cex = 2, data = multi,main = "Clasificación multiclase",xlab = "x", ylab = "y")
abline(multi.int.1, multi.slope.1) 
abline(multi.int.2, multi.slope.2)
*Figura 4-15. Dos líneas que separan tres clases de datos*

Figura 4-15. Dos líneas que separan tres clases de datos

Regresión logística con caret

El paquete caret hace que los problemas de regresión logística sean muy sencillos para ejemplos más complejos que los que hemos estado haciendo hasta ahora. El uso de caret es bastante sencillo, aunque para algunos métodos de aprendizaje automático en particular, pueden estar justificadas algunas optimizaciones y ajustes para lograr los mejores resultados posibles. A continuación se muestra un ejemplo de cómo se puede realizar una operación con caret:

require(caret)
data("GermanCredit")
Train<-createDataPartition(GermanCredit$Class, p = 0.6, list = FALSE) 
training<- GermanCredit[Train, ]
testing <- GermanCredit[-Train, ]
mod_fit<- train(Class ~ Age + ForeignWorker+Property.RealEstate + Housing.Own + CreditHistory.Critical, data = training, method = "glm",family = "binomial")
predictions<-predict(mod_fit, testing[, -10]) 
knitr::kable(table(predictions, testing[, 10]))
Bad Good
Bad 19 31
Good 101 249

Este sencillo ejemplo utiliza datos del conjunto de datos GermanCredit y muestra cómo se puede puede construir una matriz de confusión a partir de un objeto de entrenamiento caret. En este caso, el ajuste no parece no parece muy bueno, porque alrededor del 50% de los datos parecen estar clasificados incorrectamente. Aunque caret ofrece algunas formas estupendas de afinar cualquier método de aprendizaje automático método de aprendizaje automático en el que esté interesado, también es bastante flexible a la hora de cambiar los métodos de aprendizaje automático. Simplemente editando la opción de método, puede especificar uno de los otros 12 algoritmos de regresión logística para pasar al modelo, como se muestra aquí:

mod_fit<- train(Class ~ Age + ForeignWorker + Property.RealEstate +Housing.Own + CreditHistory.Critical, data = training,
method = "LogitBoost",
family = "binomial")
predictions <- predict(mod_fit, testing[, -10])
knitr::kable(table(predictions, testing[, 10]))
Bad Good
Bad 18 27
Good 102 253